有参考基因组的转录组生物信息分析模板

一、生物信息分析流程

获得原始测序序列(Sequenced Reads)后，在有相关物种参考序列或参考基因组的情况下，通过如下流程进行生物信息分析：

二、项目结果说明

1 原始序列数据

高通量测序(如illumina HiSeqTM2000/MiSeq等测序平台)测序得到的原始图像数据文件经碱基识别(Base Calling)分析转化为原始测序序列(Sequenced Reads)，我们称之为Raw Data或Raw Reads，结果以FASTQ(简称为fq)文件格式存储，其中包含测序序列(reads)的序列信息以及其对应的测序质量信息。

FASTQ格式文件中每个read由四行描述，如下：

@EAS139:136:FC706VJ:2:2104:15343:197393 1:Y:18:ATCACG GCTCTTTGCCCTTCTCGTCGAAAATTGTCTCCTCATTCGAAACTTCTCTGT +

@@CFFFDEHHHHFIJJJ@FHGIIIEHIIJBHHHIJJEGIIJJIGHIGHCCF

其中第一行以“@”开头，随后为illumina 测序标识符(Sequence

Identifiers)和描述文字(选择性部分)；第二行是碱基序列；第三行以“+”开头，随后为illumina 测序标识符(选择性部分)；第四行是对应序列的测序质量(Cock et al.)。

illumina 测序标识符详细信息如下：

EAS139 136 FC706VJ 2 2104 15343 197393 1 Y 18 ATCACG Unique instrument name Run ID Flowcell ID Flowcell lane Tile number within the flowcell lane 'x'-coordinate of the cluster within the tile 'y'-coordinate of the cluster within the tile Member of a pair, 1 or 2 (paired-end or mate-pair reads only) Y if the read fails filter (read is bad), N otherwise 0 when none of the control bits are on, otherwise it is an even number Index sequence 第四行中每个字符对应的ASCII值减去33，即为对应第二行碱基的测序质量值。如果测序错误率用e表示，illumina HiSeqTM2000/MiSeq的碱基质量值用Qphred表示，则有下列关系：

公式一： Qphred = -10log10(e)

illumina Casava 1.8版本测序错误率与测序质量值简明对应关系如下：

测序错误率 5% 1% 0.1% 0.01% 测序质量值 13 20 30 40 对应字符 . 5 ? I

2 测序数据质量评估

2.1 测序错误率分布检查

每个碱基测序错误率是通过测序Phred数值(Phred score, Qphred)通过公式1转化得到，而Phred 数值是在碱基识别(Base Calling)过程中通过一种预测碱基判别发生错误概率模型计算得到的，对应关系如下表所显示：

illumina Casava 1.8版本碱基识别与Phred分值之间的简明对应关系

Phred分值 10 20 30 40 不正确的碱基识别 1/10 1/100 1/1000 1/10000 碱基正确识别率 90% 99% 99.9% 99.99% Q-sorce Q10 Q20 Q30 Q40 测序错误率与碱基质量有关，受测序仪本身、测序试剂、样品等多个因素共同影响。对于RNA-seq技术，测序错误率分布具有两个特点：

(1)测序错误率会随着测序序列(Sequenced Reads)的长度的增加而升高，这是由于测序过程中化学试剂的消耗而导致的，并且为illumina高通量测序平台都具有的特征(Erlich and Mitra, 2008; Jiang et al.)。

(2)前6个碱基的位置也会发生较高的测序错误率，而这个长度也正好等于在RNA-seq建库过程中反转录所需要的随机引物的长度。所以推测前6个碱基测序错误率较高的原因为随机引物和RNA模版的不完全结合(Jiang et al.)。测序错误率分布检查用于检测在测序长度范围内，有无异常的碱基位置存在高错误率，比如中间位置的碱基测序错误率显著高于其他位置。一般情况下，每个碱基位置的测序错误率都应该低于0.5%。

图2.1 测序错误率分布图

横坐标为reads的碱基位置，纵坐标为单碱基错误率

2.2 GC含量分布检查

GC含量分布检查用于检测有无AT、GC 分离现象，而这种现象可能是测序或者建库所带来的，并且会影响后续的定量分析。

在illumina测序平台的转录组测序中，反转录成cDNA时所用的6bp 的随机引物会引起前几个位置的核苷酸组成存在一定的偏好性。而这种偏好性与测序的物种和实验室环境无关，但会影响转录组测序的均一化程度(Hansen et al.)。除此之外，理论上G和C碱基及A和T碱基含量每个测序循环上应分别相等，且整个测序过程稳定不变，呈水平线。对于DGE测序来说，由于随机引物扩增偏差等原因，常常会导致在测序得到的每个read前6-7个碱基有较大的波动，这种波动属于正常情况。

图2.2 GC含量分布图

图3.2 Reads在参考基因组不同区域的分布情况

3.3 Reads在染色体上的密度分布情况

对Total mapped reads的比对到基因组上的各个染色体（分正负链）的密度进行统计，如下图所示，具体作图的方法为用滑动窗口(window size)为1K，计算窗口内部比对到碱基位置上的reads的中位数，并转化成 log2 。正常情况下，整个染色体长度越长，该染色体内部定位的reads总数会越多(Marquez et al.)。从定位到染色体上的reads数与染色体长度的关系图中，可以更加直观看

出染色体长度和reads总数的关系。

图3.3 Reads在染色体上的密度分布图

上图：横坐标为染色体的长度信息(以百万碱基为单位)，纵坐标为log2(reads的密度的中位数)，绿色为正链，红色为负链 下图：横坐标为染色体的长度信息(单位为Mb)，纵

坐标为mapped到染色体上的reads数(单位为M)

3.4 Reads比对结果可视化

我们提供RNA-seq Reads在基因组上比对结果的bam格式文件，部分物种还提供相应的参考基因组和注释文件，并推荐使用IGV (Integrative Genomics Viewer) 浏览器对bam文件进行可视化浏览。IGV浏览器具有以下特点：(1)能在不同尺度下显示单个或多个读段在基因组上的位置，包括读段在各个染色体上的分布情况和在注释的外显子、内含子、剪接接合区、基因间区的分布情况等；(2)能在不同尺度下显示不同区域的读段丰度，以反映不同区域的转录水平；(3)能显示基因及其剪接异构体的注释信息；(4)能显示其他注释信息；(5)既可以从远程服务器端下载各种注释信息，又可以从本地加载注释信息。IGV浏览器使用方法可参考我们提供的使用说明文档(IGVQuickStart.pdf)。

图3.4 IGV浏览器界面

4 可变剪切分析

用ASprofile软件对Cufflinks (Trapnell et al.)预测出的基因模对每个样品的可变剪切事件分别进行分类和表达量统计。分析流程及ASprofile中的可变剪切事件分类如下图所示：

12类可变剪切事件定义如下:

(1) TSS: Alternative 5' first exon (transcription start site) 第一个外显子可变剪切

(2) TTS: Alternative 3' last exon (transcription terminal site) 最后一个外显子可变剪切

(3) SKIP: Skipped exon (SKIP_ON,SKIP_OFF pair) 单外显子跳跃 (4) XSKIP: Approximate SKIP (XSKIP_ON,XSKIP_OFF pair) 单外显子跳跃（模糊边界）

(5) MSKIP: Multi-exon SKIP (MSKIP_ON,MSKIP_OFF pair) 多外显子跳跃 (6) XMSKIP: Approximate MSKIP (XMSKIP_ON,XMSKIP_OFF pair) 多外显子跳跃（模糊边界）

(7) IR: Intron retention (IR_ON, IR_OFF pair) 单内含子滞留

(8) XIR: Approximate IR (XIR_ON, XIR_OFF pair) 单内含子滞留（模糊边界）

(9) MIR: Multi-IR (MIR_ON, MIR_OFF pair) 多内含子滞留

(10) XMIR: Approximate MIR (XMIR_ON, XMIR_OFF pair) 多内含子滞留（模糊边界）

(11) AE: Alternative exon ends (5', 3', or both) 可变 5'或3'端剪切

(12) XAE: Approximate AE 可变 5'或3'端剪切（模糊边界） 4.1 可变剪切事件分类和数量统计

图4.1 AS分类和数量统计

纵轴为可变剪切事件的分类缩写，横轴为该种事件下可变剪切的数量，不同样品用不同子

图和颜色区分

4.2 可变剪切事件结构和表达量统计

表4.2 AS结构和表达量统计

event_ievent_tyd 1000001 100000pe TSS TSS chroevent_staevent_enevent_pattestranm 1 1 rt d rn d + + gene_id CUFF.100 CUFF.10fpkm ref_id 3438277 3438330 3438330 3450218 3450253 3450253 1.00000000ENSGALT0000001000 225 3.00000000ENSGALT000000102 1000003 1000004 TSS TTS 0 CUFF.100 CUFF.100 1 1 3456744 3457165 3457165 3474806 3478178 3474806 + + 00 00 00 225 225 225 2.00000000ENSGALT000000105.00000000ENSGALT00000010(1) event_id: AS事件编号 (2) event_type: AS事件类型 (TSS, TTS, SKIP_{ON,OFF}, XSKIP_{ON,OFF}, MSKIP_{ON,OFF}, XMSKIP_{ON,OFF}, IR_{ON ,OFF}, XIR_{ON,OFF}, AE, XAE) (3) gene_id: cufflink组装结果中的基因编号 (4) chrom: 染色体编号

(5) event_start: AS事件起始位置 (6) event_end: AS事件结束位置

(7) event_signature: AS事件特征 (for TSS, TTS - inside boundary of alternative marginal exon; for *SKIP_ON,the coordinates of the skipped exon(s); for *SKIP_OFF, the coordinates of the enclosing introns; for *IR_ON, the end coordinates of the long,

intron-containing exon; for *IR_OFF, the listing of coordinates of all the exons along the path containing the retained intron; for *AE, the coordinates of the exon variant) (8) strand: 基因正负链信息

(9) fpkm: 此AS类型该基因表达量

(10) ref_id: 此基因在参考注释文件中的编号

5 新转录本预测

将所有测序reads数据的基因组定位结果放到一起，用 Cufflinks 进行组装，然后用Cuffcompare和已知的基因模型进行比较，可以:(1)发现新的未知基因（相对于原有基因注释文件）；(2)发现已知基因新的外显子区域；(3)对已知基因的起始和终止位置进行优化。新基因和新外显子区域预测结果为GTF格式的注释文件。GTF格式的详细说明可参考(http://mblab.wustl.edu/GTF22.html)

表5.1 新转录本结构注释结果

seqname source feature start end score strand frame 1 1 1 1 novelGene exon 18531 19499 . novelGene exon 20813 21813 . novelGene exon 23917 24402 . novelGene exon 25189 26100 . + + + + . . . . attributes gene_id \\gene_id \\gene_id \\gene_id \\(1) seqname：染色体编号

(2) source：来源标签，这里的novelGene指新基因 (3) feature：区域类型，目前我们预测外显子区域 (4) start：起始坐标 (5) end：终止坐标 (6) score：不必关注 (7) strand：正负链信息 (8) frame：不必关注

(9) attributes：属性，包括基因编号、转录本编号等信息

表5.2 已知基因结构优化

Gene_id ENSGALG00000000003 ENSGALG00000000004 ENSGALG00000000011 ENSGALG00000000013 Chromosome Strand 1 Z 6 22 + - - + Original_span 19947199~19967662 34293201~34299523 30928671~30932616 2783575~2787337 Assembled_span 19947199~19967662 34293201~34299554 30928671~30932616 2783575~2787453 (1) Gene_id：原注释文件中基因命名编号 (2) Chromosome：染色体编号 (3) Strand：正负链信息

(4) Original_span：原注释文件中基因起始位置~终止位置

(5) Assembled_span：转录组拼接结果中基因起始位置~终止位置

6 SNP和Indel分析

SNP全称Single Nucleotide Polymorphisms，是指在基因组上由单个核苷酸变异形成的遗传标记，其数量很多，多态性丰富。从理论上来看每一个SNP 位点都可以有4 种不同的变异形式，但实际上发生的只有两种，即转换和颠换，二者之比为1:2。SNP在CG序列上出现最为频繁，而且多是C转换为T，原因是CG中的C常为甲基化的，自发地脱氨后即成为胸腺嘧啶。一般而言，SNP是指变异频率大于1%的单核苷酸变异。Indel(insertion-deletion)是指相对于参考基因组，样本中发生的小片段的插入缺失，该插入缺失可能含一个或多个碱基。 我们通过samtools和picard-tools等工具对比对结果进行染色体坐标排序、去掉重复的reads等处理，最后通过变异检测软件GATK(McKenna et al., 2010)分别进行SNP Calling和Indel Calling，并对原始结果进行过滤，得到如下表形式的分析结果。其中Indel分析结果每列的含义和SNP结果是一致的。

表6 SNP分析结果

#CHROM 1 POS 502 REF A ALT G HS1 .. HS2 .. HT1 11 HT2 .. HW1 00 HW2 .. 1 1 1 563 1213 1316 C A G A G A .. .. 00 .. .. 00 11 11 01 11 .. .. 00 .. 00 .. .. 00 #CHROM：SNP位点所在染色体 POS：SNP位点坐标

REF：参考序列在该位点的基因型 ALT：该位点的其它基因型 other coloums：每个个体该位点的基因型(0 与REF一致；1 与ALT一致；. 缺少数据支持)

7 基因表达水平分析

一个基因表达水平的直接体现就是其转录本的丰度情况，转录本丰度程度越高，则基因表达水平越高。在RNA-seq分析中，我们可以通过定位到基因组区域或基因外显子区的测序序列(reads)的计数来估计基因的表达水平。Reads计数除了与基因的真实表达水平成正比外，还与基因的长度和测序深度成正相关。为了使不同基因、不同实验间估计的基因表达水平具有可比性，人们引入了RPKM的概念，RPKM(Reads Per Kilo bases per Million reads)是每百万reads中来自某一基因每千碱基长度的reads数目。RPKM同时考虑了测序深度和基因长度对reads计数的影响，是目前最为常用的基因表达水平估算方法 (Mortazavi et al., 2008)。

结果文件分别统计了不同表达水平下基因的数量以及单个基因的表达水平。一般情况下，RPKM数值0.1或者1作为判断基因是否表达的阈值，不同的文献所采用的阈值不同。

表7.1 不同表达水平区间的基因数量统计表

RPKM Interval 0-1 1-3 3-15 15-60 >60 HS1 HS2 HT1 HT2 HW1 HW2 11678(44.62%) 11157(42.63%) 11644(44.49%) 11552(44.14%) 11663(44.56%) 11652(44.52%) 3416(13.05%) 3829(14.63%) 3497(13.36%) 3622(13.84%) 3359(12.83%) 3503(13.39%) 6586(25.17%) 6741(25.76%) 6719(25.67%) 6731(25.72%) 6441(24.61%) 6522(24.92%) 3436(13.13%) 3421(13.07%) 3278(12.53%) 3277(12.52%) 3612(13.80%) 3442(13.15%) 1055(4.03%) 1023(3.91%) 1033(3.95%) 989(3.78%) 1096(4.19%) 1052(4.02%) 表7.2 基因表达水平统计表

geneID 000003 HS1 01677 HS2 40486 HT1 85256 HT2 70462 HW1 19692 HW2 09366 ENSGALG000000.1774627480.57509924280.18207677010.32875021680.22426855480.2831623894ENSGALG000009.4448453447.21882264978.93638040289.694005382810.4063844079.5631290234000004 000011 000013 41505 60835 18205 4696 8464 7827 376 6865 0897 1103 6002 7573 8955 413 2031 8064 986 3859 ENSGALG0000014.7023552914.23548457111.32025669512.29436572213.35758690510.060029486ENSGALG000006.7488688417.48873293736.66826144957.41338358787.53637022217.50828610128 RNA-seq整体质量评估

8.1 表达水平的饱和曲线检查

定量饱和曲线检查反映了基因表达水平定量对数据量的要求。表达量越高的基因，就越容易被准确定量；反之，表达量低的基因，需要较大的测序数据量才能被准确定量。

表达水平的饱和曲线的具体算法描述如下：分别对10%、20%、30%??90%的总体测序数据单独进行基因定量分析，并把所有数据条件下得到的基因的表达水平作为最终的数值。用每个百分比条件下求出的单个基因的RPKM数值和最终对应基因的表达水平数值进行比较，如果差异小于15%，则认为这个基因在这个条件下定量是准确的。

图8.1 定量饱和曲线检查分布图

横坐标代表定位到基因组上的reads数占总reads数的百分比，纵坐标代表定量误差在

15%以内的基因的比例

8.2 RNA-Seq相关性检查

生物学重复是任何生物学实验所必须的，高通量测序技术也不例外(Hansen et al.)。生物学重复主要有两个用途：一个是证明所涉及的生物学实验操作是可以重复的且变异不大，另一个为后续的差异基因分析所需要的。样品间基因表达水平相关性是检验实验可靠性和样本选择是否合理性的重要指标。相关系数越接近1，表明样品之间表达模式的相似度越高。Encode计划建议皮尔逊相关系数的平方(R2)大于0.92(理想的取样和实验条件下)。具体的项目操作中，我们要求R2至少要大于0.8，否则需要对样品做出合适的解释，或者重新进行实验。此部

分，我们同时计算了spearman相关系数和kendall-tau相关系数作为参考，这两个主要是针对顺序变量的相关系数，即秩相关。

图8.2 RNA-Seq相关性检查

R^2:pearson相关系数的平方; rho:spearman相关系数; tau:kendall-tau相关系数

8.3 均一性分布检查

理想条件下，对于RNA-seq技术来说，测序序列(reads)之间为独立抽样并且reads在所有表达的转录本上的分布应该呈现均一化分布。然而很多研究表明，很多偏好型的因素都会影响这种均一化的分布(Dohm et al., 2008)。例如，在RNA-seq建库过程中，片段破碎和RNA反转录的顺序不一样会导致RNA-seq最终的数据呈现严重的3’偏好性。其他因素还包括转录区域的GC含量不同、随机引物等等，并且生物体内从5’或者3’的降解过程同样会导致不均一性分布。

图8.3 不同表达水平的转录本的reads密度分布图

High：高表达量转录本；Medium：中度表达量转录本；Low：低表达量转录本；横坐标为

距离转录本5’端的相对位置(以百分比表示)，纵坐标为覆盖深度的平均值

9 基因差异表达分析

9.1 基因表达水平对比

通过所有基因的RPKM的分布图以及盒形图对不同实验条件下的基因表达水平进行比较。对于同一实验条件下的重复样品，最终的RPKM为所有重复数据的平均值。

图9.1 不同实验条件下基因表达水平比对图

RPKM分布图(左图)的横坐标为log10(RPKM), 纵坐标为基因的密度。RPKM盒形图(右图)的横坐标为样品名称，纵坐标为log10(RPKM)，每个区域的盒形图对五个统计量(至上而下分

别为最大值，上四分位数，中值，下四分位数和最小值)

9.2 差异表达基因列表

基因差异表达的输入数据为基因表达水平分析中得到的readcount数据。对于有生物学重复的样品，分析我们采用DESeq（Anders et al, 2010）进行分析： 该分析方法基于的模型是负二项分布，第 i 个基因在第 j 个样本中的 read count 值为Kij，则有

Kij ～ NB(μij,σ

2ij

)

对于无生物学重复的样品，先采用TMM对read count数据进行标准化处理，之后用DEGseq进行差异分析。差异表达基因列表如下：

表9.2 差异基因列表

gene_id readcount_HS readcount_HT 2835.13040784994 3.08090973423561 872.135924228591 143.509885523625 log2FoldChange -9.4224 9.0233 -4.5455 3.295 pval padj Novel05868 4.13191705589116 Novel05608 1603.13964277834 Novel08190 37.3477675965147 Novel05435 1408.60468719178 4.6547e-140 1.139e-135 8.0171e-124 9.8089e-120 1.1332e-40 1.0151e-31 9.2435e-37 6.21e-28 差异基因列表主要包括的内容：

(1) Gene_id: 基因编号

(2) readcount_Sample1：校正后样品1的readcount值 (3) readcount_Sample2：校正后样品2的readcount值 (4) log2FoldChange: log2(Sample1/Sample2) (5) pvalue(pval): 统计学差异显著性检验指标

(6) qvalue(padj): 校正后的pvalue。qvalue越小，表示基因表达差异越显著

9.3 差异表达基因筛选

用火山图可以推断差异基因的整体分布情况，对于无生物学重复的实验，为消除生物学变异，我们从差异倍数和显著水平两个水平进行评估，对差异基因进行筛选，阈值设定一般为: |log2(FoldChange)| > 1 且 qvalue < 0.005。对于有生物学重复的实验，由于DESeq已经进行了生物学变异的消除，我们对差异基因筛选的标准一般为: padj < 0.05。

图9.3 差异基因火山图

有显著性差异表达的基因用红色点表示；横坐标代表基因在不同样本中表达倍数变化；纵

坐标代表基因表达量变化差异的统计学显著性

9.4 差异基因聚类分析

聚类分析用于判断差异基因在不同实验条件下的表达模式；通过将表达模式相同或相近的基因聚集成类，从而识别未知基因的功能或已知基因的未知功能；因为这些同类的基因可能具有相似的功能，或是共同参与同一代谢过程或细胞通路。以不同实验条件下的差异基因的RPKM值为表达水平，做层次聚类(hierarchical clustering)分析，不同的颜色的区域代表不同的聚类分组信息，同组内的基因表达模式相近，可能具有相似的功能或参与相同的生物学过程。 除了差异基因表达量rpkm层次聚类分析，我们对还分别用H-cluster、

K-means和SOM等三种方法对差异基因的相对表达水平值log2(ratios)进行聚类。不同的聚类算法分别将差异基因分为若干cluster，同一cluster中的基因在不同的处理条件下具有相似的表达水平变化趋势。

图6.4 差异基因聚类图

左图为整体rpkm层次聚类图，以log10RPKM值进行聚类，红色表示高表达基因，蓝色表示低表达基因。颜色从红到蓝，表示log10(RPKM)从大到小；右图为log2(ratios)折线图，每个子图中的灰色线条表示一个cluster中的基因在不同实验条件下相对表达量的折线图，蓝色线条表示这个cluster中的所有基因在不同实验条件下相对表达量的平均值的折

线图，x轴表示实验条件，y轴表示相对表达量

本文来源：https://www.bwwdw.com/article/o0xo.html