加载页面中...
毕设研究思路及结果图 | lwstkhyl

毕设研究思路及结果图

other-other

研究思路概括和写进正文中的结果图，以及相应结论（byGPT）

构建Seurat对象和数据预处理

从STARsolo/stellarscope/scTE结果的计数矩阵构建Seurat对象

  • 比对参数均遵循官方说明文档

质控前：共243200个细胞，78691个基因

  • 细胞阈值：

    dataset min_features max_features min_counts max_counts max_mito max_ribo
    GSE157827 200 8000 500 40000 14 3
    GSE174367 300 12000 1000 50000 10 2
  • 基因阈值：选择至少在0.05%的细胞中有表达的基因
  • scDblFinder去双细胞

质控后：共207009个细胞，35688个基因

  • 14968个hERV位点
  • 58个hERV亚家族
  • 3个hERV超家族
> table(seu_filt$dataset)
GSE157827 GSE174367 
   151686     55323 
> table(seu_filt$group)
    NC     AD 
 90772 116237

数据预处理：

  • Seurat官方的经典流程：NormalizeData+FindVariableFeatures+ScaleData+RunPCA
  • 使用harmony去批次效应
  • 选前30个pc
  • 分辨率设置为0.1

01qc

  • A：QC阈值分布图（过滤前），GSE174367这个数据集测序深度较深
  • B：各样本QC前后细胞数对比，GSE174367这个数据集测序质量较好
  • C：去批次效应后UMAP图按数据集分组，展示两个数据集经过去批次效应后整合程度较高

02herv_counts

  • A：stellarscope计数结果，展示hERV在单细胞中整体占多少比例，两个dataset之间大概是什么水平，重点是HERV_fraction，该值在0.2%-0.3%左右，符合stellarscope官方的预测比例
  • B/C：scTE计数结果，展示各超家族和各亚家族的计数占比，ERVL和ERV1最多，HERVH最多

细胞类型注释：marker基因参照这两个数据集的原论文，共分为7个细胞类型

  • Astro：AQP4/GFAP/SLC1A2/GPC5/ADGRV1
  • Endo：CLDN5/FLT1/EBF1
  • Excit：SLC17A7/CAMK2A/SATB2/SNAP25/CBLN2
  • Inhib：GAD1/GAD2/LHFPL3/PCDH15
  • Mic：P2RY12/CSF1R/CX3CR1/CD74/C3
  • Oligo：MBP/PLP1/MOBP/OPALIN/MOG
  • OPC：PDGFRA/CSPG4/VCAN

03celltype

  • A：UMAP图按细胞类型分组
  • B：各细胞类型的细胞数占比
  • C：marker基因在各细胞类型中的表达量

超家族/亚家族层面的hERV重调情况

使用Seurat包的FindMarkers函数，test.use使用MAST方法，min.pct设置为0.05（至少在5%的细胞中表达）

04HERV_subfamily_scTE_celltype_avgexp_heatmap

  • 各细胞类型的hERV表达量和差异表达分析的显著性：按行标准化，每个格子代表该hERV家族在各细胞类型及AD-NC中的相对表达量
  • 可以看到hERV家族的表达量在各细胞类型中以及AD-NC中均有差异

05class_family_DE_summary

  • A：class层级的bubble plot，HERV-total指HERV的scTE计数总量，颜色是avg_log2FC，点大小是-log10(FDR)。对应结论中“herv全局有变化，但class层级幅度很小，细胞类型之间不一致”
  • B：部分亚家族的bubble plot，在总的DE结果中选出差异较显著、表达量较高的有研究价值的亚家族。对应结论中“family-level变化不是统一方向的，不同celltype的模式不同，能看到一些重复出现的family，但幅度仍不算大”
  • C：effect size-significance散点图，每个点是一个family-celltype组合，横轴是avg_log2FC，纵轴是-log10(FDR)，按细胞类型分面。表示“很多点虽然统计显著，但横向偏移并不大，说明这是大样本下的小效应，而不是强烈的family-level重调”，直观表明“p值很小但意义有限”
  • D：部分亚家族的donor-level箱线图，每个donor在某celltype中的平均family表达量，说明“组间趋势有时存在，但donor间离散度很大，所以family-level差异很容易被个体差异覆盖”

在family/class层级，AD相关hERV改变并非完全缺失。无论在全体细胞还是部分细胞类型中，均可检测到若干family/class达到统计学显著，说明AD与hERV转录改变之间仍存在一定关联。然而，这些差异的效应值整体较小，多数表现为小幅偏移而非强烈重调。进一步分细胞类型分析后，可以看到family-level变化具有明显的细胞类型依赖性：Excit、Inhib、Oligo和Mic中可见相对更清晰的变化，而Endo和OPC中信号较弱。与此同时，同一class或family在不同细胞类型中的方向并不完全一致，提示将多个位点聚合到family/class层级后，位点间异质性会被平均，从而削弱真正具有生物学意义的位点特异信号

  • 全局层面不是“强重编程”，而是“很弱的小幅偏移”：ERVK、ERVL、HERV-total在全体细胞里都能显著，但avg_log2FC都很小。在这么大的细胞数下，小偏移也会有很小的p值，所以这里更该强调效应值小，而不是只看显著性
  • family层级的变化明显是细胞类型依赖的：Endo最弱、OPC也较弱、Mic/Oligo/Excit/Inhib更能看到变化，把所有细胞混在一起时，family/class层面的差异会被不同细胞类型中方向不一致的变化互相抵消，只有分细胞类型后，才能看到较清晰但仍偏小的family-level效应
  • family聚合会掩盖位点异质性：HERVK-int在All/Astro/Excit里偏下调，但更粗的ERVKclass层级，在某些细胞类型里又是轻度上调或接近不变，说明同属一个family/class的不同位点，并不一定同向变化，因此family聚合后，真实存在的位点特异信号会被平均掉
  • 结合后续位点级DE、WGCNA及hERV位点-gene相关性结果，说明更有解释力的AD相关hERV信号主要存在于特定位点及特定细胞状态，而非整个family的一致性上调或下调。

WGCNA

对表达比例≥10%的基因进行hdWGCNA分析，在每个细胞类型中分别构建出5~10个基因模块

亚家族

06wgcna_mic

  • A：每个cluster内的hERV家族表达量
  • B：每个cluster内的基因模块表达量
  • C：基因模块与hERV家族表达量的相关性
  • D：根据ABC图得到的结果，以及各cluster的细胞数量和AD-NC比例，选择细胞数较多（约为1k，AD细胞占比较大）的cluster2作为主要分析亚群，又根据对各基因模块的GO富集分析结果，选择blue和pink模块作为主要分析模块。将cluster2与其它亚群进行差异表达分析，得到差异显著（p值<0.05，|logFC|>0.1）的基因，将DEG与各模块基因取交集后进行GO富集分析

在小胶质细胞中，hERV高表达状态主要与广义先天免疫/炎症样激活相关，同时伴随稳态和神经互作相关程序下降。具体来看，Mic的cluster2中AD比例相对升高，并表现出较明显的hERV相关免疫激活特征。与hERV正相关的pink模块在该亚群中上调，其功能富集集中在defense response to virus、response to virus和defense response to symbiont等条目，同时相关差异基因包括TLR2、MX1和MX2，提示这一亚群更接近先天免疫和抗病毒样激活状态。相反，与hERV负相关的blue模块在该亚群中下调，其GO条目涉及learning、cognition、trans-synaptic signaling和neuron projection morphogenesis，说明hERV高表达可能伴随小胶质细胞稳态或神经互作相关程序的减弱。总体来看，Mic中的结果最能体现“具体hERV家族—特定细胞类型—功能状态”的对应关系，即HERVH/HERVL/HERVK9主要关联小胶质细胞的先天免疫、炎症激活和抗病毒样反应

（这几个补充图不会出现在正文中，只会作为附录，因此这段结果不需提到cluster层面）

补充01wgcna_astro

  • 在星形胶质细胞中，hERV相关模块主要富集到gliogenesis、cell adhesion、cell-cell signaling和synapse maturation等方向，说明hERV相关状态也可能参与星形胶质细胞的神经支持、细胞通讯和轻度炎症响应

补充02wgcna_excit

  • 兴奋性神经元提供了与胶质细胞不同的模式，hERV相关模块主要富集到synaptic signaling、synaptic vesicle cycle、neurotransmitter release和vesicle-mediated transport in synapse等条目，而这些突触功能模块整体与hERV呈负相关，提示在兴奋性神经元中，hERV升高更可能对应突触传递和神经元功能程序的下降，而不是类似胶质细胞中的免疫激活或重塑状态

补充03wgcna_oligo

  • 在少突胶质细胞中，不同hERV家族与少突状态的关系并不完全一致。HERVK3、HERVK11与高代谢和高翻译活性状态存在关联；而HERVH、HERVL和HERVK等多个家族则更可能共同参与膜脂代谢、膜转运和信号重塑相关状态以及结构重塑和黏附变化趋势。这几个细胞类型的GO证据相对较弱，需要进一步验证和分析

OPC、Inhib和Endo：OPC中存在protein folding/stress、synaptic projection和STAT1/external stimulus相关线索；Inhib中主要表现为较弱的protein folding/stress和神经突触相关信号；Endo中更多体现vascular/circulation和DNA replication/cell-cycle相关背景变化。这些结果整体相关性和GO富集结果较弱，只作为简单补充

07wgcna_summary

hERV重调在AD脑内具有明显的细胞类型和细胞状态依赖性。当前结果中，最主要的发现集中在小胶质细胞：以HERVH、HERVL和HERVK9为代表的hERV家族与小胶质细胞的先天免疫/炎症样激活、抗病毒样反应和免疫重塑相关，同时与稳态和神经互作相关程序下降相联系。在少突胶质细胞中，HERVK3、HERVK11选择性关联高代谢和高翻译活性状态，而HERVH/HERVL/HERVK相关信号更多表现为结构重塑或膜脂/信号重塑趋势。在兴奋性神经元中，hERV升高更可能与突触传递程序减弱相关。在星形胶质细胞中，hERV主要关联胶质支持、神经互作和轻度炎症响应。整体来看，hERV关联不是统一的，而是明显细胞类型依赖，而且AD和hERV不是简单同向变化，在很多模块里都能看到hERV和模块相关明显，但与AD-NC状态的相关性较弱，甚至方向相反，说明hERV更像是反映细胞状态轴的标志，而不是简单的AD/NC二分类标志

位点

对hERV位点进行差异表达分析，使用的方法同亚家族的差异表达分析，在得到的结果中选出在各细胞类型中较显著且表达比例较高的位点，与前面得到的WGCNA基因模块进行相关性分析

  • 相关性分析的方法——线性拟合：先计算每个donor在对应细胞类型内的hERV平均表达和基因集平均表达，然后用线性模型校正log10_nCount_RNA和percent_mito的影响，提取残差后再进行Spearman相关性检验

在与hERV位点相关性较大的模块中选hub gene，然后算样本水平上的位点~gene相关性，确定“该module里哪些具体基因和这个位点关系最稳定”

  • 每个“细胞类型-位点”保留相关性绝对值最大的前2个module
  • 提取这些module的hub gene
  • 在对应celltype内、对应module里，找与该位点最相关的具体基因
  • 找“既是hub gene，又和该位点显著相关”的交集基因，作为样本级相关性的输入

08locus_DE_WGCNA_geneCorr

  • A：选出来的位点在各细胞类型中的差异表达分析结果，点大小是表达比例（AD-NC中取最大值），颜色是avg_log2FC，星号说明显著性。展示这些位点在DE里为什么值得看
  • B：选取的与位点作相关性分析的基因在哪个模块里，就画位点和哪个模块的相关性。展示这些位点在WGCNA里主要连到哪些module
  • C：位点-基因样本级相关性图。展示这些module最后落到了哪些具体gene，并且donor-level支持

位点的差异表达分析结论：

  • excit/inhib神经元内：HERV3-1p34.3这个位点logFC很大，细胞比例也是AD更多
  • 多个细胞类型一致下调、效应量不小：
    • HML3-1q32.1b
    • HML1-1q32.1
    • HML4-16p13.3
  • Oligo内特异下调的位点比较多，不是零散几个位点变化，而是一组位点整体偏低，更像是一个系统性状态改变
  • Mic和Inhib有一个上调的HARLEQUIN-6p21.31，而且它在所有细胞水平也上调，像是一个跨细胞类型但偏向非少突的上调位点

位点wgcna结论：在此基础上，我们进一步结合位点级WGCNA中显著module的key genes，并要求cell-level与donor-level相关方向一致、样本层面能够稳定复现，最终筛选出更适合展示和解释的位点-gene组合。结果显示，真正稳定的位点-gene关联主要集中在Mic、Oligo、Astro和Excit四类细胞中

  • Mic中，HARLEQUIN-6p21.31与IL4R、ADGRE2、ALOX5和HCK等免疫相关基因在cell-level和donor-level均呈一致正相关，提示该位点更可能标记一条microglial免疫活化状态轴
  • Oligo中，HERVW-15q21.2与GLDN和PIK3C2B呈正相关，提示该位点与少突细胞特定状态或结构程序相关；MER41-13q13.3与ATP11A呈负相关，则提示Oligo中可能还存在一条与膜状态调节相关的反向位点-gene轴
  • Astro中，HML1-1q32.1和HML3-1q32.1b分别与ABCC9和COLEC12相关程序呈稳定正相关，其中ABCC9在两个位点上均获得支持，提示这两个HML位点并非彼此独立，而更可能属于同一条共享的astrocytic hERV状态轴
  • Excit中，HML4-16p13.3与MTR和MICU2呈稳定正相关，提示该位点可能关联神经元的一碳代谢及线粒体功能/钙稳态背景
  • 相比之下，Inhib、OPC和Endo中虽可见一定cell-level位点-gene关联，但多数未在donor层面稳定复现，提示这些关联更可能反映细胞内局部异质性，而非稳健的样本级协同变化
  • 总体来看，位点级hERV并非简单重复AD/NC差异，而是在不同细胞类型中对应若干相对独立的转录状态轴

基因组浏览器和GWAS分析

使用pyGenomeTracks画基因组上hERV位点位置，使用的gtf注释只保留编码基因

09genome_pos_noNCRNA

10genome_gene

  • A：根据该结果，又选了几个与目标hERV位点距离较近的基因作相关性分析
  • B：上述提到的各基因的DE结果汇总，点大小是表达比例（在该细胞类型中的所有细胞中），颜色是avg_log2FC，星号说明显著性

基因组浏览器和基因相关性结论：结合genome browser局部注释、局部候选基因共表达分析以及前文位点级WGCNA结果，我们将代表性locus-level hERV位点大致分为三类

  • 第一类为宿主基因相关型，包括HARLEQUIN-6p21.31和HERVW-15q21.2：前者位于ANKS1A基因内部并与ANKS1A呈显著正相关，后者位于GLDN内部且与GLDN稳定正相关，提示这两个位点的表达信号更可能与局部宿主转录本或特定isoform相关
  • 第二类为共享非编码转录单元型，包括HML1-1q32.1和HML3-1q32.1b：二者位于同一局部区域，但均未获得PM20D1或SLC26A9的明确宿主支持，更可能共同代表一条astrocytic非编码hERV位点轴
  • 第三类为独立转录/调控相关型，包括MER41-13q13.3和HML4-16p13.3：二者均缺乏典型宿主基因同链结构，其中MER41-13q13.3与NHLRC3呈区域性正相关，而HML4-16p13.3与KCTD5呈显著负相关，提示其更可能反映独立hERV转录，或与局部调控环境相关的表达信号
  • 根据参考文献《Integrating human endogenous retroviruses into transcriptome-wide association studies highlights novel risk factors for major psychiatric conditions》位于宿主基因内部且同链的位点更容易被解释为宿主转录本相关信号，而反义链或基因间区位点则更可能对应独立的ncRNA样转录；同时，短读长数据只能提供支持性证据，尚不足以完全确定精确转录本结构

之后在GWAS和ADVP数据库中查询以上基因

基因DE与ADVP/GWAS结论：位点级hERV-gene关系具有明显的细胞类型特异性，其中Oligo和Astro的证据链最完整，既能在局部基因组环境中获得较清晰的结构解释，也能在GWAS/ADVP中获得较好的神经认知、脑结构、衰老或痴呆背景支持；Excit具有一定的代谢和脑结构相关支持，但强度略弱；Mic则更多反映免疫程序相关的状态轴，而非直接的AD特异遗传命中。整体上，这些结果更支持位点级hERV所对应的是若干细胞类型特异的转录状态轴，这些状态轴与神经认知、脑结构及衰老遗传背景存在不同程度联系，而不是简单重复经典AD风险基因或单一的AD/NC均值差异

  • 这些与位点级hERV稳定相关的候选基因，并不普遍对应经典AD核心风险基因。总体上，ADVP中的直接命中相对有限，而GWAS Catalog中则能看到更多与认知、脑结构、衰老及痴呆相关的条目，说明这些候选基因并非随机筛出，而是落在更广义的神经生物学遗传背景中
    • GLDN在browser中属于典型宿主基因内同链位点，且与HERVW-15q21.2稳定正相关；在GWAS侧虽然不是最直接的AD命中，但可见cognitive performance、educational attainment、intelligence、memory test以及neurofibrillary tangles等相关条目，提示HERVW-15q21.2所对应的少突状态轴，可能连接到更广义的脑功能和认知遗传背景
    • ABCC9虽然不是经典AD核心基因，但在GWAS中可见hippocampal sclerosis of aging和entorhinal cortical thickness等与脑老化、神经退行背景贴近的表型，同时又与astrocytic hERV轴稳定相关，因此更适合被解释为“脑老化背景中的astrocytic状态相关基因”
    • COLEC12则更强一些，因为它既与HML3相关程序稳定耦联，又可见dementia直接条目及cortical surface area等脑结构相关记录，说明这条astrocytic hERV轴并不是简单重复疾病均值差，而是落在更广义的神经退行与脑结构背景中
    • NHLRC3对MER41-13q13.3这条Oligo位点尤其有价值：MER41在browser上本来就更偏独立/局部区域耦联型，而NHLRC3又可见late-onset Alzheimer’s disease相关条目，因此可作为“区域层面的AD相关支持”
    • MTR具有hippocampal volume相关记录，可为Excit中的代谢/神经元功能解释提供中等强度支持
    • PIK3C2B虽有Alzheimer’s disease or educational attainment (pleiotropy)等线索，但更适合放在“认知/状态协同背景”而非直接AD输出的框架中
    • MICU2、KCTD5、ADGRE2和IL4R则分别更偏泛脑结构量化表型、局部调控环境补充支持或精神/免疫背景支持
  • 需要特别说明的是，部分基因的DE方向与预期相反（图B标红框的基因），例如Oligo中HERVW-15q21.2在AD中下调，但其正相关基因PIK3C2B却在AD中上调；类似地，COLEC12和MICU2也存在与对应hERV位点整体方向不完全一致的情况。这并不否定位点-gene相关性本身，而是说明“相关性”和“组间差异”描述的是两种不同层面的信号。前者刻画的是细胞类型内部或不同样本之间的协同变化轴，即某个位点升高时，哪些基因倾向于一起升高或降低；后者刻画的是AD与NC之间的组间均值差异。二者可以同向，也可以不同向。换句话说，某个基因往往同时受到两股力量影响：一股是与hERV绑定的局部状态轴，另一股是AD背景下相对独立的疾病效应。当后一股力量更强，或者当AD组中某类亚状态/细胞组成比例发生变化时，就可能出现“hERV与gene在组内稳定正相关，但gene在AD/NC均值比较中却呈相反方向”的现象。结合结果来看，这些方向不一致的基因通常avg_log2FC都比较小，更像是大样本检出了一个幅度有限的组间偏移，而不是特别强的生物学输出，说明位点级hERV更可能对应细胞状态协同变化轴，而非简单等同于疾病组间均值变化

选几个最关键的基因-位点pair看看在其它细胞类型中的情况

11supply_gene_herv_DE_corr

  • 跨细胞类型比较显示，不同代表性位点-gene pair具有不同的泛化模式
    • GLDN ~ HERVW-15q21.2在Oligo中最强，但在Excit、Mic及Astro中仍可见一定一致性，提示其并非严格局限于单一细胞类型
    • ABCC9 ~ HML1-1q32.1则明显以Astro为主，离开Astro后相关性与DE一致性显著减弱，支持其更偏astrocytic状态轴
    • MTR ~ HML4-16p13.3除Excit外，在Oligo、Astro和Inhib中亦保留一定一致模式，提示该位点可能对应更广义的代谢/线粒体功能背景
    • 总体上，这些结果进一步说明位点级hERV-gene关系并不是统一的泛细胞类型规则，而是由“主细胞类型特异性”与“有限跨细胞类型延伸”共同构成的

结合WGCNA结果可以看到，这些位点-基因pair与对应细胞类型的功能模块具有一定一致性。GLDN-HERVW-15q21.2主要出现在Oligo中，而Oligo的hERV相关模块本身就涉及膜状态、结构重塑和细胞功能调节；ABCC9-HML1/HML3主要出现在Astro中，与Astro中胶质支持、细胞黏附和细胞通讯相关模块相吻合；MTR-HML4在Excit中较稳定，而Excit的hERV相关模块与突触传递和神经元功能下降有关，MTR所代表的代谢背景可能参与这种神经元状态变化。总体而言，位点级hERV-gene关系既有主细胞类型特异性，也有有限的跨细胞类型延伸。这说明AD脑组织中的hERV信号更可能对应若干细胞状态轴，而不是一种适用于所有细胞类型的统一规则。

局限性

样本数较少，很难用伪bulk的方法，只能用细胞级分析，容易出现p值虚低的情况，也很难排除单个样本特异性对结果的影响

使用的数据的建库方法是10x Chromium Single Cell 3′ v3——“用poly(dT)引物捕获带poly(A)尾的转录本”，可能会低估或者漏掉不带polyA的ncRNA（非编码RNA），同时TE里面有相当一部分是non-polyA的，所以在hERV计数上大概率会偏低。不过我这个分析主要是AD-NC/不同细胞类型之间的组间比较，或者是共表达这种定性的找趋势，所以在这类不需要精确计数的分析上应该影响不大。当然也不可能说完全没有影响：

  • 之前别人的bulk RNA-seq分析用的是total RNA的测序方法（而不是我这种polyA+的测序方法），这可能是结论有差异的原因之一
  • 我这里能被polyA+的测序方法测出来的hERV读段，可能是hERV中有一定特殊性的片段，比如是剪接到了正常gene中的、或是借用了周围的polyA