RNA-Seq测序数据分析服务流程(试运行)(精)

1、北京大学生科院/CLS生物信息平台RNA-Seq测序数据分析服务流程（试运行2015.3平台联系人：李程（lch3000文档撰写：张超Table of Contents1 .测序质量评估（31.1 测序数据过滤（31.2 质量值分布（31.3 GC含量分布（42 .参考序列比对（43 .基因表达水平（63.1 基因表达水平定量（63.2 基因表达水平分步（63.3 生物学重复相关性分析（63.4 样本间层次聚类及PCA分析（74 .差异基因分析（74.1 基因表达标准化（74.2 差异基因列表（84.3 差异基因可视化（84.4 差异基因聚类（95 .差异表达基因功能分析（105.1 GO富集

2、分析（105.2 信号通路富集分析（105.3 癌基因功能注释（116 .基因结构差异分析（116.1 可变剪切分析（117 . SNP 分析（127.1 SNP 检测（127.2 SNP 筛选（127.3 GO/KEGG 富集（121 .测序质量评估通过测序的数据进行进行质控，保证数据质量适合下游分析。这里我们使用 fastqc和RNA-SeQC来对数据进行质量评定。1.1 测序数据过滤测序得到的原始下机数据往往有许多问题，不能直接使用，通常会经过以下过滤 尽量保证测序数据的质量。a去除带测序接头的测序序列（reads;b.去除低质重的reads1.2 质量值分布 按照现有的测序技术（ill

3、umina平台单碱基的错误率应控制在1%以下，即质量值在20以上。q it 0 8 鼻 I bn- I fwigw I kuni- 1 S it1 tiding I4 二.1*5及l m 3 *、n 了引slh 6 心 x / 熊，州型力函k m h 及 n re w jo v a 阿 用 酊便由6%陇nd Ibp横坐标为reads的碱基位置，纵坐标为单碱基质量值质量值与错误率的关系：Q=-10log10（e;其中Q phred为测序碱基质量值，e为测phred序错误率。1.3 GC含量分布对于RNA测序,鉴于序列通过超声随机打断，所以理论上每个测序循环上的C、 G及A、T含量应分布相等，并且

4、CG-content对于每个物种应大致相同。横坐标为reads的碱基位置，纵坐标为各种碱基的不同比例2 .参考序列比对对于通过质量控制的数据，可以进行后续分析。首先需要将 clean reads比对到 参考基因组上。由于测序时reads是随机的，只有这些reads的碱基信息和质量信息 没有其在基因组上的位置信息，比对这一步就是给所有reads一个在基因组上位置的 信息。在RNA测序中，其实测的是cDNA的序列，由于内含子的存在，所以会较常出现 一条read跨内含子的情况，tophat2可以较好的处理这种情况，所以我彳门选用tophat2 来做比对。比对率间接反应了测序的质量和建库的质量，若比对

5、率低，很可能建库时混入了 其他物种的序列，导致无法比对到研究的物种参考基因组上。Left reads:Input: 11607353Mapped; 11607353of input)q网电(g.S%】 have nultiple aligiments82 howe 26) Rightinput i11607353Hipped :11607352(1加.被of input)of ttese:60103(have nuttplealtgimentshcnreads比对到基因上的位置统计Sample IntragenicRateExonicRateIntronicRateIntergenicRate

6、SplitReadsExpressionProfilingEfficiencyTranscriptsDetectedGenesDetected1BJ 0.885 0.738 0.147 0.114 9,910,010 0.738 32,796 15,434(Sample才羊本名(2IntragenicRate比对到基因内的reads比例(3ExonicRate:比对至U外显子的reads比例(4IntronicRate:比对到内含子的reads比例(5IntergenicRate比对到基因间区的reads比例(6SplitReads:比对至U两外显子交接处的 reads数(7Expressio

7、nProfilingEfficiency:比对到外显子上的reads占总体的比例(8TranscriptsDetected比对上reads数大于5的转录本数(9GenesDetectect匕对上reads数大于5的基因数Left reads;InpMtt； 11 即7353Mapped: 11667353 aee+0K of input)o6010 ( e,SMJ have nultiple olignnertts Cfi2 Iwe 20 Right reads: input ： 11687353Mapped : 11 瑰7352 Q掘I 限 of iwO of tfww; 60103 ( A

8、L5S0 have nulttple a1 ignments Cq haw* 2) IMOK owrall isud nqppin$ rate.AligMd pdifE： L1W3UoF tK&e:601B3 C 0泉D have raultile oligmerrt上8 C or* d1 scorttont atignr帕nt型 lOO.fl% concordant pair allQ-nent rate.3 .基因表达水平3.1 基因表达水平定量在RNA-seq分析中，我们可以通过定位到基因组区域或基因外显子区的reads的计数来估计基因的表达水平。Reads计数除了与基因的真实表达水平成

9、正比外，还 与基因的长度和测序深度成正相关。为了使不同基因、不同实验间估计的基因表达 水平具有可比性，人们引入了 RPKM的概念,RPKM(Reads Per Kilo bases per Million reads是每百万reads中来自某一基因每千碱基长度的 reads数目。RPKM同时考虑 了测序深度和基因长度对reads计数的影响，是目前最为常用的基因表达水平估算方 法(Mortazavi et al., 2008。Gene_ID Sample1 Sample2 Sample3 Sample4 Sample5 Sample6 ENSG00000000003 49.32 46.94 48

10、.91 22.51 20.60 22.95 ENSG00000000419 35.92 34.58 33.69 32.80 35.65 32.73 ENSG00000000457 1.34 0.94 1.19 2.06 2.13 2.26ENSG00000000460 1.19 1.20 1.22 3.00 3.33 3.06(1 Gene_ID:Ensembl基因 ID(2 Other columns各样本中该基因的表达水平(RPKM3.2 基因表达水平分步每个样本所有基因的RPKM盒形图可以展示出不同实验条件下基因表达水平 的分布情况。图3.2.1不同条件下的基因表达水平分布图3.3 生物

11、学重复相关性分析生物学重复主要有两个用途：一个是证明所涉及的生物学实验可重复性强、差 异小，另一个用于估计生物学变异进行差异基因检测。样品间基因表达水平相关性是检验实验可靠性和样本选择是否合理的重要指标。相关系数越接 近1,表明样品之间表达模式的相似度越高。1 EDon EK图3.3.1生物学重复散点图3.4 样本间层次聚类及PCA分析当样本数目较多时，可以利用基因的表达量进行样本间聚类分析及 PCA分析，对 样本间关系进行探究或者对实验设计进行验证。样本聚类距离或者 PCA距离越近， 说明样本越相似。4 .差异基因分析4.1 基因表达标准化对于有生物学重复的样品，我们采用DESeq2提出的s

12、caling factor的方法对原始 的readcount行标准化（normalization。以消除非生物学引起的readcount的差异（最主要消除各个文库测序数据量不同带来的差异。对于标 准化的结果，我们采用MA-plot或box-plot来评价。图4.1.1 MA-plot横坐标为表达量，纵坐标为log后的表达差异倍数基于大部分基因都是非差异表达的，所以大多点应在log fold change=0左右，并 且不随表达量的变化而变化。4.2 差异基因列表对于有生物学重复的的样品，我们采用DESeq2来分析差异表达基因。该方法 基于负二项分布模型(K ij NB( pij,来检测差异表达

13、基因。Gene baseMean log2FoldChange pvalue padjFBgn0000370 31324.379200 -1.3665378519 5.6393206e-176 2.9843284e-OI6J 印1psfi曾ie Gi tc4Ci lerOJ tetO5ftiean e(ortssion172FBgn0033913 17544.483454 -1.1571536021 6.3177309e-90 1.3373372e-87(1Gene:基因 ID(2baseMean所有样本矫正后的平均reads数(3log2FoldChange:log2后的表达量差异(4pva

14、lue统计学差异显著性检验指标(5padj:校正后的pvalue。padj越小，表示基因表达差异越显著4.3 差异基因可视化用火山图可以推断差异基因的整体分布情况。图4.3.1显著性差异表达基因用红色点表示；横坐标表示基因在不同样本中的表达倍数变化；纵坐标表示统计学上基因表达量变化差异的显著性对于特定基因在不同实验中的表达情况，和此基因的不同转录本在不同实验中 的表达情况中白，U七1图4.3.2左图为regucalcin基因在两个样本中的表达差异情况；右图为此基因在不同转录本中的表达差异情况4.4 差异基因聚类聚类分析用于判断差异基因在不同实验条件下的表达模式。通过将表达模式相 同或相近的基因

15、聚集成类，从而识别未知基因的功能或已知基因的未知功能。5 .差异表达基因功能分析5.1 GO富集分析Gene Ontology箭称 GO, /是基因功能标准分类体 系。研究差异基因在Gene Ontology中的分布情况将阐明差异基因富集的生物学功 能。M口, Z Scwfl5.2 信号通路富集分析在生物体内，不同基因相互协调实现其生物学功能，通过Pathway显著性富集能 确定差异表达基因参与的最主要信号通路。KEGG(Kyoto Encyclopedia of Genes andGenomes,http:/www.kegg.jp层有关 P

16、athway的主要公共数据库(Kanehisa,200& Pathway显著性富集分析以KEGG Pathway为单位，应用统计检验找出差异表达基因 中显著性富集的Pathway。Color Ker啊dl HSg.raip5.3 癌基因功能注释原癌基因(Proto-oncogene是参与细胞生长、细胞分裂和细 胞分化的正 常基因，当其发生突变后(如基因序列被改变就会变成致癌基因 (Oncogene通常在月中瘤或恶性细胞系中某些特异性癌基因会上调表达，通过了解 癌基因在实验不同组的表达情况有助于深入认识疾病的发病机理。Cosmic( https:cancer.sanger.ac.uk/cosmi

17、是英国 Sanger实验室开发并维护的癌基 因及相关注释数据库，有较高的权威性及可信度，通过与数据库进行比对，可对差异表达基因中的癌基因部分进行鉴别和注释。6.基因结构差异分析6.1可变剪切分析对于RNA-seq,除了 gene水平的差异分析外，还可以进行 exon水平的差异 分析。不用的exon表达，表明了有着不同的剪切方式。这时可以使用 Bioconductor的DEXSeq软件包。该分析可以给出每个基因在不同的实验条件下， 外显子的使用情况。比 如上图的10号外显子在control和knockdown两组中的表 达差别较大，此外显子的表达量情况，也反映到了在两组中此基因的剪切形式有 差异。7. SNP 分析 7.1 SNP 检测 SNP 全称 Single Nucleotide Polymorphisms.,是指在 基因组上由单个核甘酸变异形成的遗传标记，其数量很多，多态性丰富