1.1 文库基本情况

每一个文库的基本情况见表1。

1.2 数据整理

样品经过上机测序，得到图像文件，由测序平台自带软件进行转化，生成 FASTQ 的原始数据（Raw Data），即下机数据。我们对每个样品的下机数据（Raw Data）分别进行统计，包括样品名、 Q30 、模糊碱基所占百分比、以及 Q20(%) 和 Q30(%) 。统计结果见表2。

1.3 数据过滤

测序数据包含一些带接头 、低质量的 Reads，这些序列会对后续的信息分析造成很大的干扰，因此需要对测序数据进行进一步过滤。数据过滤的基本情况见表3。数据过滤的标准主要包括：

• 1）采用 Cutadapt 去除 3' 端带接头的序列；
• 2）去除平均质量分数低于 Q20 的 Reads。

1.4 碱基质量分布

测序错误率受测序仪本身、测序试剂、样品等多个因素共同影响。 对于 RNASeq 技术，测序错误率分布具有两个特点：

• 1）测序错误率会随着测序序列的长度增加而升高，这是由测序过程中化学试剂的消耗导致的，是 Illumina 高通量测序平台都具有的特征；
• 2）前6个碱基的位置（即建库过程中反转录所需要的随机引物的长度）也会发生较高的测序错误率，这种错误是由随机引物和 RNA 模板的不完全结合引起的。

我们用测序数据的单碱基质量分布图评价单个位置的碱基质量。一般而言，Reads 的 5’端和 3’端的碱基质量较低，中间部分的碱基质量较高。大部分序列的碱基质量在20以上，代表测序质量较好。

samples_R1samples_R2

图1 单碱基质量分布图

注：横坐标是 Reads 中碱基位置（5’->3’），纵坐标是对应位点碱基 Q 值。红线代表中位数，蓝线代表平均数，黄色区域代表 25%-75%区间（按照四分位数划分），触须表示 10%-90% 区间。

1.5 碱基含量分布

碱基含量分布一般用于检测有无AT、GC分离现象。对于RNASeq来说，鉴于序列打断的随机性和G/C、A/T含量分别相等的原则，理论上每个测序循环中的GC含量相等、AT含量相等（如果是链特异性建库，可能会出现AT分离和/或GC分离），且在整个测序过程基本稳定不变，呈水平线。但在现有的高通量测序技术中，反转录合成 cDNA 时所用的6bp的随机引物会引起前几个位置的核苷酸组成存在一定的偏好性，这种波动属于正常情况。碱基含量分布结果见图2。

samples_R1samples_R2

图2 碱基含量分布图

注：横坐标是 Reads 中碱基位置(5’->3’)，纵坐标是该位点某碱基所占的比例统计。

1.6 Reads 平均质量分布

Reads 平均质量分布主要用来检测测序数据的平均质量分布情况。峰尖代表主体 Reads 测序质量，峰宽表示整体测序质量分布，较大峰宽或者前拖尾峰表示测序数据的部分数据质量偏低，峰尖对应值较低表明整体测序结果较差，如果没有峰尖，则表示 Reads 整体质量很好。Reads 平均质量分布结果见图3。

samples_R1samples_R2

图3 Reads 平均质量分布图

注：Reads质量分布图，横坐标表示 Reads 的平均质量，纵坐标为对应平均质量值的 Reads 数目。

2.1 比对结果基本统计

通过Bowtie2建立参考基因组索引，然后使用Bowtie2（http://bowtie-bio.sourceforge.net/index.shtml）将过滤后的Reads比对到参考基因组上。序列比对的基本信息统计见下表。

2.2 比对区域分布统计

比对到基因上不同结构区域的Reads数目，除了与样品有关，还与该物种基因组注释的情况有关。所使用的参考基因组注释的越全面，比对统计结果就越能反映样品实际情况。

samples

图4 比对结果统计（基因区域）

2.3 基因覆盖均一度

测序 Reads 在基因上覆盖度的分布情况，见图6。展示每个样品所有基因的5'到3'区域上序列覆盖情况，用于评估测序结果的均一性（或是否有偏向性）。理想条件下，Reads在所有表达的基因上的分布应该呈现均一化分布。

samples

图5 基因上测序深度分布

注：横坐标为单个基因的碱基长度占总碱基长度的百分比，0表示基因的5'端,100表示基因的3'端；纵坐标为比对到所有基因的横轴位置上相应区间内的序列条数的总和。图中体现了所有基因覆盖情况的叠加结果，曲线中每个点的纵坐标表示所有基因在该相对比例位置上所有序列的数量；曲线反映了测序所得序列是否在基因上均匀分布。若无明显偏向锋，则说明测序无偏向性。

1.1 表达量分析

我们使用HTSeq 0.6.1p2（http://www-huber.embl.de/users/anders/HTSeq）统计比对到每一个基因上Read Count值，作为基因的原始表达量。统计比对到每一个基因上Read Count值，作为基因的原始表达量。Reads计数与基因的真实表达水平，以及基因的长度和测序深度成正相关。为了使不同基因、不同样本间的基因表达水平具有可比性，我们采用FPKM对表达量进行标准化（Nosrmalization），FPKM（Fragments Per Kilo bases per Million fragments）是每百万片段中来自某一基因每千碱基长度的片段数目，对于 Pair-End 测序，每个 Fragments 会有两个 Reads，FPKM 只计算两个 Reads 能比对到同一个转录本的 Fragments 数量。在有参转录组当中，我们一般认为FPKM>1的基因是表达的。这个阈值是主流杂志推荐的，也能够很好的反应基因的表达水平。 统计方法如下：首先读取基因结构注释信息（GTF文件），然后将比对结果与基因结构进行比较并统计结果。HTSeq 有三种统计方案，如图1，其区别在于当一个 Read 仅有一部分覆盖在基因区域上或有一部分覆盖在基因的内含子区域时， Union 方案和 Intersection_nonempty 方案认定 Read 属于该基因，而 Intersection_strict 方案认定 Read 不属于任何基因；当一个 Read 全部覆盖在一个基因上，并且部分覆盖在另一个基因上时，Union 方案认定 Read 同属于两个基因，Intersection_strict 方案和 Intersection_nonempty 方案则认定 Read 属于第一个基因。如无特殊要求，应按 Union 方案统计，该方案较为稳健（如果是链特异性建库，则还需要判别是否和注释中的 Feature 方向一致）。

图1 Read Count统计模式图

根据表达量的计算结果表格，将表达量分为不同的区间，对各样本在不同表达量区间内的基因的数目进行统计，结果见下图。

samples

图2 mRNA表达量区间统计图

注：横坐标表示表达量FPKM值范围，纵坐标表示该表达量区间基因的个数。

根据表达量的计算结果表格，对各样本中鉴定到的mRNA进行统计，并对各样本间共有特有的mRNA统计结果进行upset图展示，结果见下图。

图3 每个样本鉴定得到的mRNA upset图

注：number in each set 表示每个样本鉴定到的全部mRNA的数目；number of each intersection表示多个样本鉴定到的共有mRNA的数目；横坐标一个点表示该样本鉴定到的特有mRNA的数目；横坐标多个点连线表示连线的多个样本鉴定到的共有mRNA的数目。

1.2 FPKM 密度分布

FPKM密度分布能整体考察样品的基因表达模式，一般来说中等表达的基因占绝大多数，低表达和高表达的基因占一小部分。

图4 FPKM密度分布

注：横坐标为基因的log10(FPKM)值，纵坐标为对应表达量的基因分布密度。

图5 FPKM密度分布与密度分布统计

注：中间的横线是中位数，盒型的上下边缘为75%，上下限为90%。外部形状为核密度估计。

1.3 饱和度分析

我们使用RSeQC分析测序饱和度，即：对测序结果进行抽样，对不同的抽样比例分别获得所有基因的RPKM值（标准化的基因表达量），然后与全部测序数据得到的RPKM值行比较，获得相对误差。如果随抽样比例的增加，相对误差减少，说明测序结果趋于饱和。

samples

图6 RPKM饱和图

注：该图为RPKM饱和图，横坐标是重采样的比率，纵坐标为相对误差。

1.4 样品相关性检验

样品间基因表达水平相关性是检验实验可靠性和样本选择是否合理的重要指标，在做差异表达分析之前，应先检查样品间基因的表达水平相关性。我们用皮尔逊相关系数表示样品间基因的表达水平相关性，相关系数越接近1，表明样品之间表达模式的相似度越高。一般来说，相关系数在0.8-1之间属于极强相关，如果生物学重复的样本之间相关系数低于0.8，表示样品之间的相关性较低，见图7。

图7 样品相关性检验

注：左侧和上侧为样本聚类情况，图中右侧和下侧为样本名称，不同颜色的方块代表两个样本的相关性高低情况。

1.5 PCA分析

PCA 主成分分析（Principal Components Analysis），通过线性变换，将高维数据降低至二维或三维，同时保持各方差贡献最大的特征，即降低数据复杂度。当有多个样品时，我们使用R语言的 DESeq软件包，根据表达量对各样品进行 PCA 主成分分析。PCA 分析可以把相似的样本聚到一起，距离越近表明样本间相似性越高。结果见图8， PCA图中Replication与sample的对应关系见分析结果中的PCA.xlsx。

图8 PCA分析

注：横坐标为第一主成分，纵坐标为第二主成分。在图中的不同形状表示不同的样品，不同的颜色表示不同的分组。

2.1 两样品差异表达检测

我们采用 DESeq 对基因表达进行差异分析，筛选差异表达基因条件为：表达差异倍数 |log2FoldChange| > 1 ，显著性P-value<0.05 。部分差异分析结果见表2，不同分组之间的差异表达基因统计结果见表3。

图9 表达差异分析结果统计

注：横坐标表示进行差异分析的比较组，纵坐标表示差异基因的个数，颜色中红色表示上调基因，绿色表示下调基因。

采用R语言 ggplots2 软件包绘制差异表达基因的火山图和 MA 图，见图10。火山图展示的是基因分布情况，基因的表达倍数差异和显著性结果，正常状况下，该图左右差异基因分布应大致对称，左侧为 Case 相比于 Control 下调基因，右侧为 Case 相比于 Control 上调基因。MA 图常用来展示标准化后的基因分布情况，一般标准化后，基因的表达量呈对称分布，即表达差异趋势不随基因表达量变化而发生偏向。

mRNA_de_pairs 火山图

图10 差异表达基因的火山图

注：火山图：横坐标为 log2FoldChange，纵坐标为 -log10（p-value）。图中两条竖虚线为2倍表达差异阈值；横虚线为 P-value=0.05 阈值。红点表示该组上调基因，蓝点表示该组下调基因，灰点表示非显著差异表达基因。

mRNA_de_pairs MA图

图11 差异表达基因MA图

注：横坐标为两样品基因表达量之积对2取对数的值，即log2(A*B)，A和B分别表示基因在两样本的表达量，纵坐标为表达量之商对2取对数的值，即 log2(A/B)。红点表示该组上调基因，蓝点表示该组下调基因，灰点表示非显著差异表达基因。

2.2 基因组圈图

我们使用R语言Circlize包，根据基因组信息和RNA差异表达分析结果，在基因组上标记差异表达的RNA。

图12 基因组圈图

注：最外圈是染色体条带，从外向内为不同差异分析的差异表达分析结果。红色和绿色分别为上调和下调基因的 log2FoldChange 值的柱状图，灰色为无差异表达基因的 log2FoldChange 值的散点图。

2.3 多组差异表达分析比较

根据差异分析的结果，统计各比较组之间的共有特有差异基因数量，如图13。维恩图（仅提供两组、 三组、 四组和五组比较的维恩图）展示了各比较组间差异基因的个数，以及比较组间的重叠关系，Upset图（仅提供两组以上比较的矩阵图）展示了两两比较组间共有的差异基因数目。

Upset图维恩图

图13 各组差异分析之间共有特有的差异基因

注：Upset图：number in each set 表示每个比较组鉴定到的全部差异基因的数目，number of each intersection表示多个比较组鉴定到的共有差异基因的数目，横坐标一个点表示该比较组鉴定到的特有差异基因的数目，横坐标多个点连线表示连线的多个比较组鉴定到的共有差异基因的数目；维恩图：各个圈中的数字之和代表该比较组合的差异基因总个数，圆圈交叠的部分表示两个比较组之间共有的差异基因。

2.4 聚类分析

聚类分析用于判断差异表达基因在不同实验条件下的表达模式；在样品间表达量相关性高的基因被归为一类，通常这些基因在某些生物学过程，或者某个代谢、信号通路中存在实际的联系。因此通过表达量聚类我们可以发现基因间未知的生物学联系。

我们使用R语言 Pheatmap 软件包对所有比较组的差异基因的并集和样品进行双向聚类分析，根据同一基因在不同样品中的表达水平和同一样品中不同基因的表达模式进行聚类，采用 Euclidean 方法计算距离，层次聚类最长距离法（Complete Linkage）进行聚类。结果见图14。

图14 差异表达基因聚类

注：横向表示基因，每一列为一个样本，红色表示高表达基因，绿色表示低表达基因。

2.5 趋势分析

根据上一节聚类分析中层次聚类的结果，将差异基因按照表达模式的不同，划分到不同的Cluster（同一个Cluster中的基因表达趋势相近），从而获得基因的聚类结果。结果展示见图15。

图15 趋势分析

注：图中灰色线展示每个Cluster中基因的表达模式，蓝色的线表示Cluster中的所有基因在样品中表达量的平均值。

2.6 蛋白网络互作

STRING (Search3 Tool for the Retrieval of Interacting Genes/Proteins)是EMBL开发的蛋白质互作数据库，该数据库从最有力的实验证据到数据挖掘、同源预测的蛋白质互作关系都有收录。 当STRING数据库中收录该物种的PPI信息时，我们根据基因差异表达分析结果，直接数据库中筛选双端节点皆为差异基因并且Score>0.95的PPI作用对。当STRING数据库中没有该物种的PPI信息时，我们选择相近物种与该物种的序列进行比对，进而得到该物种的蛋白质之间的相互关系。 分析结果见图17。

mRNA_de_pairs

图16 PPI网络图

注：一条连线代表一组相互关系；节点大小表示该节点的度。

3.1 GO富集分析

使用topGO进行GO富集分析，分析时利用GO term 注释的差异基因对每个term 的基因列表和基因数目进行计算，然后通过超几何分布方法计算P-value（显著富集的标准为P-value < 0.05），找出与整个基因组背景相比，差异基因显著富集的GO term ，从而确定差异基因行使的主要生物学功能。

对差异表达的基因的GO富集分析结果，按照分子功能MF、生物过程BP 和细胞组分CC 进行GO 分类，挑选每个GO 分类中挑选p-value最小即富集最显著的前10个GO term条目进行展示，结果见图18。

mRNA_de_pairsmRNA_de_pairs.UPmRNA_de_pairs.DOWN

图18 GO富集分析柱状图

注：横坐标为Go level2等级的term，纵坐标为每个term富集的-log10（p-value）。

根据GO富集结果，通过Rich factor、FDR值和富集到此GO Term 上的基因个数来衡量富集的程度。其中，Rich factor 指该GO Term 中富集到的差异基因个数与注释到的差异基因个数的比值。Rich factor越大，表示富集的程度越大。FDR一般取值范围为0-1，越接近于零，表示富集越显著。挑选FDR值最小的即富集最显著的前20个GO Term 条目进行展示，结果见图19。

mRNA_de_pairsmRNA_de_pairs.UPmRNA_de_pairs.DOWN

图19 GO富集分析气泡图

3.2 GO有向无环图

富集分析结果会分别给出GO三个 ontology（细胞组分、分子功能和生物过程）的有向无环图，在这个图中，越接近根结点的 GO term越概括，往下分支的 GO term为注释到更细层级的 term。 程序默认把显著性最高前10个GO term设置为方形，其他的GO term为圆形。颜色越深，代表该GO term越显著，颜色由浅到深为无色-浅黄-深黄-红色，见图20。

图20 Gene Ontology DAG富集图（示例）

注：每个节点代表一个 GO 术语，分支代表包含关系，从上至下所定义的功能范围越来越小，方框代表富集程度 top10 的 GO 术语，颜色越深代表富集程度越高。

3.3 KEGG富集分析

根据差异表达的基因的 KEGG 富集分析结果，挑选 p-value 最小即富集最显著的前20个 Pathway 进行展示，结果见图21。

mRNA_de_pairsmRNA_de_pairs.UPmRNA_de_pairs.DOWN

图21 KEGG Pathway 富集结果柱状图

注：横坐标为pathway名称，纵坐标为每个pathway富集的-log10（p-value）。

根据 KEGG 富集结果，通过 Rich factor、FDR 值和富集到此 pathway上的基因个数来衡量富集的程度。其中，Rich factor 指该pathway 中富集到的差异基因个数与注释到的差异基因个数的比值。Rich factor 越大，表示富集的程度越大。FDR 一般取值范围为0-1，越接近于零，表示富集越显著。挑选 FDR 值最小的即富集最显著的前20条 KEGG pathway 进行展示，结果见图22。

mRNA_de_pairsmRNA_de_pairs.UPmRNA_de_pairs.DOWN

图22 KEGG 富集分析

3.4 KO分析

ko 编号表示一个通路，这个通路是不分物种的，相当于所有物种的这一通路的并集。K 编号表示一个基因，是 ko 通路中的基本单位，某一 K 编号代表的不是某一具体物种的基因，而是所有物种的某一同源基因的统称。在 KEGG_enrichment.xls（具体位置见下方结果文件）中 URL 列为该行 Pathway 在 KEGG 数据库中对应的代谢通路图的链接，打开 URL 列的网页，可查看实验组和对照组差异表达基因在 Pathway 通路上的不同。下图为一个代谢通路图示例（仅为图片示例，不能演示 KEGG 网站中的动态效果），标红色的节点包含实验组上调基因；标绿色的节点包含对照组上调基因；标紫色的节点既包含实验组上调基因，也包含对照组上调基因。将鼠标指在节点上将会弹出节点的详细描述，点击各个节点可以跳转到 KEGG 数据库中对应的具体信息页。右键点击图片可选择保存该代谢通路图。

图23 代谢通路图示例

注：方框一般表示酶，小圆圈表示代谢物，圆角框表示另一个代谢通路图。

4.1 cSNP 和 InDel分析

cSNP 包括转换和颠换两种。SNP 在 CG 序列上出现的最为频繁，而且多是 C 转换为 T，原因是 CG 中的 C 常为甲基化的，自发地脱氨后即为胸腺嘧啶。SNP 出现的原因有很多，有的是遗传背景的单核苷酸多态性，有的是建库技术上造成的突变，有的则可能是测序中读取错误。InDel可能造成移码突变，导致 mRNA 翻译时遇上一个错误的终止密码子。一般 InDel 不是3的倍数的情况在编码区不常发生，在非编码区相对频繁地发生。除了在高度重复区域附近，InDel 发生的频率一般会低于 SNP。

我们采用 Varscan 程序获取 SNP 和 InDel 位点，过滤标准为：

• 1）SNP 位点碱基 Q >20；
• 2）覆盖该位点的 Reads 数目> 8 ；
• 3）支持突变位点的 Reads 数目> 2；
• 4）SNP 位点的 p-value < 0.01。

SNP 分析结果见表6，InDel 分析结果格式同 SNP 分析结果，SNP、InDel 和转换/颠换统计见图24。

SNP 统计InDel 统计转换/颠换统计

图24 SNP 和 InDel 统计

注：Homo（homozygous-variant） 表示纯合子变异体，即这一位点的等位基因都突变了，并且突变相同，Hete（heterozygous-variant）表示杂合子变异体，即这一位点的等位基因至少有一个突变了，且突变后等位基因不同。

4.2 操纵子结构分析

完整的操纵子（Operon）包括操纵子基因以及与之相关的调节基因（Regulator），如下图是典型的乳糖操纵子。我们通过Rockhopper预测操纵子转录起始位点（Transcription Start Site, TSS）、转录终止位点（Transcription Termination Site, TTS）。

4.3 UTR

我们根据基因转录起始位点（转录终止位点）和翻译起始位点（翻译终止位点）信息，提取5'UTR（3'UTR）序列，并对其长度分布情况进行统计。

针对5'UTR，用ELFH（http://ccb.jhu.edu/software/ELPH/index.shtml）软件对预测基因上游25bp序列按6bp窗口进行分析，查找可能存在的RBS序列（ribosomebinding site，核糖体结合位点，也称为SD序列）。

针对3'UTR，用TransTermHP（http://transterm.cbcb.umd.edu/）软件对不依赖σ因子的终止子进行预测。