鸿源韬生物

RNA-Seq生信分析报告

1.项目简介

1.1 样本信息

1.2 实验原理及流程

转录组测序是基于二代测序平台，研究特定组织或细胞在某个时期转录出来的所有mRNA，对理解生物体的发育和疾病的发生具有重要作用。随着基因测序技术的发展以及测序成本的降低，RNA-seq凭借高通量、高灵敏度、应⽤范围⼴等优势，已成为转录组研究的主要方法。转录组研究能够从整体水平研究基因的表达和功能，揭示特定生物学过程以及疾病发生过程中的转录层面的分子机理(Wang Z et al., 2009)。

1.3 实验流程

RNA-Seq实验主要步骤有七步:
1.样品采集与RNA提取，使用 TRIzol提取总 RNA；
2.RNA 质量检测，使用Agilent 2100 Bioanalyzer测定RNA完整性；
3.mRNA富集，通过Oligo(dT)磁珠富集带有polyA 尾的mRNA；
4.片段化，在Fragmentation Buffer中用二价阳离子将得到的mRNA随机打断；
5.文库构建，以mRNA片段为模板通过反转录、纯化、末端修复、加A尾以及接头连接、PCR扩增以及纯化产物得到测序文库；
6.上级测序，使用Qubit2.0Fluorometer进行初步定量以及Agilent 2100 Bioanalyzer测定文库插入片段分布，质控合格后使用illumina测序，产生PE 150bp的测序数据；
7.生信分析。

RNA-Seq 实验原理(Wang Z et al., 2009)

1.4 分析流程

获得测序原始数据（raw data）后，首先对原始数据进行过滤，获得高质量的测序数据（clean data），将测序数据（clean data）比对到项目物种的参考基因组上。根据比对结果对各样本基因表达量进行统计，并进行差异基因分析。最后对差异表达基因（Differential Expression Genes, DEGs）进行GO、KEGG富集分析。

RNA-Seq 生物信息学分析流程

2. 数据质控

我们交付的原始数据为fastq（简称fq）格式文件的压缩包，文件名后缀通常为 “.fq.gz”。交付数据前我们会计算每个压缩文件的md5值。在您拿到数据之后，请您先校>验每个压缩文件的md5值，Linux下可以在数据目录使用“md5sum -c <*md5.txt>”命令进行校验，Windows下可使用hashmyfiles等校验工具，如发现压缩文件md5值与附在数据文件目录下的md5文档中的不一致则说明文件可能在传输的过程中被损坏。数据文件大小为文件占用磁盘空间的大小，文件的大小通常与磁盘格式、压缩比例等因素有关，与测序数据量（碱基数）的多少无对应关系，因此对应PE测序的 read1和read2两个文件大小也可能不相同。

将高通量测序得到的原始图像数据经过Base Calling 转化为序列数据，即FASTQ格式，得到最原始的测序数据文件。FASTQ 格式文件可记录所测读段（read）的碱基及其质量分数。FASTQ 格式以测序读段为单位进行存储，每条读段占 4 行，第一行是序列标识（read ID）以及相关的描述信息，以“@” 开头；第二行即为碱基序列，长度由测序策略决定；第三行以“+”开头，后面是序列标示符、描述信息，或者什么也不加； 第四行是测序质量值（phred），与第二行一一对应，phred值以ASCII码标记，对应的 ASCII 值减去33，即为第二行对应碱基的测序质量值，示例如下：

@HWI-ST1276:71:C1162ACXX:1:1101:1208:2458 1:N:0:CGATGT
NAAGAACACGTTCGGTCACCTCAGCACACTTGTGAATGTCATGGGATCCAT
+
#55???BBBBB?BA@DEEFFCFFHHFFCFFHHHHHHHFAE0ECFFD/AEHH

测序错误率用e表示， 平台测得数据的碱基质量值用Qphred表示，则有：Qphred=-10log10(e)。软件中碱基识别正确率与Phred分值之间的简明对应关系见下表：

测序Reads的错误率往往会随着测序接近尾声而升高，这是由测序过程中化学试剂的消耗造成共有的特征。

2.1 原始数据质控

RNA-Seq实验基于第二代测序（NGS）平台完成，采用双端测序文库构建策略（插入片段~300 bp）。我们需要对原始测序数据进行质量评估与过滤，以确保后续分析的可靠性。首先，我们使用FastQC（version 0.12.1）(Andrews, 2010)对原始测序数据（raw data）进行全局质量分析，包括碱基质量分布（Phred score）、碱基组成平衡性（base content uniformity）、重复序列比例（duplication level）及GC含量偏差等指标，以全面评估测序质量。
我们使用Fastp（version 0.24.0）(Chen et al., 2018)对原始测序数据进行以下过滤操作。
接头序列去除：识别并切除双端reads中的接头序列；
低复杂度序列过滤：剔除含模糊碱基（N碱基占比≥10%）的reads；
动态质量修剪：通过滑动窗口法（5 bp窗口步长）评估局部序列质量，当窗口平均Phred score小于20时，执行3'端截断；
长度筛选：保留长度≥25 bp的paired-end reads，长度不足的reads及其匹配reads（R1/R2）均被排除。
原始和过滤后质控结果请详见result/1.qc文件夹，raw为原始数据质控结果，clean为过滤后质控结果。

panicle_1.R1panicle_1.R2panicle_2.R1panicle_2.R2yleaves_1.R1yleaves_1.R2yleaves_2.R1yleaves_2.R2

❮

❯

图2.1 各个样本平均测序碱基质量分数，横坐标代表150 bp长度序列中各个位置，纵坐标为该位置平均的碱基质量值Q；盒形图中间的红线表示中位数（median value）；黄色部分代表四分位距（25-75%）；上下分割线代表 90%和 10%的上下临界值；蓝色的线代表碱基质量的平均值。

panicle_1.R1panicle_1.R2panicle_2.R1panicle_2.R2yleaves_1.R1yleaves_1.R2yleaves_2.R1yleaves_2.R2

❮

❯

图2.2 各个样本碱基平衡性，图中四条线代表A T C G在每个位置平均含量。理论上，A和T应该相等，G和C应该相等，且4种碱基平行且接近分布。正常情况下四种碱基的出现频率应该是接近的，而且没有位置差异。因此好的样本中四条线应该平行且接近。当部分位置碱基的比例出现 bias 时，即四条线波动较大时可能存在测序数据或者文库污染。如果所有位置的碱基比例一致的表现出bias 时，即四条线平行但分开，往往代表文库有 bias (建库过程或本身特点)，或者是测序中的系统误差。测序刚开始由于测序仪状态不稳定，在15bp之前很可能出现波动。

panicle_1.R1panicle_1.R2panicle_2.R1panicle_2.R2yleaves_1.R1yleaves_1.R2yleaves_2.R1yleaves_2.R2

❮

❯

图2.3 各个样本重复序列水平，测序深度越高，越容易产生一定程度的重复（duplication），这属于正常的现象。但如果duplication 的程度很高，就提示我们可能有 bias 的存在（如建库过程中由于 PCR 扩增引起的duplication）。横坐标为 reads 重复的次数，纵坐标为重复次数对应的 reads 占 unique reads 的比例，以unique reads 的总数作为 100%。这里，我们仅对文件前 2000000 个reads 进行统计：对长度小于75bp 的reads 将其截短为 50bp，用于统计重复。

panicle_1.R1panicle_1.R2panicle_2.R1panicle_2.R2yleaves_1.R1yleaves_1.R2yleaves_2.R1yleaves_2.R2

❮

❯

图2.4 各个样本序列平均GC含量分布图，横坐标是GC含量，纵坐标是每条序列GC含量对应的reads数量。蓝色的线是程序根据经验分布给出的理论值，红色是真实值，理想情况下两者应该接近，当红色的线出现双峰，说明可能混入了其他物种的DNA序列。

2.2 过滤后数据质控

这里展示Fastp过滤后的数据质控结果，图片内容与上面raw data类似。

panicle_1.R1panicle_1.R2panicle_2.R1panicle_2.R2yleaves_1.R1yleaves_1.R2yleaves_2.R1yleaves_2.R2

❮

❯

图2.5 各个样本平均测序碱基质量分数，横坐标代表150 bp长度序列中各个位置，纵坐标为该位置平均的碱基质量值Q；盒形图中间的红线表示中位数（median value）；黄色部分代表四分位距（25-75%）；上下分割线代表 90%和 10%的上下临界值；蓝色的线代表碱基质量的平均值。

panicle_1.R1panicle_1.R2panicle_2.R1panicle_2.R2yleaves_1.R1yleaves_1.R2yleaves_2.R1yleaves_2.R2

❮

❯

图2.6 各个样本碱基平衡性，图中四条线代表A T C G在每个位置平均含量。理论上，A和T应该相等，G和C应该相等，且4种碱基平行且接近分布。正常情况下四种碱基的出现频率应该是接近的，而且没有位置差异。因此好的样本中四条线应该平行且接近。当部分位置碱基的比例出现 bias 时，即四条线在某些位置波动较大时，可能测序数据或者文库存在污染。当所有位置的碱基比例一致的表现出bias 时，即四条线平行但分开，往往代表文库有 bias (建库过程或本身特点)，或者是测序中的系统误差。一般测序的时候，刚开始测序仪状态不稳定，在15bp之前很可能出现波动。

panicle_1.R1panicle_1.R2panicle_2.R1panicle_2.R2yleaves_1.R1yleaves_1.R2yleaves_2.R1yleaves_2.R2

❮

❯

图2.7 各个样本重复序列水平，测序深度越高，越容易产生一定程度的重复（duplication），这属于正常的现象。但如果duplication 的程度很高，就提示我们可能有 bias 的存在（如建库过程中由于 PCR 扩增引起的duplication）。横坐标为 reads 重复的次数，纵坐标为重复次数对应的 reads 占 unique reads 的比例，以unique reads 的总数作为 100%。这里，我们仅对文件前 2000000 个reads 进行统计：对长度小于75bp 的reads 将其截短为 50bp，用于统计重复。

panicle_1.R1panicle_1.R2panicle_2.R1panicle_2.R2yleaves_1.R1yleaves_1.R2yleaves_2.R1yleaves_2.R2

❮

❯

图2.8 各个样本序列平均GC含量分布图，横坐标是GC含量，纵坐标是每条序列GC含量对应的reads数量。蓝色的线是程序根据经验分布给出的理论值，红色是真实值，理想情况下两者应该接近，当红色的线出现双峰，说明可能混入了其他物种的DNA序列。

2.3 数据过滤结果统计

我们对数据过滤结果进行统计，如下表所示：

表 2.1数据过滤结果统计：
Sample：样品名称；
Raw_Total_Reads/Clean_Total_Reads：过滤前后样本总reads数量，单位为百万；
Raw_Total_Bases/Clean_Total_Bases：过滤前后样本总碱基数量，单位为百万；
Raw_Q20_Rate/Clean_Q20_Rate：过滤前后样本Q20碱基比例；
Raw_Q30_Rate/Clean_Q30_Rate：过滤前后样本Q30碱基比例；
Raw_GC_Content/Clean_GC_Content:过滤前后样本GC含量。

3. 比对参考基因组

我们将各样品过滤后的clean data的reads与参考基因组进行比对，获取Reads在参考基因组上的定位信息，这里使用的软件是Hisat2(version 2.2.1)(Kim D et al., 2015)。

3.1 比对参考基因组情况

表 3.1比对结果统计：
Sample：样品名称；
clean_reads：clean后reads总数；
prop_map_reads：完美比对的reads总数，PE两端reads比对到同一条序列，且根据比对结果推断的插入片段大小符合设置的阈值；
prop_map %：完美比对reads百分比；
MAPQ30：MAPQ值大于30的reads数。

3.2 比对可视化

软件比对所得结果为bam格式文件(位于report/result/2.map文件夹中)，bam文件是压缩的⼆进制⽂件，无法直接作为文本打开查看。由于bam文件数据较大，我们将其转为较小的bw格式文件。客户可以结合物种参考基因组和注释文件使用IGV (Integrative Genomics Viewer) 浏览器对bam、bw、bed等文件进行可视化浏览。IGV浏览器使用方法可参考我们提供的使用说明文档IGV快速上手

4. 基因表达定量

根据基因比对在参考基因组上的位置信息，从而统计每个基因从起始到终止范围内覆盖的reads数。分别过滤掉比对质量值低于20的reads、非成对比对上的reads、比对到基因组多个区域的reads。我们使用featureCounts(version 2.0.6)(Liao Y et al., 2014)对各样本比对结果进行基因表达量统计。本节内容请详见report/result/3.exp文件夹，all_counts.txt为各个样本的基因表达量统计结果，allexp_fpkm.tsv为各个样本的FPKM矩阵。

4.1 FPKM矩阵

FPKM（Fragments Per Kilobase of exon model per Million mapped fragments）是转录组测序中用于校正基因长度及测序深度的相对定量指标，计算公式定义为FPKM=（比对到基因i的Fragment数/基因i外显子总长度）/（所有基因总外显子长度/10^6）。该公式通过同时校正基因长度和测序深度，消除不同样本间测序数据量的系统性偏差。我们将featurecount得到的各个基因对应reads数量进行计算得到FPKM矩阵。

表 4.1 FPKM矩阵：
geneid：基因ID；
gene_name：基因名称；
样本名：样本对应的各个基因的FPKM值。

4.2 样本表达量

allexp_boxplot

❮

❯

图4.1 各个样本各个基因FPKM值箱线图：
横坐标为样本名称，纵坐标为各个基因的FPKM值。箱线图最顶端的横线代表上边缘，最底端的横线代表下边缘，箱线图中间的横线代表中位数，箱线图的箱体代表25%和75%分位数。

allexp_pca

❮

❯

图4.2 各个样本各个基因FPKM值PCA图：
主成分分析是将原来较多维度的指标 (各样本基因FPKM值)，降维到较低的维度（二维），来研究样品间的主成分关系。二维PCA分析结果中，会展示主成分1(PC1) 和主成分2(PC2)分别作为 X 轴和 Y 轴的散点图，每个点代表 1 个样本。坐标轴上百分比代表主成分的贡献率，贡献率越大，说明该主成分对样本差异的解释能力越强。如果两个样本距离越远，则说明样本基因量整体分布的差异越大。 反之，则说明相应样本基因量整体分布模式越接近。所以，PCA 分析常用于评估样本重复性的好坏。理想情况下，生物学重复的样本应该相近，而不同组间应该可以区分开。

5. 差异基因分析

在存在多个分组且组内有生物学重复的情况下，可以对组间进行差异表达基因(Differential Expression Genes， DEGs)分析，以确定哪些基因在组间存在显著差异。如果没有差异分析则本节内容为空。

5.1 差异基因分析结果

存在组内生物学重复时，我们使用软件DESeq2(version 1.48.1)(Love MI et al., 2014)对DEGs进行分析。DESeq2计算流程： 1.归一化（Normalization）：DESeq2使用通过median-of-ratios方法计算每个样本的size factor，用它来校正不同样本的测序深度和文库差异，使得各基因间的表达水平更具可比性，注意它并不等于FPKM。
2.离散度估计（Dispersion Estimation）：DESeq2 假设每个基因的read count服从负二项分布，并估计其离散度参数。
3.差异检验（Statistical Testing）：通过对处理组与对照组的计数数据进行广义线性模型拟合，计算每个基因的log2 fold change（差异倍数）和显著性检验P值。
4.多重检验校正：使用 Benjamini–Hochberg 方法控制假阳性发现率，得到调整后的P值（padj）。
本节结果详见report/result/4.diff。差异分析用到的判断显著性阈值记录在火山图图例中。
DiffGene:“sampA_sampB_res.tsv”为各组样本差异分析结果，sampA代表实验组，sampB代表对照组； “sampA_sampB_up.csv”为sampA和sampB比较，表达量上调的基因;“sampA_sampB_down.csv”为sampA和sampB比较，表达量下调的基因。

表 5.1 各比较组差异分析结果：
geneid：基因ID；
baseMean：各组样本的平均标准化表达量，越大表示基因整体表达水平越高；
log2FoldChange：实验组vs对照组的对数2倍数变化，正值表示上调，负值表示下调，绝对值越大差异越明显；
lfcSE：log2 fold change的标准误，表示倍数变化的估计精度；
stat：Wald 检验统计量，用于计算P值；
pvalue：原始P值，反映差异显著性；
padj：校正后的P值；
样本名：各样本归一化后表达量;
diff:是否有显著差异，Up代表实验组中上调，Down代表实验组中下调，Not Sign.代表无显著差异的基因。

5.2 差异基因可视化

panicle_yleaves_volcano

❮

❯

图5.1 差异基因火山图。横坐标为log2(Fold Change)，纵坐标为-log10(p.adjust)，蓝色为显著性下调的基因，红色为显著性上调的基因，灰色为非显著性差异的基因。

panicle_yleaves_heatmappanicle_yleaves_names_heatmap

❮

❯

图5.2 差异基因热图，行代表基因列代表样本，颜色由每一行标准化后表达量觉得，Up/Down后括号标注了基因数量。"*_names_heatmap"代表上下调log2FC前20的基因，标注出基因id。

6. 基因富集分析

我们使用clusterprofiler(version 4.14.6)(Wu T. et al., 2021)进行GO和KEGG通路富集分析。富集分析结果表格未使用阈值过滤，您可以在表格中查看所有可能富集的通路。

6.1 GO富集分析

GO (Gene Ontology, http://www.geneontology.org) 是基因本体论联合会建立的将全世界所有与基因有关的研究结果进行分类汇总的综合数据库。该数据库标准化了不同数据库中关于基因和基因产物的生物学术语，适用于各物种，对基因和蛋白功能进行限定和描述。利用GO 数据库，可以对差异表达基因进行富集分析，可以找到不同条件下的差异表达基因按照其参与的BP(Biological Process, 生物过程)、MF(Molecular Function, 分子功能) 及CC(Cellular Component, 细胞组分) 三个方面进行分类注释。GO 注释有助于理解基因背后所代表的生物学意义。GO功能显著性富集分析给出与基因组背景相比，在相关基因中显著富集的GO功能条目，从而给出与差异表达基因与哪些生物学功能显著相关。该分析首先把所有相关向Gene Ontology数据库的各个term映射，计算每个term的基因数目，然后应用超几何检验，找出与整个基因组背景相比，在与差异表达基因中显著富集的GO条目。

下面展示差异表达的基因富集GO富集分析部分结果，完整结果请见/result/5.gokegg/sampA_sampB_up/down_GO_res.csv。GO富集分析完整结果请详见位于report/result/5.gokegg文件夹的*_GO_res.csv表格文件。

表6.1 GO富集分析部分结果：
ONTOLOGY:GO方面，细胞成分，生物过程或分子功能之一;
ID:GO标识符，GO ID;
Description:GO术语的文字描述;
GeneRatio:该条目基因比例，分子是富集到这个GO条目上的基因的数目，分母是所有差异表达基因的数目;
BgRatio:背景比例，分母是物种全部有GO注释的基因的数目，分子是这些基因中注释到这个GO条目上面的基因的数目;
RichFactor​​:富集因子（Enrichment Factor）= GeneRatio / BgRatio;
​​FoldEnrichment​:富集倍数（Fold Enrichment）= (富集通路基因数 / 输入基因数) / (背景通路基因数 / 背景总基因数);
​​zScore​:标准化富集得分（基于超几何分布的 Z 值）;
pvalue:富集的p值;
p.adjust:使用BH校正之后的p值;
qvalue:q值,使用FDR校正之后的p值，q-value相比于p-value更加严格,表示p-value产生假阳性的概率;
geneID:富集到这个GO条目上面的具体的基因ID;
Count:富集到这个GO条目上面的基因的数目。

panicle_yleaves_down_GO_dotplotpanicle_yleaves_up_GO_dotplot

❮

❯

图6.2 Peak关联基因GO气泡图。纵坐标是GO Term 名称，横坐标是对应GO Term 中检出的基因占背景基因的个数，颜色代表显著性，气泡大小代表该条目基因比例。

panicle_yleaves_down_GO_barplotpanicle_yleaves_up_GO_barplot

❮

❯

图6.3 Peak关联基因GO条状图。按照BP、MF、CC三个方面分别展示GO富集结果。纵坐标是GO Term 名称，横坐标值越大显著性越高，如果为0代表qvalue等于1。

6.2 KEGG富集分析

KEGG (Kyoto Encyclopedia of Genes and Genomes, http://www.genome.jp/kegg/) 是日本京都大学构建的基因组信息数据库，它将基因组序列信息与功能信息相结合，提供了一个全面的基因组功能信息资源。在PATHWAY数据库里，包括图解的细胞生化过程如代谢、膜转运、信号传递、细胞周期，还包括同系保守的子通路等信息。KEGG富集分析可以对差异表达基因进行KEGG通路富集分析。

下面展示差异表达的基因富集KEGG富集分析部分结果，完整结果请见/result/5.gokegg/sampA_sampB_up/down_KEGG_res.csv。GO富集分析完整结果请详见位于report/result/5.gokegg文件夹的*_KEGG_res.csv表格文件。