返回顶部

项目概况

1 项目概况

2 实验流程

3 参考基因组

我们在分析有参考基因组的转录组数据时,需要将测序数据与参考基因组进行比对,本次项目使用的基因组信息如下表所示。

表1 基因组组装统计

将数据库中获得的基因组信息整理如下表所示:

表2 参考基因组注释信息统计

从不同的数据库中整理该物种的注释信息包括:1)基因的编号(ENSEMBL ID);2)染色体上位点信息(所在染色体序列编号、起始终止位点、长度、转录方向);3)基因的命名(HGVS Symbol、NCBI Gene ID、UniProtKB ID);4)基因的分类(GO、KO、EC分类);5)基因的文字注释信息。基因注释情况见下表。

表3 参考基因组基因注释信息统计

4 原核分析流程

首先对原始下机数据进行过滤以获得高质量序列,将过滤后的序列比对到该物种的参考基因组上。根据比对结果,计算每个基因的表达量。在此基础上,进一步对样本进行表达差异分析、富集分析和聚类分析;同时,对样本的转录本结构进行分析;获得样本基因的操纵子结构,UTR、cSNP和InDel等信息;对数据库中尚无注释的表达区域进行注释;最后将结果进行可视化展示。

基本分析

1 原始数据整理、过滤及质量评估

1.1 文库基本情况

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

表1 文库基本情况

注:Sample:样品名称;
Lib. Name:文库名;
Lib. Insert Size:文库插入片段长度;
Sequencing platform:测序平台;
Sequencing Mode:测序模式。

1.2 数据整理

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

表2 下机数据统计

注:Sample:样品名;
Reads No.:Reads总数;
Bases(bp): 碱基总数;
Q30 (bp):碱基识别准确率在99.9%以上的碱基总数;
N(%):模糊碱基所占百分比;
Q20(%):碱基识别准确率在99%以上的碱基所占百分比;
Q30(%):碱基识别准确率在99.9%以上的碱基所占百分比。

samples Read1 FASTQ文件:./result/01_Rawdata/samples_R1.fastq.gz

samples Read2 FASTQ文件:./result/01_Rawdata/samples_R2.fastq.gz

1.3 数据过滤

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

  • 1)采用 Cutadapt 去除 3' 端带接头的序列;
  • 2)去除平均质量分数低于 Q20 的 Reads。
表3 数据过滤统计

注:Sample:样品名;
Clean Reads No:高质量序列read数;
Clean Data (bp):高质量序列碱基数;
Clean Reads %: 高质量序列 reads 占测序 reads 的百分比;
Clean Data %: 高质量序列碱基占测序碱基的百分比。

1.4 碱基质量分布

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

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

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

图1 单碱基质量分布图

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

samples Read1 单碱基质量分布图:./result/02_FastQC/samples_R1_fastqc/per_base_quality.png

samples Read2 单碱基质量分布图:./result/02_FastQC/samples_R2_fastqc/per_base_quality.png

1.5 碱基含量分布

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

图2 碱基含量分布图

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

1.6 Reads 平均质量分布

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

图3 Reads 平均质量分布图

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

2 比对分析

2.1 比对结果基本统计

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

表4 RNASeq Map 统计

注:Sample:样品;
Useful Reads:用于比对的序列总数;
Total Mapped Reads(%):比对上参考基因组的序列总数及比例;
Uniquely Mapped Reads(%):只比对到一个位置的序列总数及比例;
Multiple Mapped Reads(%):比对到多个位置的序列总数及比例。

2.2 比对区域分布统计

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

表5 RNASeq Mapped Events

注:Sample:样品名;
Total Mapped Reads:比对上参考基因组的序列总数;
InterGene(%):比对到基因间区域的Reads总数及比例;
Gene(%):比对到基因区域的Reads总数及比例;
mRNA(%):比对到蛋白编码基因上的Reads总数及比例;
rRNA(%):比对到rRNA基因间上的Reads总数及比例。

图4 比对结果统计(基因区域)

samples Reads 在基因区域分布图(png格式):./images/Read_Distribution/samples.read_distribution.png

samples Reads 在基因区域分布图:./result/04_MapQC/read_distribution/samples.read_distribution.pdf

2.3 基因覆盖均一度

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

图5 基因上测序深度分布

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

samples 基因上测序深度分布图(png格式):./images/Gene_Body_Coverage/samples.geneBodyCoverage.curves.png

samples 基因上测序深度分布图:./result/04_MapQC/genebody_coverage/samples.geneBodyCoverage.curves.pdf

mRNA

1 表达量分析

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统计模式图
表1 基因表达量(部分)

注:Gene_ID:基因ID
*read count:样本*的基因表达量(Reads数);
*.fpkm:样本*的基因表达量(FPKM值)。

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

图2 mRNA表达量区间统计图

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

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

单样品无此图图3 每个样本鉴定得到的mRNA upset图

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

samples 表达量区间统计图(png格式):./images/mRNA/1_Expression/samples.fpkm_distribution.png

样本间UpSet图(png格式):./images/mRNA/1_Expression/UpSet.png

所有样品的Reads count 和 FPKM:./result/mRNA/1_Expression/Expression.xlsx

样本间UpSet图:./result/mRNA/1_Expression/UpSet.pdf

samples 表达量区间统计图:./result/mRNA/1_Expression/samples.fpkm_distribution.pdf

样品表达量区间分布统计表:./result/mRNA/1_Expression/expression.rangecount.xlsx

1.2 FPKM 密度分布

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

图4 FPKM密度分布

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

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

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

FPKM 密度分布图(png格式):./images/mRNA/1_Expression/density.png

FPKM 密度分布小提琴图(png格式):./images/mRNA/1_Expression/violin_plot.png

FPKM 密度分布小提琴图:./result/mRNA/1_Expression/violin_plot.pdf

FPKM 密度分布图:./result/mRNA/1_Expression/density.pdf

1.3 饱和度分析

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

图6 RPKM饱和图

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

samples RPKM 饱和度图(png格式):./images/Saturation/samples.saturation.png

samples RPKM 饱和度图:./result/04_MapQC/RPKM_saturation/samples.saturation.pdf

1.4 样品相关性检验

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

图7 样品相关性检验

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

样品相关性检验图(png格式):./images/mRNA/1_Expression/cor_pearson.png

样品相关性检验图:./result/mRNA/1_Expression/cor_pearson.pdf

样品相关性表:./result/mRNA/1_Expression/cor_pearson.xlsx

1.5 PCA分析

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

图8 PCA分析

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

2 表达差异分析

2.1 两样品差异表达检测

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

表2 表达差异分析结果(部分)

注:id:基因编号;
baseMean:表示进行差异比较的两组全部样品该基因read count的均一化结果;
baseMean(sample):表示该组全部样品该基因read count的均一化结果;
foldChange(Case/Control):表达差异倍数;
log2FoldChange:表达差异倍数对2取对数的值;
pval:显著性 p-value;
padj:校正后的显著性 p-value。

表3 表达差异分析结果统计

注:Case:实验组样本;
Control:对照组样本;
Up-regulated Genes:Case 相比于Control 上调表达基因
Down-regulated Genes:Case 相比于Control 下调表达的基因
Total DEGs:Case 相比于Control 差异表达基因

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

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

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

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

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

图11 差异表达基因MA图

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

表达差异分析结果统计图(png格式):./images/mRNA/2_DEG/DEG_stat.png

mRNA_de_pairs 火山图(png格式):./images/mRNA/2_DEG/mRNA_de_pairs.Volcano.png

mRNA_de_pairs MA图(png格式):./images/mRNA/2_DEG/mRNA_de_pairs.MA.png

差异分析结果统计:./result/mRNA/2_DEG/DESeq/DEG_stat.xlsx

差异分析结果统计图:./result/mRNA/2_DEG/DESeq/DEG_stat.pdf

mRNA_de_pairs 差异分析结果:./result/mRNA/2_DEG/DESeq/mRNA_de_pairs/mRNA_de_pairs.DESeq.xlsx

mRNA_de_pairs 中实验组下调基因:./result/mRNA/2_DEG/DESeq/mRNA_de_pairs/mRNA_de_pairs.DESeq.Down.xlsx

mRNA_de_pairs 中实验组上调基因:./result/mRNA/2_DEG/DESeq/mRNA_de_pairs/mRNA_de_pairs.DESeq.Up.xlsx

mRNA_de_pairs 火山图:./result/mRNA/2_DEG/V_MA/mRNA_de_pairs.Volcano.pdf

mRNA_de_pairs MA图:./result/mRNA/2_DEG/V_MA/mRNA_de_pairs.MA.pdf

2.2 基因组圈图

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

染色体过多无法画此图图12 基因组圈图

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

2.3 多组差异表达分析比较

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

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

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

各组差异分析之间共有特有的差异基因UpSet图(png格式):./images/mRNA/2_DEG/UpSet.png

各组差异分析之间共有特有差异基因维恩图(png格式):./images/mRNA/2_DEG/venn.png

各组差异分析之间共有特有的差异基因维恩图:./result/mRNA/2_DEG/Venn/venn.pdf

各组差异分析之间共有特有的差异基因Upset图:./result/mRNA/2_DEG/Venn/UpSet.pdf

各组差异分析之间共有特有的差异基因列表:./result/mRNA/2_DEG/Venn/diff_gene_of_all_groups.xlsx

2.4 聚类分析

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

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

图14 差异表达基因聚类

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

差异表达基因聚类(png格式):./images/mRNA/2_DEG/ALL.png

mRNA_de_pairs差异表达基因聚类图:./result/mRNA/2_DEG/Heatmap/mRNA_de_pairs.pdf

差异表达基因聚类表:./result/mRNA/2_DEG/Heatmap/ALL.heatmap.xlsx

差异表达基因聚类图:./result/mRNA/2_DEG/Heatmap/ALL.pdf

mRNA_de_pairs差异表达基因聚类表:./result/mRNA/2_DEG/Heatmap/mRNA_de_pairs.heatmap.xlsx

2.5 趋势分析

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

图15 趋势分析

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

趋势分析图(png格式):./images/mRNA/2_DEG/ALL_cluster_plots.png

2.6 蛋白网络互作

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

图16 PPI网络图

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

mRNA_de_pairs PPI网络图(png格式):./images/mRNA/2_DEG/mRNA_de_pairs.network.png

趋势分析图:./result/mRNA/2_DEG/trend/ALL_cluster_plots.pdf

趋势分析结果表:./result/mRNA/2_DEG/trend/ALL_cluster.xlsx

mRNA_de_pairs PPI网络图:./result/mRNA/2_DEG/PPI/mRNA_de_pairs.network.pdf

mRNA_de_pairs 节点信息表:./result/mRNA/2_DEG/PPI/mRNA_de_pairs.attributes.xlsx

mRNA_de_pairs 互作关系对信息表:./result/mRNA/2_DEG/PPI/mRNA_de_pairs.network.xlsx

3 差异表达基因功能富集分析

3.1 GO富集分析

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

表4 GO富集分析(部分)

注:Catergory:富集到的GO Term所处的分类;
GO_Term:富集到的GO条目;
Up/Down:富集到该GO条目的上/下调基因;
Pvalue:富集显著性P值;
FDR:P值校正值。

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

图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_pairs差异基因GO富集分析柱状图(png格式):./images/mRNA/3_Enrichment/mRNA_de_pairs_GO_enrichment_pvalue_barplot.png

mRNA_de_pairs差异基因GO富集气泡图(png格式):./images/mRNA/3_Enrichment/mRNA_de_pairs_GO_richfactor.png

mRNA_de_pairs上调基因GO富集分析柱状图(png格式):./images/mRNA/3_Enrichment/mRNA_de_pairs_GO_enrichment_pvalue_barplot_UP.png

mRNA_de_pairs下调基因GO富集分析柱状图(png格式):./images/mRNA/3_Enrichment/mRNA_de_pairs_GO_enrichment_pvalue_barplot_DOWN.png

mRNA_de_pairs上调基因GO富集气泡图(png格式):./images/mRNA/3_Enrichment/mRNA_de_pairs_GO_richfactor_UP.png

mRNA_de_pairs下调基因GO富集气泡图(png格式):./images/mRNA/3_Enrichment/mRNA_de_pairs_GO_richfactor_DOWN.png

mRNA_de_pairs 差异基因 GO 富集分析结果:./result/mRNA/3_Enrichment/mRNA_de_pairs/GO_enrichment.xlsx

mRNA_de_pairs差异基因GO富集气泡图:./result/mRNA/3_Enrichment/mRNA_de_pairs/GO.richfactor.pdf

mRNA_de_pairs差异基因GO富集柱状图:./result/mRNA/3_Enrichment/mRNA_de_pairs/GO_enrichment_pvalue_barplot.pdf

mRNA_de_pairs上调基因GO富集分析柱状图:./result/mRNA/3_Enrichment/mRNA_de_pairs.UP/GO_enrichment_pvalue_barplot.pdf

mRNA_de_pairs上调基因GO富集分析结果:./result/mRNA/3_Enrichment/mRNA_de_pairs.UP/GO_enrichment.xlsx

mRNA_de_pairs上调基因GO富集气泡图:./result/mRNA/3_Enrichment/mRNA_de_pairs.UP/GO.richfactor.pdf

mRNA_de_pairs下调基因GO富集气泡图:./result/mRNA/3_Enrichment/mRNA_de_pairs.DOWN/GO.richfactor.pdf

mRNA_de_pairs下调基因GO富集分析柱状图:./result/mRNA/3_Enrichment/mRNA_de_pairs.DOWN/GO_enrichment_pvalue_barplot.pdf

mRNA_de_pairs下调基因GO富集分析结果:./result/mRNA/3_Enrichment/mRNA_de_pairs.DOWN/GO_enrichment.xls

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 术语,颜色越深代表富集程度越高。

mRNA_de_pairs 差异基因在生物学过程分类下的 GO 富集 DAG图:./result/mRNA/3_Enrichment/mRNA_de_pairs/topGO_BP_top10.pdf

mRNA_de_pairs 差异基因在细胞的组件作用分类下的 GO 富集 DAG图:./result/mRNA/3_Enrichment/mRNA_de_pairs/topGO_CC_top10.pdf

mRNA_de_pairs 差异基因在基因的分子功能分类下的 GO 富集 DAG图:./result/mRNA/3_Enrichment/mRNA_de_pairs/topGO_MF_top10.pdf

mRNA_de_pairs 上调基因在生物学过程分类下的 GO 富集 DAG图:./result/mRNA/3_Enrichment/mRNA_de_pairs.UP/topGO_BP_top10.pdf

mRNA_de_pairs 上调基因在细胞的组件作用分类下的 GO 富集 DAG图:./result/mRNA/3_Enrichment/mRNA_de_pairs.UP/topGO_CC_top10.pdf

mRNA_de_pairs 上调基因在基因的分子功能分类下的 GO 富集 DAG图:./result/mRNA/3_Enrichment/mRNA_de_pairs.UP/topGO_MF_top10.pdf

mRNA_de_pairs 下调基因在生物学过程分类下的 GO 富集 DAG图:./result/mRNA/3_Enrichment/mRNA_de_pairs.DOWN/topGO_BP_top10.pdf

mRNA_de_pairs 下调基因在细胞的组件作用分类下的 GO 富集 DAG图:./result/mRNA/3_Enrichment/mRNA_de_pairs.DOWN/topGO_CC_top10.pdf

mRNA_de_pairs 下调基因在基因的分子功能分类下的 GO 富集 DAG图:./result/mRNA/3_Enrichment/mRNA_de_pairs.DOWN/topGO_MF_top10.pdf

3.3 KEGG富集分析

表5 KEGG富集分析(部分)

注:PathwayID:通路ID;
Pathway:通路名称;
Up/Down_number:富集到该条通路下的上/下调基因数量;
DEG_number:富集到该条通路下的差异基因总数;
total_number:该条通路下的基因总数;
Up/Down_gene:富集到该条通路下的上/下调基因ID(名称)。

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

图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_pairs差异基因KEGG富集分析柱状图(png格式):./images/mRNA/3_Enrichment/mRNA_de_pairs_KEGG_enrichment_pvalue_barplot.png

mRNA_de_pairs 差异基因 KEGG 富集气泡图(png格式):./images/mRNA/3_Enrichment/mRNA_de_pairs_KEGG_richfactor.png

mRNA_de_pairs上调基因KEGG富集分析柱状图(png格式):./images/mRNA/3_Enrichment/mRNA_de_pairs_KEGG_enrichment_pvalue_barplot_UP.png

mRNA_de_pairs上调基因KEGG富集气泡图(png格式):./images/mRNA/3_Enrichment/mRNA_de_pairs_KEGG_richfactor_UP.png

下调基因KEGG富集分析柱状图(png格式):./images/mRNA/3_Enrichment/mRNA_de_pairs_KEGG_enrichment_pvalue_barplot_DOWN.png

下调基因KEGG富集气泡图(png格式):./images/mRNA/3_Enrichment/mRNA_de_pairs_KEGG_richfactor_DOWN.png

mRNA_de_pairs 差异基因 KEGG 富集分析结果:./result/mRNA/3_Enrichment/mRNA_de_pairs/KEGG_enrichment.xlsx

mRNA_de_pairs差异基因KEGG富集气泡图:./result/mRNA/3_Enrichment/mRNA_de_pairs/KEGG.richfactor.pdf

mRNA_de_pairs差异基因KEGG富集柱状图:./result/mRNA/3_Enrichment/mRNA_de_pairs/KEGG_enrichment_pvalue_barplot.pdf

mRNA_de_pairs上调基因KEGG富集分析柱状图:./result/mRNA/3_Enrichment/mRNA_de_pairs.UP/KEGG_enrichment_pvalue_barplot.pdf

mRNA_de_pairs上调基因KEGG富集分析结果:./result/mRNA/3_Enrichment/mRNA_de_pairs.UP/KEGG_enrichment.xls

mRNA_de_pairs上调基因KEGG富集气泡图:./result/mRNA/3_Enrichment/mRNA_de_pairs.UP/KEGG.richfactor.pdf

mRNA_de_pairs下调基因KEGG富集分析柱状图:./result/mRNA/3_Enrichment/mRNA_de_pairs.DOWN/KEGG_enrichment_pvalue_barplot.pdf

mRNA_de_pairs下调基因KEGG富集分析结果:./result/mRNA/3_Enrichment/mRNA_de_pairs.DOWN/KEGG_enrichment.xls

mRNA_de_pairs下调基因KEGG富集气泡图:./result/mRNA/3_Enrichment/mRNA_de_pairs.DOWN/KEGG.richfactor.pdf

3.4 KO分析

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

图23 代谢通路图示例

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

mRNA_de_pairsKEGG代谢通路结果网页版:./result/mRNA/3_Enrichment/mRNA_de_pairs/KEGG_enrichment_map.html

mRNA_de_pairsKEGG Pathway通路图:./result/mRNA/3_Enrichment/mRNA_de_pairs/KEGG_Pathway_map

mRNA_de_pairs上调基因KEGG代谢通路结果网页版:./result/mRNA/3_Enrichment/mRNA_de_pairs.UP/KEGG_enrichment_map.html

mRNA_de_pairs上调基因KEGG Pathway通路图:./result/mRNA/3_Enrichment/mRNA_de_pairs.UP/KEGG_Pathway_map

mRNA_de_pairs下调基因KEGG代谢通路结果网页版:./result/mRNA/3_Enrichment/mRNA_de_pairs.DOWN/KEGG_enrichment_map.html

mRNA_de_pairs下调基因KEGG Pathway通路图:./result/mRNA/3_Enrichment/mRNA_de_pairs.DOWN/KEGG_Pathway_map

4 结构分析

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

表6 SNP/InDel 分析结果(部分)

注:CHROM:SNP位点所在染色体;
POS:SNP位点在染色体上的位置;
REF/ALT:参考序列在该位点的基因型/突变基因型;
样品名下的列:支持各基因型的 Reads 数目。

图24 SNP 和 InDel 统计

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

4.2 操纵子结构分析

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

表7 操纵子预测结果(部分)

注:Start:第一个基因起始位置;
Stop.:最后一个基因的终止位置;
Strand.:方向;
Number of Genes:基因数量;
Genes:基因列表。

表8 TSS和TTS预测结果(部分)

注:TSS:转录起始位点;
TSS:转录终止位点;
Strand.:方向;
Gene name:基因名;
Product:基因产物。

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/)软件对不依赖σ因子的终止子进行预测。

表9 UTR信息(部分)

注:Chromosome:染色体编号,;
Start:UTR起始位置;
End:UTR终止位置;
Name:UTR类型|基因ID|基因名;
Length:UTR长度。

表10 RBS序列预测结果(部分)

注:Seq.no:序列顺序;
Pos:在染色体上的位置;
*****,Motif,*****:RBS序列(及邻近的序列);
Prob:概率;
Seq.Id:序列名。

表11 UTR不依赖σ因子的终止子预测(部分)

注:Gene:终止子所在基因;
TermStart:终止子起始位点;
Term End:终止子起始位点;
Strand:方向;
5'tail:茎环结构的 5’尾序列;
5'stem-loop-3'stem:茎环结构的 5’茎-环-3’茎序列;
3'tail:茎环结构的 3’尾序列。

RBS序列预测结果:./result/mRNA/4_Structure/UTR/ELPH_report.xlsx

UTR分析结果:./result/mRNA/4_Structure/UTR/UTR.xlsx

3’UTR不依赖σ因子的终止子预测:./result/mRNA/4_Structure/UTR/terminators.xlsx

附录

1 数据库介绍

GO

基因本体论联合会建立的数据库(Gene Ontology,http://geneontology.org/)。GO的产生主要是为了解决同一基因在不同数据库定义的混乱性以及不同物种的同一基因在功能定义上的混乱性。它是一个国际标准化的基因功能分类体系,提供了一套动态更新的标准词汇表(Controlled Vocabulary)来全面描述生物体中基因和基因产物的属性。GO 涵盖三个方面,分别描述基因的分子功能(Molecular Function)、细胞的组件作用(Cellular Component)、参与的生物学过程(Biological Process)。基因或蛋白质可以通过 ID 对应或者序列注释的方法找到与之对应的 GO 编号,而 GO 编号可用于对应到 GO Term,即功能类别或者细胞定位。

GO 的基本单元是 Term,每个 Term 有一个唯一的标示符(由 “GO:” 加上7个数字组成,例如 GO:0072669);每类 Ontology 的 Term 通过它们之间的联系( is_a, part_of, regulate)构成一个有向无环的拓扑结构。GOSlim 是缩减版的 GO 术语,它提供了 GO 注释的概述性结果。

京都基因与基因组百科全书(Kyoto Encyclopedia of Genes and Genomes,http://www.kegg.jp/)是一个整合了基因组、化学和系统功能信息的数据库。把从已经完整测序的基因组中得到的基因目录与更高级别的细胞、物种和生态系统水平的系统功能关联起来是KEGG数据库的特色之一。KEGG 注释主要包括:(1)KO (KEGG Ortholog)注释,即将分子网络的相关信息进行跨物种注释; (2)KEGG Pathway 注释,即代谢通路注释,获得物种内分子间相互作用和反应的网络。

UniProt 知识库(UniProt Knowledgebase,http://www.uniprot.org/help/uniprotkb)的子数据库,是经过有经验的分子生物学家和蛋白质化学家仔细核实的高质量的、手工注释的、非冗余的蛋白数据集。SwissProt数据库的每个条目都有详细的注释,包括结构域、功能位点、跨膜区域、二硫键位置、翻译后修饰、突变体等。该数据库中还包括了与核酸序列数据库EMBL/GenBank/DDBJ、蛋白质结构数据库PDB以及Prosite、PRINTTS等十多个二次数据库的交叉引用代码。

EC

国际生物化学会酶学委员会(Enzyme Commission,http://enzyme.expasy.org/),根据酶所催化的反应类型和机理,把酶分成6大类:氧化还原酶、转移酶、水解酶、裂合酶、异构酶及合成酶。

真核生物直系同源蛋白质聚类(Evolutionary Genealogy of Genes: Non-supervised Orthologous Groups,http://eggnog.embl.de/version_3.0/),具体信息请参考(http://www.ncbi.nlm.nih.gov/COG/ )。我们将列出所有基因的 eggNOG ID ,然后把这些 eggNOG ID 归入适当的 eggNOG 分类单元(Category),由此对基因组的所有基因功能做分类统计,从宏观上认识该物种的基因功能分布特征。

构成每个 eggNOG 的蛋白都是被假定为来自于一个祖先蛋白,并且因此或者是直系同源(Orthologs) 或者是旁系同源(Paralogs)。Orthologs 是指来自于不同物种的由垂直家系(物种形成)进化而来的蛋白,并且典型的保留与原始蛋白有相同的功能。Paralogs 是那些在一定物种中的来源于基因复制的蛋白,可能会进化出新的与原来有关的功能。

2 所用软件介绍

表1 所用软件介绍
软件功能参数
Cutadapt数据过滤至少 10 bp Overlap(AGATCGGAAG),允许 20% 的碱基错误率
FastQC质量控制默认参数
RSeQCRPKM 饱和度分析默认参数
HTSeq表达定量使用union方案
DESeq差异分析|log2foldchang|>1 和 pvalue<0.05
ggplot2绘制火山图、MA图默认参数
Pheatmap聚类分析默认参数
topGOGO映射到DAG默认参数
Bowtie2Reads 与转录本/基因序列比对默认参数

3 参考文献

[1] Radakovits R, Jinkerson RE, Fuerstenberg SI, Tae H, Settlage RE, Boore JL, Posewitz MC. Draft genome sequence and genetic transformation of the oleaginous alga: Nannochloropsis gaditana[J]. Nature Communications. 2012,3: 686.
[2] Tatusov RL, Fedorova ND, Jackson JD, Jacobs AR, Kiryutin B, Koonin EV, Krylov DM, Mazumder R, Mekhedov SL, Nikolskaya AN, Rao BS, Smirnov S, Sverdlov AV, Vasudevan S, Wolf YI, Yin JJ, Natale DA. The COG database: an updated version includes eukaryotes[J]. BMC Bioinformatics. 2003 Sep 11,4:41.
[3] The Gene Ontology Consortium, Michael Ashburner, Catherine A. Ball, Judith A. Blake, David Botstein, Heather Butler, J. Michael Cherry, Allan P. Davis, Kara Dolinski, Selina S. Dwight, Janan T. Eppig, Midori A. Harris, David P. Hill, Laurie Issel-Tarver, Andrew Kasarskis, Suzanna Lewis, John C. Matese, Joel E. Richardson, Martin Ringwald, Gerald M. Rubin, and Gavin Sherlock. Gene Ontology: tool for the unification of biology[J]. Nat Genet. 2000 May, 25(1): 25–29.
[4] Minoru Kanehisa,* Susumu Goto, Shuichi Kawashima, Yasushi Okuno, and Masahiro Hattori. The KEGG resource for deciphering the genome[J]. Nucleic Acids Res. 2004 January 1; 32(Database issue): D277–D280.
[5] Powell S, Szklarczyk D, Trachana K, Roth A, Kuhn M, Muller J, Arnold R, Rattei T, Letunic I, Doerks T, Jensen LJ, von Mering C, Bork P. eggNOG v3.0: orthologous groups covering 1133 organisms at 41 different taxonomic ranges[J]. Nucleic Acids Res. Epub 2011 Nov 16; PubMed 22096231.
[6] Ryan McClure, Divya Balasubramanian, Yan Sun, Maksym Bobrovskyy, Paul Sumby, Caroline A. Genco, Carin K. Vanderpoo and Brian Tjaden. Computational analysis of bacterial RNA-Seq data[J]. Nucleic Acids Res. Aug 2013; 41(14): e140.
[7] T.H. Chang, L.C Wu, C.T. Yeh, B.J. Liu, H.D. Huang and J.T. Horng. Computational Identification of Riboswitches Based on RNA Conserved Functional Sequences and Conformations[J]. RNA, 2009, 15: 1426-1430.
[8] C. Kingsford, K. Ayanbule and S.L. Salzberg. Rapid, accurate, computational discovery of Rho-independent transcription terminators illuminates their relationship to DNA uptake. Genome Biology, 2007, 8:R22.

4 名词解释

  • 模糊碱基/ N:测序中不能确定的碱基,以N表示。一条序列中 N 越多说明该序列质量越低,一般该种序列需要剔除掉。
  • 接头/ Adapter:接头是测序时在序列两端分别加上的一段人工序列,接头上含有与测序引物互补结合的序列,通过和测序引物结合来对目的片段进行测序。当加上接头后的序列片段比实际测序读长短时, 3’端会测到接头序列,接头序列在分析之前需要去除掉。
  • log2FoldChange:表达倍数差异,同一基因在两个样品中的表达量之商对2取对数,即 log2(sampleA/sampleB)。
  • P-value:显著性,统计学根据显著性检验方法所得到的P 值,一般以P < 0.05 为显著, P <0.01 为非常显著,其含义是样本间的差异由抽样误差所致的概率小于0.05 或0.01。
  • cSNP:SNP(Single Nucleotide Polymorphisms)是指在基因组上由单个核苷酸变异形成的遗传标记,其数量很多,多态性丰富。cSNP 是指在编码区出现的 SNP,这些 SNP 直接影响到氨基酸密码子。
  • SNP 转换:嘧啶变成嘧啶或嘌呤变成嘌呤,即 A、G互换,T、C互换。
  • SNP 颠换:嘧啶突变成嘌呤或者相反,即 A、T互换,A、C互换,G、T互换,G、C互换。
  • InDel:Insertion-Deletion,指相对于参考基因组,样本中发生的小片段的插入或者缺失,该插入缺失可能含有一个或多个碱基。InDel 可作为一种基因标记用于研究系统进化或物种鉴定。
  • Read count:比对到一个基因上的 Reads 数目。
  • Read / Reads:测序中每一条序列称为一个 Read。
  • Raw Data / Raw Reads:测序下机的原始数据。
  • Clean Data / Clean Reads:去除接头和低质量 Reads 后的数据,后续分析均基于Clean Data。

5 常用术语

FASTQ 格式(http://en.wikipedia.org/wiki/FASTQ_format)是一种文本格式,常用于存储生物学序列及其对应的质量分值。FASTQ 格式文件可以采用文本编辑软件(如写字板、UltraEdit 、EditPlus等工具)打开,由于文件较大,对电脑的内存要求较高。FASTQ 格式中,每个 Read 由四行信息表示。

第一行为序列名称,以 @ 开头,其后是序列描述;第二行为碱基序列;第三行为 “+” 号,不代表任何意义;第四行碱基质量,与第二行的碱基序列一一对应。示例如下:

  • @M00200:111:000000000-A6VNV:1:1101:15594:1337 1:N:06
  • ACGCGGGTATCTAATCCTGTTTGCTCCCCACGCTTTCGCGCCTCAGTGTCAGTTAC
  • +
  • ABABADBBDFFFGGGFGGGFGGHGBGHGGHGGGGGGHGGGGGGGHHGGFBGEGGEG

我们使用 Sanger 质量值来评估下机数据的测序质量。质量值,简称 Q 值,是碱基读取错误率 p 的取整映射结果,等于 Phred quality score,计算公式为:

Qphred   =   -10 log10  p

Phred quality score计算公式

测序错误率与 Q 值的简明对应方式如表2所示。

表2 错误率与Q值对应关系
测序错误率Q 值
5%13
1%20
0.1%30
0.01%40

不同的测序平台,采用不同的方案对FASTQ文件中的碱基进行质量编码, Q 值与碱基质量的对应关系为:Q 值加上一个偏移数值,得到的结果按照 ASCII 码值对照表(见表3)转换成对应的字符,参考信息如下所示:

我们的FASTQ文件采用 Illumina 1.8+ 版本编码,将所有字符的 ASCII 值减去偏移值 33,即可得到碱基的 Q 值。例如,字符 I 的ASCII值为73,减去33 后得到 40,那么该字符对应位置的碱基质量为40,测序错误率则为0.01%。

表3 ASCII码表
十进制字符Q值十进制字符Q值十进制字符Q值十进制字符Q值
32  4801564@3180P47
3304911665A3281Q48
3415021766B3382R49
35#25131867C3483S50
36$35241968D3584T51
37%45352069E3685U52
38&55462170F3786V53
3965572271G3887W54
40(75682372H3988X55
4185792473I4089Y56
42*958:2574J4190Z57
43+1059:2675K4291[58
441160<2776L439259
45-1261=2877M4493]60
46.1362>2978N4594^61
471463?3079O4695_62

四分位数是指把所有数值由小到大排列并分成四等份,处于第一和第三个分割点位置的数值就是四分位数。

四分位数示例:

Sam(sequence alignment/map format)是一种由 Sanger 制定的序列比对格式标准,以 Tab 为分割符的文本格式,可用文本编辑软件打开(如写字板、UltraEdit 、EditPlus等工具),主要应用于测序序列比对到基因组上的结果表示,当然也可以表示任意的多重比对结果。当把 fastq 文件比对到基因组上之后,我们通常会得到一个 Sam 或者 Bam 为扩展名的文件。而 Bam 就是 Sam 的二进制文件(B 取自 binary),占用空间更小,不可打开,只能用 samtools 等软件转换为 Sam 格式后打开。

Sam 分为两部分,注释信息(header section)和比对结果(alignment section)。注释信息可有可无,每一行都是以@开头,用不同的 tag 表示不同的信息,tag 包括 @HD(符合标准的版本、对比序列的排列顺序说明)、@SQ(参考序列说明)、@RG(比对上的序列说明)、@PG(使用的程序说明)、@CO(任意的说明信息)。比对结果部分的每一行表示一个片段(segment)的比对信息,包括11个必须的字段(mandatory fields)和一个可选的字段,字段之间用 Tab 分割。示例及介绍如下:

1. QNAME,比对片段的编号

2. FLAG,位标识,1 表示该 read 是 pair 中的一条(read 表示本条 read,mate 表示 pair 中的另一条 read ),2 表示 pair 一正一负地比对上参考序列,4 表示这条 read 没有比对上,8 表示 mate 没有比对上,16 表示这条 read 比对上负链,32 表示 mate 比对上负链,64 表示这条 read 是 read1,128 表示这条 read 是 read2 等,FLAG 的值是符合情况的数字相加总和,即83=(64+16+2+1)表示该 read 为 read1,比对到负链上,其 mate 比对到正链上

3. RNAME,参考序列的编号

4. POS,比对上的位置,注意是从1开始计数,如果没有比对上,此处为0

5. MAPQ,比对的质量,越高则位点越独特,计算方法:Q = -10 log 10p,p 是该序列不来自这个位点的估计值

6. CIGAR(Compact Idiosyncratic Gapped Alignment Report),使用数字加字母表示比对结果,如 M 表示 match/mismatch,I 表示 insertion,D 表示 deletion等,数字表示碱基个数,即 42M4I5M 为该序列42个碱基匹配,4个 insertion,5个碱基匹配

7. RNEXT,mate 的名称,如果没有 mate,用 * 表示

8. PNEXT,mate 的位置,如果没有 mate,用 0 表示

9. TLEN, paired reads 间的距离,当 mate 序列位于本序列上游时该值为负值,如果比对区域仅有一个区段,或者不可用时,此处为0

10.SEQ,read 序列

11.QUAL,read 质量

12.Optional Fields,可选字段,格式如:TAG:TYPE:VALUE,其中 TAG 由两个大写字母组成,每个 TAG 代表一类信息,如 AS 表示匹配的得分,XS 表示第二好的匹配得分,YS表示 mate 序列匹配的得分等,TYPE 表示 TAG 对应值的类型,可以是字符串(Z)、整数(i)等

不同样品过滤后获得的数据量不可能完全一致,不同基因长度也有很大差异。为了能够在样品内(不同基因)以及样品间(不同分组)比较基因的表达量, 需要采用 FPKM 对表达量进行标准化(Normalization)。

FPKM(Fragments Per Kilobase Million),为每百万 Reads 中来自某一基因每千碱基长度的 Reads 数目,是一种普遍采用的基因表达量标准化方法,这种方法同时考虑了测序深度和基因长度对基因表达量计数的影响。其计算公式如下所示:

图1 FPKM 计算公式

目前以FPKM为基准的表达量标准化方法逐渐被其他统计学方法所取代,但是作为一种绝对化的标准化方法,其生物学意义明确,利于不同项目之间比较。在有参转录组当中, 我们一般认为 FPKM>1 的基因是表达的。这个阈值是主流杂志推荐的, 也能够很好的反应基因的表达水平。

此外,也可通过 RPKM 值描述基因表达量,FPKM 与 RPKM 计算方法基本一致,不同点在于: 对于 Pair-End 测序,每个 Fragments 会有两个 Reads,FPKM 只计算两个 Reads 能比对到同一个转录本的 Fragments 数量, RPKM 计算是比对到转录本的 Reads 数量。

gff 格式是一种 Sanger 研究所定义的,可以简单方便地描述 DNA、RNA 以及蛋白质序列的特征的数据格式,已经成为序列注释的通用格式,许多软件都支持输入或者输出 gff 格式。每一行代表一个特征条目(如基因、转录本、CDS、exon 等),每行有9列,以 Tab 为分割符,每列分别列出该特征条目的一些信息。 gff 可用文本编辑软件打开(如写字板、UltraEdit 、EditPlus 等工具)。目前 gff 的最新版本是版本3。示例如下:

ctg123PFAMgene10005000.+.ID=gene001;Name=EDEN
ctg123PFAMTF_binding_site10001012.+.Parent=gene001
ctg123PFAMmRNA10505000.+.ID=mRNA001;Parent=gene001
ctg123PFAMmRNA10505000.+.ID=mRNA002;Parent=gene001
ctg123PFAMexon13001500.+.Parent=mRNA001
ctg123PFAMexon10501500.+.Parent=mRNA001,mRNA002
ctg123PFAMCDS12013902.+0ID=cds001;Parent=mRNA001
ctg123PFAMCDS30004600.+2ID=cds001;Parent=mRNA001
ctg123PFAMCDS12011500.+1ID=cds002;Parent=mRNA002
  • 1. 序列编号,可能是染色体或者 scaffold 的名称
  • 2. 来源,产生这一特征条目的程序、数据库或者项目
  • 3. 类型,如 gene,transcript,CDS,mRNA,exon,five/three_prime_utr,start/stop_codon等
  • 4. 起始位点,这一特征条目在序列上的起始位置,从1开始计数
  • 5. 终止位点,这一特征条目在序列上的终止位置,不能大于序列的长度
  • 6. 得分,是注释信息可能性的说明,可以是序列相似性比对时的 E-values 值或者基因预测是的 P-values 值。“.”表示为空
  • 7. 序列的方向, + 表示正义链, - 反义链 , ? 表示未知
  • 8. 相位,仅对类型为 “CDS”的条目有效,有效值为0、1、2,0 表示这一特征条目的第一个碱基是一个密码子的第一个碱基,1 表示这一特征条目的第二个碱基是一个密码子的第一个碱基,以此类推
  • 9. 属性,以多个键值对组成的注释信息描述,键与值之间用“=”,不同的键值对用“;”隔开,一个键可以有多个值,不同值用“,”分割。键可以为ID(该特征条目的编号,在一个gff文件中必须唯一),Name(该特征条目的名称,可以重复),Parent(该特征条目的父级特征条目,值为父级特征条目的编号,比如外显子所属的转录本编号,转录本所属的基因的编号。值可以为多个)等

gtf 格式与 gff 格式前8列基本相同,不同之处在于第9列,虽然同样是标签与值配对的情况,但 gtf 格式的标签与值之间以空格分开,且每个属性之后都要有分号;(包括最后一个属性),而且第9列必须以 gene_id 以及 transcript_id 开头。

在生物信息学中,FASTA格式(又称为Pearson格式),是一种基于文本用于表示核苷酸序列或氨基酸序列的格式,可用文本编辑软件打开(如写字板、UltraEdit 、EditPlus等工具)。序列文件的第一行是由大于号“>”或分号“;”打头的任意文字说明(习惯常用“>”作为起始),用于序列标记。从第二行开始为序列本身,只允许使用既定的核苷酸或氨基酸编码符号。通常核苷酸符号大小写均可,而氨基酸常用大写字母。文件每行的字母一般不应超过80个字符。示例如下:

>Seq1

ADQLTEEQIAEFKEAFSLFDKDGDGTITTKELGTVMRSLGQNPTEAELQDMINEVDAD*