有参考基因组的转录组生物信息分析Word文档下载推荐.docx

1、1Member of a pair, 1 or 2 （paired-end or mate-pair reads only）YY if the read fails filter （read is bad）, N otherwise180 when none of the control bits are on, otherwise it is an even numberATCACGIndex sequence第四行中每个字符对应的ASCII值减去33，即为对应第二行碱基的测序质量值。如果测序错误率用e表示，illumina HiSeqTM2000/MiSeq的碱基质量值用Qphred表示，

2、则有下列关系：公式一：Qphred=-10log10（e）illumina Casava 1.8版本测序错误率与测序质量值简明对应关系如下：测序错误率测序质量值对应字符5%13.1%2050.1%30?0.01%40I2测序数据质量评估2.1测序错误率分布检查每个碱基测序错误率是通过测序Phred数值（Phred score, Qphred）通过公式1转化得到，而Phred 数值是在碱基识别（Base Calling）过程中通过一种预测碱基判别发生错误概率模型计算得到的，对应关系如下表所显示：illumina Casava 1.8版本碱基识别与Phred分值之间的简明对应关系Phred分值不正

3、确的碱基识别碱基正确识别率Q-sorce101/1090%Q101/10099%Q201/100099.9%Q301/1000099.99%Q40测序错误率与碱基质量有关，受测序仪本身、测序试剂、样品等多个因素共同影响。对于RNA-seq技术，测序错误率分布具有两个特点：（1）测序错误率会随着测序序列（Sequenced Reads）的长度的增加而升高，这是由于测序过程中化学试剂的消耗而导致的，并且为illumina高通量测序平台都具有的特征（Erlich and Mitra, 2008; Jiang et al.）。（2）前6个碱基的位置也会发生较高的测序错误率，而这个长度也正好等于在RNA

4、-seq建库过程中反转录所需要的随机引物的长度。所以推测前6个碱基测序错误率较高的原因为随机引物和RNA模版的不完全结合（Jiang et al.）。测序错误率分布检查用于检测在测序长度范围内，有无异常的碱基位置存在高错误率，比如中间位置的碱基测序错误率显着高于其他位置。一般情况下，每个碱基位置的测序错误率都应该低于0.5%。图2.1测序错误率分布图横坐标为reads的碱基位置，纵坐标为单碱基错误率2.2GC含量分布检查GC含量分布检查用于检测有无AT、GC 分离现象，而这种现象可能是测序或者建库所带来的，并且会影响后续的定量分析。在illumina测序平台的转录组测序中，反转录成cDNA时所

5、用的6bp 的随机引物会引起前几个位置的核苷酸组成存在一定的偏好性。而这种偏好性与测序的物种和实验室环境无关，但会影响转录组测序的均一化程度（Hansen et al.）。除此之外，理论上G和C碱基及A和T碱基含量每个测序循环上应分别相等，且整个测序过程稳定不变，呈水平线。对于DGE测序来说，由于随机引物扩增偏差等原因，常常会导致在测序得到的每个read前6-7个碱基有较大的波动，这种波动属于正常情况。图2.2GC含量分布图横坐标为reads的碱基位置，纵坐标为单碱基所占的比例；不同颜色代表不同的碱基类型2.3测序数据过滤测序得到的原始测序序列，里面含有带接头的、低质量的reads，为了保证信

6、息分析质量，必须对raw reads进行过滤，得到clean reads，后续分析都基于clean reads。数据处理的步骤如下：（1） 去除带接头（adapter）的reads；（2） 去除N（N表示无法确定碱基信息）的比例大于10%的reads；（3） 去除低质量reads。RNA-seq 的接头（Adapter, Oligonucleotide sequences for TruSeqTM RNA and DNA Sample Prep Kits） 信息：RNA 55-AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT-

7、3RNA 35-GATCGGAAGAGCACACGTCTGAACTCCAGTCAC（6位index）ATCTCGTATGCCGTCTTCTGCTTG-3图2.3原始数据过滤结果2.4测序数据质量情况汇总表2.4数据产出质量情况一览表Sample nameRaw readsClean readsclean basesError rate（%）Q20（%）Q30（%）GC content（%）HS1_13.52G0.0397.8892.8849.39HS1_296.5090.3849.59HS2_13.51G97.8592.8149.53数据质量情况详细内容如下：（1） Raw reads：统计原

8、始序列数据，以四行为一个单位，统计每个文件的测序序列的个数。（2） Clean reads：计算方法同 Raw Reads，只是统计的文件为过滤后的测序数据。后续的生物信息分析都是基于Clean reads。（3） Clean bases：测序序列的个数乘以测序序列的长度，并转化为以G为单位。（4） Error rate：通过公式1计算得到。（5） Q20、Q30：分别计算 Phred 数值大于20、30的碱基占总体碱基的百分比。（6） GC content：计算碱基G和C的数量总和占总的碱基数量的百分比。3参考序列比对分析测序序列定位算法：根据不同的基因组的特征，我们选取相对合适的软件（动植

9、物用TopHat（Trapnell et al., 2009）、真菌或者基因密度较高的物种用Bowtie），合适的参数设置（如最大的内含子长度，会根据已知的该物种的基因模型来进行统计分析），将过滤后的测序序列进行基因组定位分析。下图为TopHat的算法示意图： Tophat的算法主要分为两个部分：（1） 将测序序列整段比对到外显子上。（2） 将测序序列分段比对到两个外显子上。我们统计了实验所产生的测序序列的定位个数（Total Mapped Reads）及其占clean reads的百分比，其中包括多个定位的测序序列个数（Multiple Mapped Reads）及其占总体（clean re

10、ads）的百分比，以及单个定位的测序序列个数（Uniquely Mapped Reads）及其占总体（clean reads）的百分比。3.1Reads与参考基因组比对情况统计表3.1Reads与参考基因组比对情况一览表HS1HS2HT1HT2HW1HW2Total readsTotal mappedMultiple mapped606556 （0.86%）633575 （0.9%）714678 （0.94%）450156 （0.89%）389470 （0.84%）335509 （0.83%）Uniquely mappedRead-1Read-2Reads map to +-Non-splic

11、e readsSplice reads9910559 （24.26%）Reads mapped in proper pairs比对结果统计详细内容如下：（1） Total reads：测序序列经过测序数据过滤后的数量统计（Clean data）。（2） Total mapped：能定位到基因组上的测序序列的数量的统计；一般情况下，如果不存在污染并且参考基因组选择合适的情况下，这部分数据的百分比大于 70%。（3） Multiple mapped：在参考序列上有多个比对位置的测序序列的数量统计；这部分数据的百分比一般会小于10%。（4） Uniquely mapped：在参考序列上有唯一比对位置

12、的测序序列的数量统计。（5） Reads map to ，Reads map to ：测序序列比对到基因组上正链和负链的统计。（6） Splice reads：（2）中，分段比对到两个外显子上的测序序列（也称为Junction reads）的统计，Non-splice reads为整段比对到外显子的将测序序列的统计，Splice reads的百分比取决于测序片段的长度。3.2Reads在参考基因组不同区域的分布情况对Total mapped reads的比对到基因组上的各个部分的情况进行统计，定位区域分为Exon（外显子）、Intron（内含子）和Intergenic（基因间隔区域）。正常情况

13、下，Exon （外显子）区域的测序序列定位的百分比含量应该最高，定位到Intron （内含子） 区域的测序序列可能是由于非成熟的mRNA的污染或者基因组注释不完全导致的，而定位到Intergenic（基因间隔区域）的测序序列可能是因为基因组注释不完全以及背景噪音。图3.2Reads在参考基因组不同区域的分布情况3.3Reads在染色体上的密度分布情况对Total mapped reads的比对到基因组上的各个染色体（分正负链）的密度进行统计，如下图所示，具体作图的方法为用滑动窗口（window size）为1K，计算窗口内部比对到碱基位置上的reads的中位数，并转化成 log2 。正常情况下

14、，整个染色体长度越长，该染色体内部定位的reads总数会越多（Marquez et al.）。从定位到染色体上的reads数与染色体长度的关系图中，可以更加直观看出染色体长度和reads总数的关系。图3.3Reads在染色体上的密度分布图上图：横坐标为染色体的长度信息（以百万碱基为单位），纵坐标为log2（reads的密度的中位数），绿色为正链，红色为负链 下图：横坐标为染色体的长度信息（单位为Mb），纵坐标为mapped到染色体上的reads数（单位为M）3.4Reads比对结果可视化我们提供RNA-seq Reads在基因组上比对结果的bam格式文件，部分物种还提供相应的参考基因组和注释文

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

16、浏览器界面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

17、_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,

18、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端剪切（模糊边界） 4.1可变剪切事件分类和数量统计图4.1AS分类和数量统计纵轴为可变剪切事件的分类缩写，横轴为该种事件下可变剪切的数量，不同样品用不同子图和颜色区

19、分4.2可变剪切事件结构和表达量统计表4.2AS结构和表达量统计event_idevent_typegene_idchromevent_startevent_endevent_patternstrandfpkmref_id1000001TSSCUFF.10034382773438330+1.0000000000ENSGALT000000102251000002345021834502533.00000000001000003345674434571652.00000000001000004TTS347480634781785.0000000000（1） event_id: AS事件编号（2）

20、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 exo

21、n; 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 c

22、oordinates of the exon variant）（8） strand: 基因正负链信息（9） fpkm: 此AS类型该基因表达量（10） ref_id: 此基因在参考注释文件中的编号 5新转录本预测将所有测序reads数据的基因组定位结果放到一起，用 Cufflinks 进行组装，然后用Cuffcompare和已知的基因模型进行比较，可以:（1）发现新的未知基因（相对于原有基因注释文件）；（2）发现已知基因新的外显子区域；（3）对已知基因的起始和终止位置进行优化。新基因和新外显子区域预测结果为GTF格式的注释文件。GTF格式的详细说明可参考（）表5.1新转录本结构注释结果seqn

23、amesourcefeaturestartendscoreframeattributesnovelGeneexon1853119499gene_id Novel00001; transcript_id Novel00001.1 exon_number 12081321813Novel00002Novel00002.12391724402Novel00003Novel00003.12518926100Novel00004Novel00004.1（1） seqname：染色体编号（2） source：来源标签，这里的novelGene指新基因（3） feature：区域类型，目前我们预测外显子区域

24、（4） start：起始坐标（5） end：终止坐标（6） score：不必关注（7） strand：正负链信息（8） frame：（9） attributes：属性，包括基因编号、转录本编号等信息表5.2已知基因结构优化Gene_idChromosomeStrandOriginal_spanAssembled_spanENSGALG00000000003ENSGALG00000000004Z-ENSGALG000000000116ENSGALG00000000013222783575278733727835752787453（1） Gene_id：原注释文件中基因命名编号（2） Chromo

25、some：（3） Strand：（4） Original_span：原注释文件中基因起始位置终止位置（5） Assembled_span：转录组拼接结果中基因起始位置终止位置6SNP和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，并对原始结果进行过滤，得到如下表形式的分析结果。