转录组测序及分析报告
项目信息
合同编号：DEMO-2021-01-29-xx
客户姓名：Client-name
客户单位：Unit-address
1. 分析流程
1.1. 建库测序流程
从RNA样品提取到最终数据获得，样品检测、建库、测序等每一环节都会直接影响数据的数量和质量，从而影响后续数据分析的结果。为从源头保证测序数据准确可靠，在数据的所有生产环节都严格把关，从根源上确保高质量数据的产出。建库测序的流程：
Total RNA 样本检测
RNA 富集
双链cDNA合成
末端修复、加A和接头
片段选择和 PCR 扩增
文库质检
Illumina测序
1.2. 信息分析流程
RNA-seq的核心是基因表达差异的显著性分析，使用统计学方法，比较两个条件或多个条件下的基因表达差异，从中找出与条件相关的特异性基因，然后进一步分析这些特异性基因的生物学意义，分析过程包括质控、比对、定量、差异显著性分析、功能富集等环节。信息分析流程如下图所示：
2. 信息分析
2.2. 测序数据质量控制
对原始测序数据及去除接头后的可用数据进行质量评估。测序数据一般为双端测序，因此，每个测序样本会有两个测序结果。
评估的具体内容见：
RawData-fastqc 文件链接： /result/qc/qc_rawdata/*.html
CleanData-fastqc 文件链接： /result/qc/qc_cleandata/*.html
Fastqc 格式补充说明： /result/qc/qc_Supplement.html
2.3. 参考基因组比对
测序片段（fragments）是mRNA随机打断的，为了确定这些一段由哪些基因转录来，需要将质控后的clean reads比对到参考基因组上。使用HISAT2软件将Clean Reads与参考基因组进行快速精确的比对，获取Reads在参考基因组上的定位信息[4]。HISAT2软件官方手册。
如果参考基因组组装的较为完善，而且所测物种与参考基因组一致，且相关实验不存在污染，那么实验所产生的测序reads成功比对到基因组的比例会高于70% (Total Mapped Reads or Fragments)。本项目所用参考基因组为 hg38 ，下载链接：Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz 。
结果文件：
各个样本的比对情况统计文件：/result/map_stat/*.flagstat.txt
2.4. 定量分析
2.4.1. 样本间相关性
生物学重复通常是任何生物学实验所必须的，目前主流期刊也基本要求生物学重复。生物学重复主要有两个用途：一个是证明所涉及的生物学实验操作不是偶然，而是可重复的。另一个是为了确保后续的差异基因分析得到更可靠的结果。样品间基因表达水平相关性是检验实验可靠性和样本选择是否合理的重要指标。相关系数越接近1，表明样品之间表达模式的相似度越高。Encode计划建议皮尔逊相关系数的平方(R2)大于0.92(理想的取样和实验条件下)。具体的项目操作中，我们要求生物学重复样品间R2至少要大于0.8，否则需要对样品做出合适的解释，或者重新进行实验。根据各样本所有基因的表达值计算组内及组间样本的相关性系数，绘制成热图，可直观显示组间样本差异及组内样本重复情况。样本间相关性系数越高，其表达模式越为接近，样本相关性热图如下图所示。
图 1 样本间相关性热图
图中横纵坐标为各样本相关系数的平方
结果文件：
样本间相关性热图结果：Quant/cor_pheatmap*
2.4.2. 主成分分析
主成分分析（PCA）也常用来评估组间差异及组内样本重复情况，PCA采用线性代数的计算方法，对数以万计的基因变量进行降维及主成分提取。我们对所有样本的基因表达值进行PCA分析，如下图所示。理想条件下，PCA图中，组间样本应该分散，组内样本应该聚在一起。
图 2 主成分分析结果图
图中横坐标为第一主成分，纵坐标为第二主成分
结果文件：
主成分分析结果：Quant/pca*
2.5. 差异分析
基因表达定量完成后，需要对其表达数据进行统计学分析，筛选样本在不同状态下表达水平显著差异的基因。差异分析主要分为三个步骤。
首先对原始的readcount进行标准化（normalization），主要是对测序深度的校正。
然后统计学模型进行假设检验概率（pvalue）的计算
最后进行多重假设检验校正，得到FDR值（错误发现率，padj是其常见形式)[1-2]。
针对不同的实验情况，我们选用合适的软件进行基因表达差异显著性分析，具体如下表所示。
表1 表达差异分析所用软件及差异基因筛选标准
类型 | 软件 | 标准化方法 | pvalue计算模型 | FDR计算方法 | 差异基因筛选标准 |
|---|---|---|---|---|---|
有生物学重复 | DESeq2(Anders et al, 2014) | DESeq | 负二项分布 | BH | |log2(FoldChange)| > 0 & padj < 0.05 |
无生物学重复 | edgeR(Robinson et al, 2010) | TMM | 负二项分布 | BH | |log2(FoldChange)| > 1 & padj < 0.05 |
若按照以上标准筛选得到的差异基因过少（低于100），很有可能导致后面的功能富集分析没有显著性结果，所以，我们会根据项目的具体情况，适当地降低筛选差异基因的阈值标准。若项目实验只关注某几个基因的表达情况（如基因敲除），不在意富集结果，从下面的差异分析表格中筛选关注的那几个基因即可。
一般来说，如果一个基因在两组样品中的表达量差异达到两倍以上，我们认为这样的基因是具有表达差异的。为了判断两个样品之间的表达量差异究竟是由于各种误差导致的还是本质差异，我们需要对所有基因在这两个样本中的表达量数据进行假设检验。而转录组分析是针对成千上万个基因进行的，这样会导致假阳性的累积，基因数目越多，假设检验的假阳性累积程度会越高，所以引入padj对假设检验的P-value进行校正，从而控制假阳性的比例[3]。
差异基因的筛选标准是非常重要的，我们给出的标准|log2(FoldChange)| > 1 & padj< 0.05是常用的经验值，在实际项目中可以根据情况灵活选择。例如，差异倍数可以选择1.5倍，也可以选择3倍，padj常用的阈值包括0.01、0.05、0.1等。若按照以上标准筛选得到的差异基因过少，很有可能导致后⾯的功能富集分析没有显著性结果。若项目实验只关注某几个基因的表达情况（如基因敲除），不在意富集结果，从下面的差异分析表格中筛选关注的那几个基因即可。反之，如果得到的差异基因数目过多，不利于后续目标基因的筛选，这个时候可使用更严格的阈值标准进行筛选，则可以使用更严格的阈值标准进行筛选。
2.5.1. 差异基因的筛选
通过Deseq2进行差异分析，我们通常采用 |log2FC|>1 & padj < 0.05 进行差异基因的筛选，随后对差异基因进行注释，得到包含注释信息的差异基因列表。
结果文件：
差异基因列表及相关注释信息（总的结果）：result/Enrichment/Allgene_anno_ALL.xls
差异基因列表及相关注释信息（筛选结果）：result/Enrichment/Allgene_anno.xls
Differential/Allgene_anno*.xls表头
表头 | 说明 |
|---|---|
ENSEMBL | 差异基因的ENSEMBL名 |
pvalue | 差异基因的置信度计算结果 |
padj | 差异基因的多重校验FDR |
log2FC | Treat组 vs Control组 差异倍数 的log2标准化结果 |
FC | Treat组 vs Control组 差异倍数 |
log2FC_abs | Treat组 vs Control组 差异倍数 的log2标准化结果的绝对值（此列便于筛选log2FC阈值） |
FC_HvsL | 高表达组 vs 低表达组 差异倍数 （此列便于筛选FC阈值） |
change | 使用本次分析的阈值，对差异基因的上下调标记 |
SYMBOL | 差异基因的SYMBOL名 |
ENTREZID | 差异基因的ENTREZID号 |
GENENAME | 差异基因的基本描述信息 |
baseMean | 差异基因的表达量标准化后的平均值 |
Samples* | 样本的原始表达矩阵表达量结果 |
Samples*_normal | 样本的表达矩阵标准化后的结果 |
2.5.2. 差异基因的热图聚类
将所有比较组的差异基因取并集之后作为差异基因集。两组以上的实验，可对差异基因集进行聚类分析，将表达模式相近的基因聚在一起。我们采用主流的层次聚类对基因的表达值进行聚类分析，对行（row）进行均一化处理（Z-score）。热图中表达模式相近的基因或样本会被聚集在一起，每个方格中的颜色反映的不是基因表达值，而是表达数据的行进行均一化处理后得到的数值（一般在-1到1之间），所以热图中的颜色只能横向比较（同一基因在不同样本中的表达情况），不能纵向比较（同一样本不同基因的表达情况）。结果文件中既有组间的聚类，也有样品间的聚类。结题报告展示了样品间的聚类，具体如下图所示。
图 3 差异表达基因聚类热图
图中横坐标为样品名，纵坐标为差异基因归一化后的数值，颜色越红，表达量越高，越蓝，表达量越低。
结果文件：
差异基因的热图结果：Differential/heatmap/
2.5.3. 差异基因的火山图分布
火山图可直观展示每个比较组合的差异基因分布情况，如下图所示。图中横坐标表示基因在处理和对照两组中的表达倍数变化(log2FoldChange)，纵坐标表示基因在处理和对照两组中表达差异的显著性水平(-log10padj或-log10pvalue)。为上调基因用红色点表示，下调基因用蓝色点表示。
图 4 差异基因火山图
图中横坐标为log2FoldChange值，纵坐标为-log10padj或-log10pvalue，蓝色的虚线表示差异基因筛选标准的阈值线
结果文件：
差异基因的火山图结果：Differential/volcano/volcano.png
2.6. 富集分析
我们根据基因表达量分析得到差异基因之后，必须进一步落到基因的功能上来。对于转录组分析而言，往往涉及到成千上万个基因，这会使分析变得很复杂。解决思路是将一个基因列表分成多个部分，从而减少分析的复杂度。为了解决怎么分成不同类，通常会对基因功能进行富集分析,
期望发现在生物学过程中起关键作用的生物通路,
从而揭示和理解生物学过程的基本分子机制。功能富集分析可以将成百上千个基因、蛋白或者其他分子分到不同的通路中，以减少分析的复杂度。另外，在两种不同实验条件下，激活的通路显然比简单的基因或蛋白列表更有说服力。基因功能富集分析首先要构建基因集（gene set，如GO和KEGG数据库等），也就是基因组注释信息进行分类。然后再把我们的目标基因集（差异基因集或者其他基因集）映射到背景基因集上，注意区分注释与富集。
我们采用clusterProfiler软件对差异基因集进行GO功能富集分析，KEGG通路富集分析等。富集分析基于超几何分布原理，其中差异基因集为差异显著分析所得差异基因并注释到GO或KEGG数据库的基因集，背景基因集为所有进行差异显著分析的基因并注释到GO或KEGG数据库的基因集。富集分析结果是对每个差异比较组合的所有差异基因集、上调差异基因集、下调差异基因集进行富集。本报告中展示的表格是选取某一个比较组合的富集分析结果，图片是所有组合的富集分析结果。
图 5 基因富集分析原理图
2.6.1. 富集分析结果文件
结果路径 | 结果说明 |
|---|---|
GO富集分析结果 | |
Results/*enrich_*/gene.ego_all-p.adjust1.00.csv | GO富集结果列表（所有结果） |
Results/*enrich_*/gene.ego_all-p.adjust0.05.csv | GO富集结果列表（按p.adj<0.05筛选后） |
Results/*enrich_*/gene.ego_ALL.csv | GO富集结果列表（MF、BP、CC所有结果） |
Results/*enrich_*/gene.GO-*-barplot.p* | GO富集分析柱状图 |
Results/*enrich_*/gene.GO-*-dotplot.p* | GO富集分析散点图 |
Results/*enrich_*/gene.GO-*-DAG.p* | GO富集分析DAG图 |
KEGG富集分析结果 | |
Results/*enrich_*/gene.KEGG.csv | KEGG富集结果列表（所有） |
Results/*enrich_*/gene.KEGG_significant.csv | KEGG富集结果列表（按p.adj<0.05筛选后） |
Results/*enrich_*/gene.KEGG-*-barplot.p* | KEGG富集分析柱状图 |
Results/*enrich_*/gene.KEGG-*-dotplot.p* | KEGG富集分析散点图 |
结果文件夹：
all分析结果文件夹：result/Enrichment/all/
up分析结果文件夹：result/Enrichment/up/
down分析结果文件夹：result/Enrichment/down/
all网页预览图：result/Enrichment/all-pdf.html
up网页预览图：result/Enrichment/up-pdf.html
down网页预览图：result/Enrichment/down-pdf.html
说明：
all/up/down分别对应总差异基因，上调差异基因，下调差异基因进行对应的富集分析。
表头说明： （Results/*enrich_*/gene.ego_*.csv GO富集结果列表）
ID | 对应GO数据库中的ID |
|---|---|
ONTOLOGY | 分子功能（Molecular Function），生物过程（biological process）和细胞组成（cellular component） |
Description | GO的描述 |
GeneRatio | 对应GO 差异基因数 / 能够对应到GO数据库中同类型的差异基因数 |
BgRatio | 对应GO包含对应物种的基因数 / GO数据库中包含对应物种的基因数 |
pvalue | 富集分析得到的p-value |
p.adjust | 校正后的p-value |
qvalue | 富集分析得到的qvalue |
Count | 富集基因数目 |
ENTREZID | 富集基因列表（ENTREZID） |
SYMBOL | 富集基因列表（SYMBOL） |
表头说明： （Results/*enrich_*/gene.KEGG*.csv KEGG富集结果列表）
ID | 对应PATHWAY数据库中的ID |
|---|---|
Description | PATHWAY的描述 |
GeneRatio | 对应PATHWAY 差异基因数 / 能够对应到PATHWAY数据库中的差异基因数 |
BgRatio | 对应PATHWAY包含对应物种的基因数 / PATHWAY数据库中包含对应物种的基因数 |
pvalue | 富集分析得到的p-value |
p.adjust | 校正后的p-value |
qvalue | 富集分析得到的qvalue |
Count | 富集基因数目 |
ENTREZID | 富集基因列表（ENTREZID） |
SYMBOL | 富集基因列表（SYMBOL） |
2.6.1. GO功能富集分析
GO(Gene Ontology)是描述基因功能的综合性数据库，可分为生物过程（biological process）和细胞组成（cellular component）分子功能（Molecular Function）三个部分。GO功能富集以padj小于0.05作为为显著性富集的阈值，富集结果见结果文件。
从GO富集分析结果中，选取最显著的30个Term绘制柱状图进行展示，若不足30个，则绘制所有Term，按生物过程、细胞组分和分子功能三大类别及差异基因上下调分类画的柱状图。
有向无环图(Directed Acyclic Graph，DAG)为差异基因GO富集分析结果的图形化展示方式。图中，分支代表包含关系，从上至下所定义的功能范围越来越小，选取每个差异比较组合的GO富集结果最显著性前5位的GO Term作为有向无环图的主节点，并通过包含关系，将相关联的GO Term一起展示，颜色的深浅代表富集程度。我们的项目中分别绘制生物过程、分子功能和细胞组分的DAG图。
图 6 GO富集分析柱状图
图中纵坐标为GO Term，横坐标为GO Term富集的显著性水平，数值越高越显著
图 7 GO富集分析散点图
图中横坐标为注释到GO Term上的差异基因数与差异基因总数的比值，纵坐标为GO Term
图 8 GO富集分析DAG图
每个节点代表一个GO术语，方框代表的是富集程度为TOP5的GO，颜色的深浅代表富集程度，颜色越深就表示富集程度越高，每个节点上展示了该TERM的名称及富集分析的padj
2.6.2. KEGG通路富集分析
KEGG(Kyoto Encyclopedia of Genes and Genomes)是整合了基因组、化学和系统功能信息的综合性数据库。KEGG通路富集以padj小于0.05作为显著性富集的阈值，富集结果见结果文件。
从KEGG富集结果中，选取最显著的20个KEGG通路绘制柱状图进行展示，若不足20个，则绘制所有通路，如下图所示。图中横坐标为通路富集的显著性水平，数值越高越显著，纵坐标为KEGG通路。
从KEGG富集结果中，选取最显著的20个KEGG通路绘制散点图进行展示，若不足20个，则绘制所有通路，如下图所示。图中横坐标为注释到KEGG通路上的差异基因数与差异基因总数的比值，纵坐标为KEGG通路，点的大小代表注释到KEGG通路上的基因数，颜色从红到紫代表富集的显著性大小。
图 9 KEGG富集分析柱状图
图中横坐标为通路富集的显著性水平，数值越高越显著，纵坐标为KEGG通路。
图 10 KEGG富集散点图
图中横坐标为注释到KEGG通路上的差异基因数与差异基因总数的比值，纵坐标为KEGG通路