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

亠、生物信息分析流程获得原始测序序列(Sequeneed Reads)后，在有相关物种参考序列或参考 基因组的情况下，通过如下流程进行生物信息分析:原始测序序別测序数据质量评佶切娈剪功分析 新转录△预测1、项目结果说明1原始序列数据高通量测序(如illumina HiSeq TM 2000/MiSeq 等测序平台)测序得到的原 始图像数据文件经碱基识别(Base Calling)分析转化为原始测序序列(Sequeneed 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测序标识符(SequeneeIdentifiers)和描述文字(选择性部分)；第二行是碱基序列；第三行以“ +”开头, 随后为illumina 测序标识符(选择性部分);第四行是对应序列的测序质量(Cockr---------------------------------------、RNA-£E 口整体质量评估 基因差异表达分桁蛋口网络互作分析 k ____________ ______ )GO 富集分析KEGCg 集分析et al.)。

illumi na 测序标识符详细信息如下：EAS139 Uni que in strume nt n ame136 Run IDFC706VJ Flowcell ID2 Flowcell la ne2104 Tile n umber with in the flowcell la ne15343 'x'-coord in ate of the cluster within the tile197393 'y'-coordi nate of the cluster within the tile1 Member of a pair, 1 or2 (paired-e nd or mate-pair reads only)Y Y if the read fails filter (read is bad), N otherwise18 0 when none of the control bits are on, otherwise it is an even number ATCACG In dex seque nee第四行中每个字符对应的ASCII值减去33 ,即为对应第二行碱基的测序质量值。

如果测序错误率用 e表示，illumina HiSeq TM2000/MiSeq 的碱基质量值用Q phred表示，则有下列关系：公式一：Q phred = -10lOg 10(e)illumi na Casava 1.8 版本测序错误率与测序质量值简明对应关系如下:测序错误率测序质量值对应字符5% 131% 20 50.1% 30 ?0.01% 40 I2测序数据质量评估2.1 测序错误率分布检查每个碱基测序错误率是通过测序 Phred数值(Phred score, Q phred)通过公式 1转化得到，而Phred数值是在碱基识别(Base Calling)过程中通过一种预测碱基判别发生错误概率模型计算得到的，对应关系如下表所显示：illumina Casava 1.8 版本碱基识别与 Phred分值之间的简明对应关系Phred分值不正确的碱基识别碱基正确识别率Q-sorce10 1/10 90% Q1020 1/100 99% Q2030 1/1000 99.9% Q3040 1/10000 99.99% Q40测序错误率与碱基质量有关，受测序仪本身、测序试剂、样品等多个因素共同影响。

对于RNA-seq技术，测序错误率分布具有两个特点：(1) 测序错误率会随着测序序列(Sequeneed Reads)的长度的增加而升高，这是由于测序过程中化学试剂的消耗而导致的，并且为illumi na 高通量测序平台都具有的特征(Erlich and Mitra, 2008; Jia ng et al.) 。

(2) 前6个碱基的位置也会发生较高的测序错误率，而这个长度也正好等于在RNA-seq建库过程中反转录所需要的随机引物的长度。

所以推测前6个碱基测序错误率较高的原因为随机引物和 RNA模版的不完全结合(Jiang et al.)。

测序错误率分布检查用于检测在测序长度范围内，有无异常的碱基位置存在高错误率，比如中间位置的碱基测序错误率显著高于其他位置。

一般情况下，每个碱基位置的测序错误率都应该低于 0.5%。

图2.1 测序错误率分布图横坐标为reads的碱基位置，纵坐标为单碱基错误率Error rate dislribuLion along reads (HSl JPosiLorr alortg reads2.2 GC 含量分布检查GC 含量分布检查用于检测有无 AT 、GC 分离现象，而这种现象可能是测序或者建库所带来的，并且会影响后续的定量分析。

在 illumina 测序平台的转录组测序中，反转录成 cDNA 时所用的 6bp 的随机引物会引起前几个位置的核苷酸组成存在一定的偏好性。

而这种偏好性与测序的物种和实验室环境无关，但会影响转录组测序的均一化程度 (Hansen et al.) 。

除此之外，理论上 G 和 C 碱基及 A 和 T 碱基含量每个测序循环上应分别相等，且整个测序过程稳定不变，呈水平线。

对于 DGE 测序来说，由于随机引物扩增偏差等原因，常常会导致在测序得到的每个 read 前 6-7 个碱基有较大的波动，这种波动属于正常情况。

T 11111II1 1fl 4060 00 WO 120 140 160 18CPosition along reads图2.2 GC 含量分布图横坐标为reads 的碱基位置，纵坐标为单碱基所占的比例；不同颜色代表不同的碱基类型2.3 测序数据过滤测序得到的原始测序序列，里面含有带接头的、低质量的 reads ，为了保证 信息分析质量，必须对raw reads 进行过滤，得到clean reads ，后续分析都 基于 clean readsBases content along reads (HS1)O■2窗BqO C ①olo d数据处理的步骤如下:⑴去除带接头（adapter）的reads ；⑵去除N（N表示无法确定碱基信息）的比例大于10%的reads ;⑶去除低质量reads。

TM RNA-seq 的接头(Adapter, Oligonucleotide sequences for TruSeqRNA and DNA Sample Prep Kits) 信息:RNA 5 ' Adapter (RA5), part # 150132055' -AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT- 3'RNA 3 ' Adapter (RA3), part # 15013207 :5' -GATCGGAAGAGCACACGTCTGAACTCCAGTCAC(6 _位index)ATCTCGTATGCCGTCTTCTGCTTG- 3'Classification of Raw Reads (HS1)图2.3 原始数据过滤结果2.4 测序数据质量情况汇总表2.4 数据产出质量情况一览表Clean Reads [35175205. 96.16%) Containing N (15020, 0.04%) Low Qualrty (1212560. 3 31 %) Adapter Related (176803, 0.48%)SampleRaw readsClean n ame reads clea n Error GCQ20(%) Q30(%)bases rate(%) con te nt(%)HS1_ _1 36579608 35175205 3.52G 0.03 97.88 92.88 49.39 HS1_ _2 36579608 35175205 3.52G 0.03 96.50 90.38 49.59 HS2 1 36547734 35119463 3.51G 0.03 97.85 92.81 49.53数据质量情况详细内容如下:(1) Raw reads :统计原始序列数据，以四行为一个单位，统计每个文件的测序序列的个数。

⑵Clean reads :计算方法同 Raw Reads，只是统计的文件为过滤后的测序数据。

后续的生物信息分析都是基于Clean reads。

⑶Clean bases :测序序列的个数乘以测序序列的长度，并转化为以G为单位。

⑷Error rate :通过公式1计算得到。

⑸Q20、Q30 :分别计算Phred数值大于20、30的碱基占总体碱基的百分比。

⑹GC content :计算碱基G和C的数量总和占总的碱基数量的百分比。

3参考序列比对分析测序序列定位算法：根据不同的基因组的特征，我们选取相对合适的软件(动植物用TopHat(Trapnell et al., 2009) 、真菌或者基因密度较高的物种用Bowtie)，合适的参数设置(如最大的内含子长度，会根据已知的该物种的基因模型来进行统计分析)，将过滤后的测序序列进行基因组定位分析。

下图为TopHatTophat 的算法主要分为两个部分：(1) 将测序序列整段比对到外显子上。

(2) 将测序序列分段比对到两个外显子上。

我们统计了实验所产生的测序序列的定位个数 (Total Map ped Reads) 及其 占clean reads 的百分比，其中包括多个定位的测序序列个数(MultipleMapped Reads)及其占总体(clean reads )的百分比，以及单个定位的测序序 列个数(Uniquely Mapped Reads) 及其占总体(clean reads )的百分比。

3.1 Reads 与参考基因组比对情况统计的算法示意图:■严生于E SL 的測序序列un 产生于测序=＞氐中1,E “胆的连按区壇測序片段将■测序序列81段比対到?卜显子上(氐彌Hud Uappin^ )将测序序列分段比刘到两个外显子上(Juncticin Keads Mapp me )TopHat 外显子优先定位算法Eson.-flist MappingRNA■ KI 口表3.1 Reads与参考基因组比对情况一览表Sample n ame HS1 HS2 HT1 HT2 HW1 HW2Total reads 70350410 70238926 76161678 50666084 46573662 4054311860529821 60232484 63555439 43461327 40246848 34971284Total mapped(86.04%) (85.75%) (83.45%) (85.78%) (86.42%) (86.26%)Multiple 606556 633575 714678 450156 389470 335509mapped (0.86%) (0.9%) (0.94%) (0.89%) (0.84%) (0.83%)Uni quely 59923265 59598909 62840761 43011171 39857378 34635775mapped (85.18%) (84.85%) (82.51%) (84.89%) (85.58%) (85.37%) 30176973 29987004 31592931 21654629 20028779 17411209 Read-1(42.9%) (42.69%) (41.48%) (42.74%) (43%) (43.02%)29746292 29611905 31247830 21356542 19828599 17224566 Read-2(42.28%) (42.16%) (41.03%) (42.15%) (42.57%) (42.35%)Reads map to 29930036 29783311 31409912 21476601 19923501 17289330 '+' (42.54%) (42.4%) (41.24%) (42.39%) (42.78%) (42.61%)Reads map to 29993229 29815598 31430849 21534570 19933877 173464451 1 (42.63%) (42.45%) (41.27%) (42.5%) (42.8%) (42.76%)Non-splice 42357242 42528691 45227757 31347392 28062847 24725216 reads (60.21%) (60.55%) (59.38%) (61.87%) (60.25%) (61.1%)17566023 17070218 17613004 11663779 11794531 9910559Splice reads(24.97%) (24.3%) (23.13%) (23.02%) (25.32%) (24.26%)Reads mapped 53795182 54428240 56181352 38524314 36101400 31246362in proper(76.47%) (77.49%) (73.77%) (76.04%) (77.51%) (77.25%)pairs比对结果统计详细内容如下:(1) Total reads :测序序列经过测序数据过滤后的数量统计 (Clean data) <⑵Total map ped :能定位到基因组上的测序序列的数量的统计；一般情况下，如果不存在污染并且参考基因组选择合适的情况下，这部分数据的百分比大于70%。