Skip to content
#转录组

转录组数据分析内容

约 4947 个字 预计阅读时间 16 分钟

专题1——转录组分析探讨

1.1 转录组研究背景

转录组学是一门在整体水平上研究细胞中基因转录的情况及转录调控规律的学科;借助mRNA测序(mRNA-Sequencing)使迅速成为研究疾病状态、生物过程等各种研究转录组变化的首选方法。除了作为一种高度敏感和准确的基因表达量化手段，mRNA测序还可以识别已知和新的转录本亚型、基因融合和其他特征以及等位基因特异性表达；mRNA测序可以提供一个不受经验限制的编码基因转录组完整视图。

1.2 转录组送样/建库/测序

  • 样本重复数：

  • 对于细胞等样品间差异较小的样品，建议生物学重复≥3个；

  • 对于动植物等样品间具有一定差异的样品，建议生物学重复≥5个；

  • 对于个体差异较大的临床样品，建议生物学重复≥10个；

  • 理论上生物学重复越多统计学意义越大，根据具体的实验设计合理安排，但也不是越多越好。

  • 送样量：

  • 样品质量

    一般建议样品的质量达到以下要求，特殊样品不参考。total RNA≥500ng，浓度＞50ng/μL，RIN≥7 ，28S/18S：1.8~2.0

  • 测序

    6G，PE150，优先推荐Illumina平台，当然国产平台也推荐，但就担心期刊的编辑可能会挑刺。

1.3 常规数据分析

1. 测序数据质控

  • 质控目的

    为了降低数据噪音，保证下游分析结果的准确性和可靠性，最直接的影响——基因组的比对率

  • 一般做哪些质控处理

  • 什么样的数据算是质控合格的

  • 原始下机数据的base有多少G?

    ==⇒ 推荐大于6G

  • 原始下机数据的Q20/Q30的比例是多少？

    RNA类别的项目推荐Q30>80%，DNA类型项目推荐Q30大于90%(Illumina平台)，ION/PGM等要求可以略低

  • clean有效的数据的占比是多少？

    推荐大于80%，比例过低，说明测序有问题，有可能是文库降解或者测序不稳定或文库不平衡等系列原因

  • rRNA占比多少？该不该去除？

    一般rRNA（核糖体RNA）通常不具有polyA结构，因此常规的polyA富集的建库方法中，rRNA理论上是不会出现的，一般推荐的他的占比<10%；针对提取total RNA建库的方案，有可能存在rRNA降解酶效率问题，rRNA可能占比很高，此时需要去除rRNA。

  • 工具选择(常用)

    software 特点 推荐度
    fastp C语言编写，分析速度快，可以自动识别adapter，自定义各项质控参数，有可视化报告，适用于各种测序平台技术和文库类型 **​
    trim galore 将cutadapt和FastQC进行了封装，可以实现和fastp类似分析，速度上较fastp慢，也可以用于各种测序平台技术和文库类型 ​
    Trimmomatic 基于Java编写，速度较前两者慢，需提供adapter序列，普遍用于Illumina的数据 -

2. 参考基因组比对

  • 目的

    为了确定这些片段由哪些基因/转录本转录而来的产物，即为了获取每一个片段(reads)的基因组坐标，进一步才能确定每一个转录本/基因的reads count(表达量)，结果往往会与上述的QC息息相关。

  • 比对结果质控指标

    理想情况下，基因组整体比对率(total_map)应该大于70%，且绝大部分reads是唯一比对率(unique_map>80%)，如果一个组中，部分样本比对率太低，在考虑重复数足够的情况下，可以酌情去除部分样本

比对率不达标？ 检查是否存在污染(取样污染，建库污染，测序污染)？参考基因组是否完整？样本类型是什么(eg：临床病人样本易出现支原体等污染)？样本做了什么处理？

  • reads基因组区域分布，exon的占比一般要大于90%，但是微量转录组70%就可以，如果内含子的占比过高，可能是建库过程中rRNA没有剔除干净，需要尝试trim rRNA

  • 其他指标

    基因组随机性分析/测序饱和度/基因body覆盖分布情况等，一般对于成熟的转录组产品，推荐的数据量都是饱和的，基因body是取决于3`端富集还是全长扩增；对于新产品和新技术，这些指标就很有必要进行分析探究。

  • 参考基因组怎么选？

    • 优先选择Ensembl/NCBI/gencode等社区比较大且共识度广的数据库，优先选者primary版本的参考基因组；但如果研究的参考基因组不在这些数据库中，推荐选择一些专有的数据库；

    • hg38?hg19/mm9/mm10等系列版本问题：优先选者较新(但不是最新)且版本稳定的版本，human目前优先hg38，小鼠优先mm10。

  • 工具选择(常用)

    software 特点 推荐度
    hisat2 C++语言编写，使用分层索引（hierarchical indexing）和全局Ferragina-Manzini (FM)索引结合多个局部FM索引，相比于Bowtie2和BWA等比对软件具有高效/低内存消耗/灵活性，且支持剪切比对 **​
    STAR C语言编写，专门针对RNA-seq数据比对工具，同样具有高效/灵活性，其对于剪切的感知和灵敏度优于hisat2，因此其可分析方向比hisat2广泛，可以用于基因融合/突变分析/circRNA分析等，此外自带定量功能，缺点是需要的内存比较大 ​
    TopHat2 早期(2010左右)转录组比对的最常用软件，后续被更快/更准确的hisat2和STAR两个软件慢慢取代

3. 基因表达定量

  • 目的

    确定每一个转录本/基因在每一个样本中的reads count(表达量)数

  • 工具选择(常用)

    software 特点 推荐度
    featureCounts featureCounts属于Subread软件中的定量工具，是一款使用于RNA-seq和DNA-seq的read summarization工具，应用了高效率的染色体哈希算法和feature区块技术。它比目前存在的工具速度都快，而且需要的内存空间少，同时可以用于单端和双端的数据 **​
    Salmon/kallisto 该两个软件并没有直接reads进行基因组比对 ，具有高效的算法、低计算成本和适用于大规模数据集。Kallisto 通过采用估计碎片相对丰度的方法，不需对整个转录组进行比对，加速了分析过程，且保证 RNA-seq 数据定量化分析与传统方法相似的精确性的基础上，极大地降低了运行时间；Salmon 采用概率模型，避免了传统比对方法的计算瓶颈，其优点包括高效的碎片量化、低计算成本和适用于大规模测序数据； ​
    HTSeq python编写，其原理于featureCounts类似，但速度奇慢 -
    其他 RSEM(使用无参)，StringTie(转录本组装)，BEDTools, Qualimap, GenomicRanges等

4. 重复性评估

  • 目的

    确认同一个处理组中，样本的重复性如何，酌情剔除样本重复不好的样本

  • PCA分析(Count矩阵）

PCA

  • 聚类热图(Count矩阵）

heatmap

  • 表达量分布(相对表达量）

expression

  • 为什么组内重复性不好？

    • Encode计划分析结果中建议皮尔逊相关系数的平方(R2 )大于0.92，基于经验和大量文献研究显示，实际的项目中要求生物学重复样品间皮尔逊相似性R2至少要均大于0.8；

    • 如果PCA图中组间样本汇合在一起(全部或部分样本混合)，且相关性系数不达标时，说明样本的组间重复性不好，此时应考虑：样本的类型，例如临床样本的高异质性的确会导致组内重复性很差，样本是否标记错误，组内样本实验处理是否一致(尤其是实验组样本)，是否组内数据存在多个批次合在一起；

    • 排除了上述原因，在重复数允许的情况下，可以酌情剔除组间混合的样本进行后续的分析

  • 工具选择(常用)

    R packages(eg：stats/Complexheatmap，Cor function等)

5. 差异分析

  • 目的

    筛选出显著的改变的基因，当然也可以筛选一致性变化或符合趋势的基因

  • 差异分析算法&阈值

其他算法：

  • limma： 适用于芯片数据或标准化后的表达谱数据(TPM、FPKM，RPKM等)，记住需要将表达谱进行log2转化；

  • t-test： 适用于两组数据定量、正态、独立、方差齐的数据，一般普遍用于于蛋白和代谢等数据，但是目前研究发现t-test的效能并不是最好的；

  • 方差分析(anova)： 适用于比较三个或三个以上样本均值是否存在显著差异的一种方法。ANOVA假设所有组的总体分布都是正态的，并且具有相同的方差。通过比较组内差异和组间差异，ANOVA可以确定不同组之间是否存在显著的均值差异。

  • 其他检验方法/工具

    Differential gene expression analysis pipelines and bioinformatic tools for the identification of specific biomarkers: A review

  • 差异基因数多少合适？

    差异基因数没有标准范围，其于实验的样本类型，实验处理等是相关的，但一般推荐差异基因数控制在2k以内是可取的，根据经验来看，重复数足够且组内重复性很好时，差异基因数也就几百个，此外基因数太多会给后续的分析带来噪音；重复性不好时，差异基因数会更多；

    差异基因数过多？可以适当将pvalue/或者padj卡严格一些；差异基因数过少？保证pvalue<0.05的前提下，适当放宽log2FC；

    差异基因数过少？特殊项目或者情况下，差异基因数会非常少，尤其是蛋白/代谢等项目，此时需要结果实验的重复性等来找原因，在条件允许的情况下，我们可以先放宽对pvalue和FC的限制(尤其是FC的限制)，尝试进行差异基因筛选，以差异基因的功能来协助验证是否可行。

    关注的基因pvalue>0.05？综合考虑各方面的因素：样本量是否足够？对于某些临床样本，当样本量较少时，检验效能不足才没有发现它们之间的真实差异；此时可以尝试计算本次研究的检验效能(1-β)，如果他们间确有差异，则本次能够发现差异的检验效能为100%，如果发现检验效能较小(如<0.5)，则很有可能是样本量不足导致的P>0.05，此时可加大样本量再继续尝试。其次还需要综合考虑样本类型，组间重复性等因素，必要时可以做一个qPCR进行验证。

  • 可视化怎么搞？

  • 工具选择(常用)

    R packages(eg：edgeR/DESeq2/limma/ggplot2/ggcoverage/ggvenn)

6. 功能富集分析

  • 目的

    探究差异基因参与的功能，筛选关键基因

  • 哪些数据库可用？

    • GO数据库富集分析==⇒找细胞相关通路(BP/MF/CC)

    • KEGG数据库富集分析==⇒找代谢相关的通路

    • 自定义数据集==⇒找自己感兴趣or符合自己预期的通路，推荐最好是自己平时整理收集的文献中通路

  • GSEA富集分析

    评估通路活性的最常用的方法之一

  • 通路结果太多怎么办？

    • 优先保证pvalue至少小于0.05，适当严格，或者卡padj；

    • 有先验知识的情况下(即有文献报道)，挑选符合自己预期或者相关的通路；无先验知识情况下，按照pvalue排序结果即可，或者依据富集通路网络来看，基于通路crosstalk挑选/排序；

    • 挑选多少通路合适？一般研究文章中，为了图片美观，挑选20个通路来展示即可，如果相关的通过的确非常多？建议进行聚类和绘制进化树类似的结果，还可以基于通路的层级关系进行聚类。

    • 当然我们也可以基于一些算法进一步进行分析缩小选择范围，例如：基于单基因GSEA分析进行筛选，基于pathfindR的蛋白互作子网活性评分算法等。

  • 工具选择(常用)

    R packages(eg：clusterProfiler/GseaVis)

7. 差异互作网络

  • 目的

    挖掘差异基因/蛋白之间的互作关系，基于互作先验信息快速筛选与研究相关的关键基因，此外还可以用于寻找更多关键基因和调控机制理解

  • 数据库有哪些？

    • STRING： STRING数据库是一个基于公共数据库和文献信息的蛋白质相互作用网络数据库。它收集了多个公共数据库，包括UniProt、KEGG、NCBI和Gene Ontology等，整合了这些数据并生成一个全面的蛋白质相互作用网络数据库。

    • BioGRID： BioGRID是一个基因相互作用的数据库，数据主要来源于文献挖掘，包括常规实验结果和高通量数据验证结果。BioGRID共收集了 2, 306, 028 种基因蛋白相互作用、29, 417 种化学相互作用和 1, 128, 339 种翻译后修饰；

    • HitPredict： 网站收集了Intact，BIOGRID和HPRD数据库中高通量实验或小规模实验得到的蛋白质互作关系，综合了三大数据库的内容，数据准确全面。此外，还有根据互作得分估计的蛋白质相互作用。

  • 网络太大一团乱，不知道怎么挑选基因？

    • 如果自己有关注的基因，可以只看关注基因的互作网络；如果没有，可以只看互作score最高的前100条关系对的网络(毕竟最可信)；其次也可以看度最大的前100个基因构成的网络；

    • 提高互作score的阈值，默认是400；

    • 网络基于子网络连通度进行聚类，关注子模块也可以。

  • 工具选择(常用)

    R packages(eg：STRINGdb/ggraph)

8. 可变剪切

  • 目的

    探索与研究相关的可变剪切分布，尤其是探究关注基因的在研究中是否发生了可变剪切，具体是哪一种可变剪切？

  • 可视化？

  • 工具选择(常用)

    R packages(eg：rMATS/ggplot2/rmats2sashimiplot)

1.4 标准分析报告必知必会

以下是层级递进的，需保证之前的分析，才能看单独部分

  1. 数据测序情况如何？===⇒ 看质控结果
  2. 数据能不能用？===⇒ 看比对结果+重复性

  3. 测序结果与实验结果是否一致？

    • 如果实验上发现某个marker有改变，看差异分析结果中相关marker/基因的FC和pvalue等信息是否达标；
    • 如果实验上发现某个现象有改变，看差异基因显著富集通路中有没有与该现象有关的通路

  4. 测序结果与实验结果相反？

    • 首先我们需要知道的是，用qPCR验证RNA-seq的结果，验证的是差异趋势，即上调还是下调，而非差异倍数。其次，目前研究显示：qPCR和RNA-seq出现一定程度的不一致（30~40%左右）是正常的也是合理的，RNA-seq与qPCR相关性在0.8左右，有15.1%-19.4%的RNA-seq结果与qPCR对应不上，non-concordant序列中有1.6%-2.8%与qPCR结果完全相反，因此两者不能一一对应；

    • 什么原因带来的不一致？

      • RNA-seq几乎覆盖基因的外显子区域，因此获得的基因表达量综合考虑了该基因所有转录本的表达水平，而qPCR实验则是设计一段引物扩增基因特定区域，通常只在80-200bp左右，不能代表所有的转录本。
      • 此外转录组由于建库和测序过程中不可避免的引入扩增和测序duplication，会对定量的准确性造成影响；因此二者在估算基因表达水平趋势时会产生不一致的情况；

    • 怎么办？哪一个结果是准确的？我该相信哪一个结果？ 当两者趋势不同时，首先我们要核实qPCR和转录组测序的设计和操作及结果是否准确，下机数据是否有效。确保数据无误后，可以按照如下方案进行甄别：

      当两者趋势不同时，确定实验和测序并无问题时，两个的结果都是正确的，如果测序和实验的RNA是同一份分出来的两部分，优先相信实验的结果；特别提醒，如果想保证测序和实验的一致性的成功率，确保样本的RNA量足够，测序和实验的RNA最好是同一份，临床样本存在很强的异质性，用其他样本来验证测序结果，一致性的概率不是很高。

    • 还有没有补救办法？

      既然已经知道了不一致的原因，此时可以看看研究关注的基因具体是哪一个转录本或者亚型，如果确实仅关注转录本，可尝试更换基因表达定量的方法，不要用基因来进行表达定量，可以尝试转录本定量。

1.5其他个性化挖掘

  • 时序/一致性聚类分析，适用于大于两个时间点/分组的数据，旨在挖掘一致性变化的基因和功能

  • 共表达网络分析(WGCNA)，适用于样本量较大(至少>12)的数据且具有临床表型的数据(分组也可以)，旨在挖掘共临床表型相关的共表达基因模块，进一步挖掘关键基因

  • APA/IPA分析(与可变剪切类似，但是研究对象关注的转录本3'端的polyA序列)，无特殊要求，目前研究发现APA异常调控与恶性肿瘤等复杂疾病密切相关

  • 融合分析/突变分析，适用于肿瘤样本或者遗传病症的样本等

  • 转录因子调控网络，无特殊要求，旨在挖掘上游的调控基因和机制

  • 多组学联合分析

    • 转录组+蛋白组
    • 转录组+表观组(Chip-seq/Cut&Tag/m6A)
    • 转录组+翻译组
    • 转录组+代谢组
    • 转录组+单细胞转录组/空间转录组
    • 转录组+WES/WGS

  • 预测模型构建，需要足够的样本量，也可以采取公共数据进行分析，多是临床样本项目

  • 拟时序分析，适用于具有时序或多分组的样本(推荐大于3组)，该方法与时序分析不同，其目的是为了实现与单细胞拟时序分析类似的目的：根据基因表达的峰值的时间点来对基因表达的先后顺序进行排序，目前该方法应用并不是很多。

1.6 后续实验怎么设计开展？

挑选基因&建议

  • 不选取表达量都很低的基因；例如：样品A，B，C，某个基因分别是0.01，0.03，0.04，这种不选，都很低，0.01和0.03虽然差了三倍，但是变差很大，随便变差一点可能就没有2倍或者3倍；宁愿选择表达量是10，20，30这样的；另外，可以选择0.03和10这样，有高表达的。

  • 可以选择一些与文章研究目标相关的基因；比如抗旱，抗虫，发育等相关基因表达确实也是预期这样的；

  • 要求选择同一批样本来做qPCR验证；最好是转录组剩余的RNA，建议转录测序时候备注，多提取一些。

  • 选择的基因数量建议不低于15个；

  • 选择稳定的，有参考文章的内参基因；提供拉丁名可以从NCBI找；

  • qPCR的实验操作和实验设计要符合规范，保证生物学重复和技术性重复的数量；

  • 一般需要提供基因的CDS序列 和内参；CDS序列的查找，如果是有参分析，通过基因ID号，在NCBI选择对应物种查找序列；如果是无参分析，转录结果里面有组装好的unigenes，直接可以搜索ID获得。

实验方案

1.7 想学数据分析？

NCBI官方分析指南

清华大学生命科学学院-本科课程