返回顶部
我们在分析有参考基因组的转录组数据时,需要将测序数据与参考基因组进行比对,本次项目使用的基因组信息如下表所示。
将数据库中获得的基因组信息整理如下表所示:
从不同的数据库中整理该物种的注释信息包括:1)基因的编号(ENSEMBL ID);2)染色体上位点信息(所在染色体序列编号、起始终止位点、长度、转录方向);3)基因的命名(HGVS Symbol、NCBI Gene ID、UniProtKB ID);4)基因的分类(GO、KO、EC分类);5)基因的文字注释信息。基因注释情况见下表。
首先对原始下机数据进行过滤以获得高质量序列,将过滤后的序列比对到该物种的参考基因组上。根据比对结果,计算每个基因的表达量。在此基础上,进一步对样本进行表达差异分析、富集分析和聚类分析;同时,对样本的转录本结构进行分析;获得样本基因的操纵子结构,UTR、cSNP和InDel等信息;对数据库中尚无注释的表达区域进行注释;最后将结果进行可视化展示。
每一个文库的基本情况见表1。
注:Sample:样品名称;
Lib. Name:文库名;
Lib. Insert Size:文库插入片段长度;
Sequencing platform:测序平台;
Sequencing Mode:测序模式。
样品经过上机测序,得到图像文件,由测序平台自带软件进行转化,生成 FASTQ 的原始数据(Raw Data),即下机数据。我们对每个样品的下机数据(Raw Data)分别进行统计,包括样品名、 Q30 、模糊碱基所占百分比、以及 Q20(%) 和 Q30(%) 。统计结果见表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
测序数据包含一些带接头 、低质量的 Reads,这些序列会对后续的信息分析造成很大的干扰,因此需要对测序数据进行进一步过滤。数据过滤的基本情况见表3。数据过滤的标准主要包括:
注:Sample:样品名;
Clean Reads No:高质量序列read数;
Clean Data (bp):高质量序列碱基数;
Clean Reads %: 高质量序列 reads 占测序 reads 的百分比;
Clean Data %: 高质量序列碱基占测序碱基的百分比。
测序错误率受测序仪本身、测序试剂、样品等多个因素共同影响。 对于 RNASeq 技术,测序错误率分布具有两个特点:
我们用测序数据的单碱基质量分布图评价单个位置的碱基质量。一般而言,Reads 的 5’端和 3’端的碱基质量较低,中间部分的碱基质量较高。大部分序列的碱基质量在20以上,代表测序质量较好。
注:横坐标是 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
碱基含量分布一般用于检测有无AT、GC分离现象。对于RNASeq来说,鉴于序列打断的随机性和G/C、A/T含量分别相等的原则,理论上每个测序循环中的GC含量相等、AT含量相等(如果是链特异性建库,可能会出现AT分离和/或GC分离),且在整个测序过程基本稳定不变,呈水平线。但在现有的高通量测序技术中,反转录合成 cDNA 时所用的6bp的随机引物会引起前几个位置的核苷酸组成存在一定的偏好性,这种波动属于正常情况。碱基含量分布结果见图2。
注:横坐标是 Reads 中碱基位置(5’->3’),纵坐标是该位点某碱基所占的比例统计。
samples Read1 碱基含量分布图:./result/02_FastQC/samples_R1_fastqc/per_base_sequence_content.png
samples Read2 碱基含量分布图:./result/02_FastQC/samples_R2_fastqc/per_base_sequence_content.png
Reads 平均质量分布主要用来检测测序数据的平均质量分布情况。峰尖代表主体 Reads 测序质量,峰宽表示整体测序质量分布,较大峰宽或者前拖尾峰表示测序数据的部分数据质量偏低,峰尖对应值较低表明整体测序结果较差,如果没有峰尖,则表示 Reads 整体质量很好。Reads 平均质量分布结果见图3。
注:Reads质量分布图,横坐标表示 Reads 的平均质量,纵坐标为对应平均质量值的 Reads 数目。
samples Read1 平均质量分布图:./result/02_FastQC/samples_R1_fastqc/per_sequence_quality.png
samples Read2 平均质量分布图:./result/02_FastQC/samples_R2_fastqc/per_sequence_quality.png
通过Bowtie2建立参考基因组索引,然后使用Bowtie2(http://bowtie-bio.sourceforge.net/index.shtml)将过滤后的Reads比对到参考基因组上。序列比对的基本信息统计见下表。
注:Sample:样品;
Useful Reads:用于比对的序列总数;
Total Mapped Reads(%):比对上参考基因组的序列总数及比例;
Uniquely Mapped Reads(%):只比对到一个位置的序列总数及比例;
Multiple Mapped Reads(%):比对到多个位置的序列总数及比例。
samples 比对结果:./result/03_Map/samples.accepted_hits.bam
比对到基因上不同结构区域的Reads数目,除了与样品有关,还与该物种基因组注释的情况有关。所使用的参考基因组注释的越全面,比对统计结果就越能反映样品实际情况。
注:Sample:样品名;
Total Mapped Reads:比对上参考基因组的序列总数;
InterGene(%):比对到基因间区域的Reads总数及比例;
Gene(%):比对到基因区域的Reads总数及比例;
mRNA(%):比对到蛋白编码基因上的Reads总数及比例;
rRNA(%):比对到rRNA基因间上的Reads总数及比例。
samples Reads 在基因区域分布图(png格式):./images/Read_Distribution/samples.read_distribution.png
samples Reads 在基因区域分布图:./result/04_MapQC/read_distribution/samples.read_distribution.pdf
测序 Reads 在基因上覆盖度的分布情况,见图6。展示每个样品所有基因的5'到3'区域上序列覆盖情况,用于评估测序结果的均一性(或是否有偏向性)。理想条件下,Reads在所有表达的基因上的分布应该呈现均一化分布。
注:横坐标为单个基因的碱基长度占总碱基长度的百分比,0表示基因的5'端,100表示基因的3'端;纵坐标为比对到所有基因的横轴位置上相应区间内的序列条数的总和。图中体现了所有基因覆盖情况的叠加结果,曲线中每个点的纵坐标表示所有基因在该相对比例位置上所有序列的数量;曲线反映了测序所得序列是否在基因上均匀分布。若无明显偏向锋,则说明测序无偏向性。
samples 基因上测序深度分布图(png格式):./images/Gene_Body_Coverage/samples.geneBodyCoverage.curves.png
samples 基因上测序深度分布图:./result/04_MapQC/genebody_coverage/samples.geneBodyCoverage.curves.pdf
我们使用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 方向一致)。
注:Gene_ID:基因ID
*read count:样本*的基因表达量(Reads数);
*.fpkm:样本*的基因表达量(FPKM值)。
根据表达量的计算结果表格,将表达量分为不同的区间,对各样本在不同表达量区间内的基因的数目进行统计,结果见下图。
注:横坐标表示表达量FPKM值范围,纵坐标表示该表达量区间基因的个数。
根据表达量的计算结果表格,对各样本中鉴定到的mRNA进行统计,并对各样本间共有特有的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
FPKM密度分布能整体考察样品的基因表达模式,一般来说中等表达的基因占绝大多数,低表达和高表达的基因占一小部分。
注:横坐标为基因的log10(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
我们使用RSeQC分析测序饱和度,即:对测序结果进行抽样,对不同的抽样比例分别获得所有基因的RPKM值(标准化的基因表达量),然后与全部测序数据得到的RPKM值行比较,获得相对误差。如果随抽样比例的增加,相对误差减少,说明测序结果趋于饱和。
注:该图为RPKM饱和图,横坐标是重采样的比率,纵坐标为相对误差。
samples RPKM 饱和度图(png格式):./images/Saturation/samples.saturation.png
samples RPKM 饱和度图:./result/04_MapQC/RPKM_saturation/samples.saturation.pdf
样品间基因表达水平相关性是检验实验可靠性和样本选择是否合理的重要指标,在做差异表达分析之前,应先检查样品间基因的表达水平相关性。我们用皮尔逊相关系数表示样品间基因的表达水平相关性,相关系数越接近1,表明样品之间表达模式的相似度越高。一般来说,相关系数在0.8-1之间属于极强相关,如果生物学重复的样本之间相关系数低于0.8,表示样品之间的相关性较低,见图7。
注:左侧和上侧为样本聚类情况,图中右侧和下侧为样本名称,不同颜色的方块代表两个样本的相关性高低情况。
样品相关性检验图(png格式):./images/mRNA/1_Expression/cor_pearson.png
PCA 主成分分析(Principal Components Analysis),通过线性变换,将高维数据降低至二维或三维,同时保持各方差贡献最大的特征,即降低数据复杂度。当有多个样品时,我们使用R语言的 DESeq软件包,根据表达量对各样品进行 PCA 主成分分析。PCA 分析可以把相似的样本聚到一起,距离越近表明样本间相似性越高。结果见图8, PCA图中Replication与sample的对应关系见分析结果中的PCA.xlsx。
注:横坐标为第一主成分,纵坐标为第二主成分。在图中的不同形状表示不同的样品,不同的颜色表示不同的分组。
PCA分析图(png格式):./images/mRNA/1_Expression/PCA.png
我们采用 DESeq 对基因表达进行差异分析,筛选差异表达基因条件为:表达差异倍数 |log2FoldChange| > 1 ,显著性P-value<0.05 。部分差异分析结果见表2,不同分组之间的差异表达基因统计结果见表3。
注:id:基因编号;
baseMean:表示进行差异比较的两组全部样品该基因read count的均一化结果;
baseMean(sample):表示该组全部样品该基因read count的均一化结果;
foldChange(Case/Control):表达差异倍数;
log2FoldChange:表达差异倍数对2取对数的值;
pval:显著性 p-value;
padj:校正后的显著性 p-value。
注:Case:实验组样本;
Control:对照组样本;
Up-regulated Genes:Case 相比于Control 上调表达基因;
Down-regulated Genes:Case 相比于Control 下调表达的基因;
Total DEGs:Case 相比于Control 差异表达基因。
注:横坐标表示进行差异分析的比较组,纵坐标表示差异基因的个数,颜色中红色表示上调基因,绿色表示下调基因。
采用R语言 ggplots2 软件包绘制差异表达基因的火山图和 MA 图,见图10。火山图展示的是基因分布情况,基因的表达倍数差异和显著性结果,正常状况下,该图左右差异基因分布应大致对称,左侧为 Case 相比于 Control 下调基因,右侧为 Case 相比于 Control 上调基因。MA 图常用来展示标准化后的基因分布情况,一般标准化后,基因的表达量呈对称分布,即表达差异趋势不随基因表达量变化而发生偏向。
注:火山图:横坐标为 log2FoldChange,纵坐标为 -log10(p-value)。图中两条竖虚线为2倍表达差异阈值;横虚线为 P-value=0.05 阈值。红点表示该组上调基因,蓝点表示该组下调基因,灰点表示非显著差异表达基因。
注:横坐标为两样品基因表达量之积对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
我们使用R语言Circlize包,根据基因组信息和RNA差异表达分析结果,在基因组上标记差异表达的RNA。
注:最外圈是染色体条带,从外向内为不同差异分析的差异表达分析结果。红色和绿色分别为上调和下调基因的 log2FoldChange 值的柱状图,灰色为无差异表达基因的 log2FoldChange 值的散点图。
根据差异分析的结果,统计各比较组之间的共有特有差异基因数量,如图13。维恩图(仅提供两组、 三组、 四组和五组比较的维恩图)展示了各比较组间差异基因的个数,以及比较组间的重叠关系,Upset图(仅提供两组以上比较的矩阵图)展示了两两比较组间共有的差异基因数目。
注: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
聚类分析用于判断差异表达基因在不同实验条件下的表达模式;在样品间表达量相关性高的基因被归为一类,通常这些基因在某些生物学过程,或者某个代谢、信号通路中存在实际的联系。因此通过表达量聚类我们可以发现基因间未知的生物学联系。
我们使用R语言 Pheatmap 软件包对所有比较组的差异基因的并集和样品进行双向聚类分析,根据同一基因在不同样品中的表达水平和同一样品中不同基因的表达模式进行聚类,采用 Euclidean 方法计算距离,层次聚类最长距离法(Complete Linkage)进行聚类。结果见图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
根据上一节聚类分析中层次聚类的结果,将差异基因按照表达模式的不同,划分到不同的Cluster(同一个Cluster中的基因表达趋势相近),从而获得基因的聚类结果。结果展示见图15。
注:图中灰色线展示每个Cluster中基因的表达模式,蓝色的线表示Cluster中的所有基因在样品中表达量的平均值。
趋势分析图(png格式):./images/mRNA/2_DEG/ALL_cluster_plots.png
STRING (Search3 Tool for the Retrieval of Interacting Genes/Proteins)是EMBL开发的蛋白质互作数据库,该数据库从最有力的实验证据到数据挖掘、同源预测的蛋白质互作关系都有收录。 当STRING数据库中收录该物种的PPI信息时,我们根据基因差异表达分析结果,直接数据库中筛选双端节点皆为差异基因并且Score>0.95的PPI作用对。当STRING数据库中没有该物种的PPI信息时,我们选择相近物种与该物种的序列进行比对,进而得到该物种的蛋白质之间的相互关系。 分析结果见图17。
注:一条连线代表一组相互关系;节点大小表示该节点的度。
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
使用topGO进行GO富集分析,分析时利用GO term 注释的差异基因对每个term 的基因列表和基因数目进行计算,然后通过超几何分布方法计算P-value(显著富集的标准为P-value < 0.05),找出与整个基因组背景相比,差异基因显著富集的GO term ,从而确定差异基因行使的主要生物学功能。
注: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。
注:横坐标为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
富集分析结果会分别给出GO三个 ontology(细胞组分、分子功能和生物过程)的有向无环图,在这个图中,越接近根结点的 GO term越概括,往下分支的 GO term为注释到更细层级的 term。 程序默认把显著性最高前10个GO term设置为方形,其他的GO term为圆形。颜色越深,代表该GO term越显著,颜色由浅到深为无色-浅黄-深黄-红色,见图20。
注:每个节点代表一个 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
注:PathwayID:通路ID;
Pathway:通路名称;
Up/Down_number:富集到该条通路下的上/下调基因数量;
DEG_number:富集到该条通路下的差异基因总数;
total_number:该条通路下的基因总数;
Up/Down_gene:富集到该条通路下的上/下调基因ID(名称)。
根据差异表达的基因的 KEGG 富集分析结果,挑选 p-value 最小即富集最显著的前20个 Pathway 进行展示,结果见图21。
注:横坐标为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
ko 编号表示一个通路,这个通路是不分物种的,相当于所有物种的这一通路的并集。K 编号表示一个基因,是 ko 通路中的基本单位,某一 K 编号代表的不是某一具体物种的基因,而是所有物种的某一同源基因的统称。在 KEGG_enrichment.xls(具体位置见下方结果文件)中 URL 列为该行 Pathway 在 KEGG 数据库中对应的代谢通路图的链接,打开 URL 列的网页,可查看实验组和对照组差异表达基因在 Pathway 通路上的不同。下图为一个代谢通路图示例(仅为图片示例,不能演示 KEGG 网站中的动态效果),标红色的节点包含实验组上调基因;标绿色的节点包含对照组上调基因;标紫色的节点既包含实验组上调基因,也包含对照组上调基因。将鼠标指在节点上将会弹出节点的详细描述,点击各个节点可以跳转到 KEGG 数据库中对应的具体信息页。右键点击图片可选择保存该代谢通路图。
注:方框一般表示酶,小圆圈表示代谢物,圆角框表示另一个代谢通路图。
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
cSNP 包括转换和颠换两种。SNP 在 CG 序列上出现的最为频繁,而且多是 C 转换为 T,原因是 CG 中的 C 常为甲基化的,自发地脱氨后即为胸腺嘧啶。SNP 出现的原因有很多,有的是遗传背景的单核苷酸多态性,有的是建库技术上造成的突变,有的则可能是测序中读取错误。InDel可能造成移码突变,导致 mRNA 翻译时遇上一个错误的终止密码子。一般 InDel 不是3的倍数的情况在编码区不常发生,在非编码区相对频繁地发生。除了在高度重复区域附近,InDel 发生的频率一般会低于 SNP。
我们采用 Varscan 程序获取 SNP 和 InDel 位点,过滤标准为:
SNP 分析结果见表6,InDel 分析结果格式同 SNP 分析结果,SNP、InDel 和转换/颠换统计见图24。
注:CHROM:SNP位点所在染色体;
POS:SNP位点在染色体上的位置;
REF/ALT:参考序列在该位点的基因型/突变基因型;
样品名下的列:支持各基因型的 Reads 数目。
注:Homo(homozygous-variant) 表示纯合子变异体,即这一位点的等位基因都突变了,并且突变相同,Hete(heterozygous-variant)表示杂合子变异体,即这一位点的等位基因至少有一个突变了,且突变后等位基因不同。
SNP 统计图(png格式):./images/mRNA/4_Structure/Snp_Stat.png
Indel 统计图(png格式):./images/mRNA/4_Structure/Indel_Stat.png
转换/颠换统计图(png格式):./images/mRNA/4_Structure/Ts_Tv_Stat.png
InDel 统计图:./result/mRNA/4_Structure/SNP/Indel_Stat.pdf
SNP 统计图:./result/mRNA/4_Structure/SNP/Snp_Stat.pdf
转换/颠换统计图:./result/mRNA/4_Structure/SNP/Ts_Tv_Stat.pdf
Indel 分析结果:./result/mRNA/4_Structure/SNP/indel.xlsx
完整的操纵子(Operon)包括操纵子基因以及与之相关的调节基因(Regulator),如下图是典型的乳糖操纵子。我们通过Rockhopper预测操纵子转录起始位点(Transcription Start Site, TSS)、转录终止位点(Transcription Termination Site, TTS)。
注:Start:第一个基因起始位置;
Stop.:最后一个基因的终止位置;
Strand.:方向;
Number of Genes:基因数量;
Genes:基因列表。
注:TSS:转录起始位点;
TSS:转录终止位点;
Strand.:方向;
Gene name:基因名;
Product:基因产物。
我们根据基因转录起始位点(转录终止位点)和翻译起始位点(翻译终止位点)信息,提取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/)软件对不依赖σ因子的终止子进行预测。
注:Chromosome:染色体编号,;
Start:UTR起始位置;
End:UTR终止位置;
Name:UTR类型|基因ID|基因名;
Length:UTR长度。
注:Seq.no:序列顺序;
Pos:在染色体上的位置;
*****,Motif,*****:RBS序列(及邻近的序列);
Prob:概率;
Seq.Id:序列名。
注: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
基因本体论联合会建立的数据库(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等十多个二次数据库的交叉引用代码。
国际生物化学会酶学委员会(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 是那些在一定物种中的来源于基因复制的蛋白,可能会进化出新的与原来有关的功能。
[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.
FASTQ 格式(http://en.wikipedia.org/wiki/FASTQ_format)是一种文本格式,常用于存储生物学序列及其对应的质量分值。FASTQ 格式文件可以采用文本编辑软件(如写字板、UltraEdit 、EditPlus等工具)打开,由于文件较大,对电脑的内存要求较高。FASTQ 格式中,每个 Read 由四行信息表示。
第一行为序列名称,以 @ 开头,其后是序列描述;第二行为碱基序列;第三行为 “+” 号,不代表任何意义;第四行碱基质量,与第二行的碱基序列一一对应。示例如下:
我们使用 Sanger 质量值来评估下机数据的测序质量。质量值,简称 Q 值,是碱基读取错误率 p 的取整映射结果,等于 Phred quality score,计算公式为:
Qphred = -10 log10 p
Phred quality score计算公式测序错误率与 Q 值的简明对应方式如表2所示。
测序错误率 | 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%。
十进制 | 字符 | Q值 | 十进制 | 字符 | Q值 | 十进制 | 字符 | Q值 | 十进制 | 字符 | Q值 |
---|---|---|---|---|---|---|---|---|---|---|---|
32 | 48 | 0 | 15 | 64 | @ | 31 | 80 | P | 47 | ||
33 | ! | 0 | 49 | 1 | 16 | 65 | A | 32 | 81 | Q | 48 |
34 | “ | 1 | 50 | 2 | 17 | 66 | B | 33 | 82 | R | 49 |
35 | # | 2 | 51 | 3 | 18 | 67 | C | 34 | 83 | S | 50 |
36 | $ | 3 | 52 | 4 | 19 | 68 | D | 35 | 84 | T | 51 |
37 | % | 4 | 53 | 5 | 20 | 69 | E | 36 | 85 | U | 52 |
38 | & | 5 | 54 | 6 | 21 | 70 | F | 37 | 86 | V | 53 |
39 | ’ | 6 | 55 | 7 | 22 | 71 | G | 38 | 87 | W | 54 |
40 | ( | 7 | 56 | 8 | 23 | 72 | H | 39 | 88 | X | 55 |
41 | ) | 8 | 57 | 9 | 24 | 73 | I | 40 | 89 | Y | 56 |
42 | * | 9 | 58 | : | 25 | 74 | J | 41 | 90 | Z | 57 |
43 | + | 10 | 59 | : | 26 | 75 | K | 42 | 91 | [ | 58 |
44 | , | 11 | 60 | < | 27 | 76 | L | 43 | 92 | 59 | |
45 | - | 12 | 61 | = | 28 | 77 | M | 44 | 93 | ] | 60 |
46 | . | 13 | 62 | > | 29 | 78 | N | 45 | 94 | ^ | 61 |
47 | / | 14 | 63 | ? | 30 | 79 | O | 46 | 95 | _ | 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 数目,是一种普遍采用的基因表达量标准化方法,这种方法同时考虑了测序深度和基因长度对基因表达量计数的影响。其计算公式如下所示:
目前以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。示例如下:
ctg123 | PFAM | gene | 1000 | 5000 | . | + | . | ID=gene001;Name=EDEN |
ctg123 | PFAM | TF_binding_site | 1000 | 1012 | . | + | . | Parent=gene001 |
ctg123 | PFAM | mRNA | 1050 | 5000 | . | + | . | ID=mRNA001;Parent=gene001 |
ctg123 | PFAM | mRNA | 1050 | 5000 | . | + | . | ID=mRNA002;Parent=gene001 |
ctg123 | PFAM | exon | 1300 | 1500 | . | + | . | Parent=mRNA001 |
ctg123 | PFAM | exon | 1050 | 1500 | . | + | . | Parent=mRNA001,mRNA002 |
ctg123 | PFAM | CDS | 1201 | 3902 | . | + | 0 | ID=cds001;Parent=mRNA001 |
ctg123 | PFAM | CDS | 3000 | 4600 | . | + | 2 | ID=cds001;Parent=mRNA001 |
ctg123 | PFAM | CDS | 1201 | 1500 | . | + | 1 | ID=cds002;Parent=mRNA002 |
gtf 格式与 gff 格式前8列基本相同,不同之处在于第9列,虽然同样是标签与值配对的情况,但 gtf 格式的标签与值之间以空格分开,且每个属性之后都要有分号;(包括最后一个属性),而且第9列必须以 gene_id 以及 transcript_id 开头。
在生物信息学中,FASTA格式(又称为Pearson格式),是一种基于文本用于表示核苷酸序列或氨基酸序列的格式,可用文本编辑软件打开(如写字板、UltraEdit 、EditPlus等工具)。序列文件的第一行是由大于号“>”或分号“;”打头的任意文字说明(习惯常用“>”作为起始),用于序列标记。从第二行开始为序列本身,只允许使用既定的核苷酸或氨基酸编码符号。通常核苷酸符号大小写均可,而氨基酸常用大写字母。文件每行的字母一般不应超过80个字符。示例如下:
>Seq1
ADQLTEEQIAEFKEAFSLFDKDGDGTITTKELGTVMRSLGQNPTEAELQDMINEVDAD*