诺禾致源dge生物信息分析结题报告

DGE生物信息分析结题报一、建库流TotalRNA上机差异GO富集分差异KEGG富集分差异蛋白互作网络分文 列结题报告PDF一、建库流从RNA样品到最终数据获得，样品检测、建库、每一个环节都会对数据质量和数量产生析正确、全面、可信的前提。为了从上保证数据的准确性、可靠性，诺禾致源对样品检测、建库、每一个生产步骤都严格把控，从根本上确保了高质量数据的产出。流程图如下TotalRNAAgilent2100精确检测RNA样品检测合格后，用带有Oligo（dT）的磁珠富集真核生物mRNA（若为原核生物，则通过试剂盒去除rRNA来富集mRNA）。随后加入fragmentationbuffer将mRNA打断成短片段，以mRNA板，用六碱基随机引物（randomhexamers）合成一链cDNA，然后加入缓冲液、dNTPs和DNApolymeraseI合成二链cDNA，随后利用AMPureXPbeads纯化双链cDNA。纯化的双链cDNA再进行末端修复、加A尾并连接接头，然后用AMPureXPbeads进行片段大小选择，最后进行PCR富集得到最终的cDNA文库。构建原理图如下：2100对文库的insertsize进行检测，insertsize符合预期后，使用Q-PCR方法对文库的有效浓度进行准确定量（文库有效浓度＞2nM），以保证文库质量。上机二、生物信息分析流程获得原始序列(SequencedReads)后，在有相关物种参考序列或参考组的情况下，通过如程进行生物信息分析：三、项目结果说明1高通量(如illuminaHiSeqTM2000/MiSeq等平台)得到的原始图像数据文件经碱基识别(BaseCalling)分析转化为原始序列(SequencedReads)，我们称之为RawData或RawReads，结果以FASTQ(简称为fq)文件格式，其中包含序列(reads)的序列信息以及其对应的质量信息。@EAS139:136:FC706VJ:2:2104:15343:1973931:Y:18:ATCACG其中第一行以“@”开头，随后为illumina标识符(SequenceIdentifiers)和描述文字(选择性部分)；第二行是碱基序列；第三行以“+”开头，随后为illumina标识符(选择性部分)；第四行是对应序列的质量(Cocketal.)。illumina标识符详细信息如下UniqueinstrumentRunFlowcell2FlowcellTilenumberwithintheflowcell'x'-coordinateoftheclusterwithinthe'y'-coordinateoftheclusterwithinthe1Memberofapair,1or2(paired-endormate-pairreadsYYifthereadfailsfilter(readisbad),N0whennoneofthecontrolbitsareon,otherwiseitisanevenIndex误率用e表示，illuminaHiSeqTM2000/MiSeq的碱基质量值用Qphred表示，则有下列关系：一： =-illuminaCasava1.8版本错误率与质量值简明对应关系如下错误率质量值对应字符.5?I2数据质量评 错误率分布检每个碱基错误率是通过Phred数值(Phredscore,Qphred)通过1转化得到，而Phred(BaseCalling)过程中通过一种预测碱基判别发生错误概率模型计算得illuminaCasava1.8版本碱基识别与PhredPhred分不正确的碱基识别碱基正确识别率Q-错误率与碱基质量有关，受仪本身、试剂、样品等多个因素共同影响。对于RNA-seq技术，错误率分布具有两个特点：(1)错误率会随着序列(Sequenced的原因为随机引物和RNA模版的不完全结合(Jiangetal.)。图2.1错误率分布横坐标为reads的碱基位置，纵坐标为单碱基错GCGC含量分布检查用于检测有无AT、GC分离现象，而这种现象可能是或者建库所带来在illumina平台的转录组中，反转录成cDNA时所用的6bp的随机引物会引起前几个位置的核苷酸组成存在一定的偏。而这种偏与的物种和环境无关，但会影响转录组的均一化程度(Hansenetal.)。除此之外，理论上G和C碱基及A和T碱基含量每个测序循环上应分别相等，且整个过程稳定不变，呈水平线。对于DGE来说，由于随机引物扩增偏差等原因，常常会导致在得到的每个read前6-7个碱基有较大的波动，这种波动属于图2.2GC横坐标为reads的碱基位置，纵坐标为单碱基所占的比例；不同颜色代表不同的碱基类 数据过必须对rawreads进行过滤，得到cleanreads，后续分析都基于cleanreads。RNA-seq(Adapter,OligonucleotidesequencesforTruSeqTMRNAandDNASamplePrepKits)信息：RNA5’Adapter(RA5),part RNA3’Adapter(RA3),part 图2.3原始数据过滤结果数据质量情况汇总

表2.4数据产出质量情况一览表SampleRawCleancleanErrorGCsamplesampleRawreads：统计原始序列数据，以四行为一个单位，统计每个文件的序列的个Cleanreads：计算方法同RawReads，只是统计的文件为过滤后的数据。后续的生物信息分析都是基于Cleanreads。Cleanbases：序列的个数乘以序列的长度，并转化为以G为单位Errorrate：通过1计算得到Q20、Q30：分别计算Phred数值大于20、30GCcontent：计算碱基G和C3序列定位算法：根据不同的组的特征，我们选取相对合适的软件(动植物用TopHat2、细菌或者密度较高的物种用Bowtie2)，合适的参数设置(如最大的内含子长度，会根据已知的该物种的模型来进行统计分析)，将过滤后的序列进行组定位分析。下图为TopHat2的算法示意图：将序列和转录组进行比对(可选将序列整段比对到组外显子将序列分段比对到组的两个外显子比正常情况下会高于70%(TotalMappedReadsorFragments)，其中具有多个定位的序列(MultipleMappedReadsorFragments)占总体的百分比通常不会超过10%。Reads与参考组比对情况统表3.1Reads与参考组比对情况一览SamplesamplesampleTotalTotalMultiple716688UniquelyReadsmaptoReadsmapto'-Non-spliceSpliceTotalreads：序列经过数据过滤后的数量统计(Cleandata)Totalmapped：能定位到组上的序列的数量的统计；一般情况下，如果不存在污染并且参考组选择合适的情况下，这部分数据的百分比大于70%。Multiplemapped：在参考序列上有多个比对位置的序列的数量统计；这部分数据Uniquelymapped：在参考序列上有唯一比对位置的序列的数量统计Readsmapto'+'，Readsmapto'-'：序列比对到组上正链和负链的统计Splicereads：(2)中，分段比对到两个外显子上的序列(也称为Junctionreads)的统计，Non-splicereads为整段比对到外显子的将序列的统计，Splicereads的百分 Reads在参考组不同区域的分布情对Totalmappedreads的比对到组上的各个部分的情况进行统计，定位区域分正常情况下，Exon(外显子)区域的序列定位的百分比含量应该最高，定位到Intron(内含子)区域的序列可能是由于非成mRNA的污染或者组注释不完全导致的，而定位到Intergenic(间隔区域)的序列可能是因为组注释不完全以及背景噪音。图3.2Reads在参考组不同区域的分布情Reads在上的密度分布情对Totalmappedreads的比对到组上的各个（分正负链）的密度进行统计，如下图所示，具体作图的方法为用滑动窗口(windowsize)为1K，计算窗口内部比对到碱基位置上的reads的中位数，并转化成log2。正常情况下，整个长度越长，该内部定位的reads总数会越多(Marquezetal.)。从定位到上的reads数与长度的关系图中，可以更加直出长度和reads总数的关系。图3.3Reads在上的密度分布图一：横坐为的度信息(以百碱基为单)纵坐标为log(reads的密的中位)，色为正链，红色为链图二：横标为长度信(单Mb)，纵标mapped到上reads单位为M)Reads我们提供RNA-seqReads在组上比对结果的bam格式文件，部分物种还提供相应的参考基因组和注释文件，并推荐使用IGV(IntegrativeGenomicsViewer)浏览器对bam文件进行可视化浏览。IGV浏览器具有以下特点：(1)能在不同尺度下显示单个或多个读段在组上的位置，包括读段在各个上的分布情况和在注释的外显子、内含子、剪接接合区、间区的分布及其剪接异构体的注释信息；(4)能显示其他注释信息；(5)既可以从服务器端各种图3.4IGV4表达水平分一个表达水平的直接体现就是其转录本的丰度情况，转录本丰度程度越高，则表达水平越高。在RNA-seq分析中，我们可以通过定位到组区域或外显子区的序列(reads)的计数来估计的表达水平。Reads计数除了与的真实表达水平成正比外，还与基因的长度和深度成正相关。为了使不同、不同实验间估计的表达水平具有可比性，人们引入了RPKM的概念，RPKM(ReadsPerKilobasesperMillionreads)是每百万reads中来自某一每千碱基长度的reads数目。RPKM同时考虑了深度和长度对reads计数的影响，是目前最为常用的表达水平估算方法(Mortazavietal.,2008)表4.1不同表达水平区间的数量统计RPKMsamplesample0-1-3-15-表4.2表达水平统计samplesampleRNA-seq表达水平的饱和曲线检查被准确定量；反之，表达量低的，需要较大的数据量才能被准确定量。比条件下求出的单个的RPKM数值和最终对应的表达水平数值进行比较，如果差异小于15%，则认为这个在这个条件下定量是准确的。图5.1定量饱和曲线检查分布图横坐标代表定位 组上的reads数占总reads数的百分比，纵坐标代表定量误差在15%以内 的比RNA-Seq生物学重复是任何生物学实验所必须的，高通量技术也不例外(Hansenetal.)。生物为后续的差异分析所需要的。样品间表达水平相关性是检验实验可靠性和样本选择是否图5.2RNA-Seq图一：样品间的相关系数散点图，R2:pearson相关系数的平方；图二：样品间相关系数热理想条件下，对于RNA-seq技术来说，序列(reads)之间为独立抽样并且reads在所有表均一化的分布(Dohmetal.,2008)。例如，在RNA-seq建库过程中，片段破碎和RNA反转录的顺序不一样会导致RNA-seq最终的数据呈现严重的3’偏。其他因素还包括转录区域的GC含量不图5.3不同表达水平的转录本的readsHighMediumLow56差异表达分 表达水平对图6.1不同实验条件下表达水平比对RPKM分布图图一)的横坐log10(RPKM),纵坐标为的密度。RPKM盒形图(二)的横坐标样品名称，纵坐标为log10(RPKM差异表达列的样品，分析我们采用DESeq（Andersetal,2010）进行分析：该分析方法基于的模型是负二项分布，第i个在第j个样本中的readcount值Kij，则Kij～对于无生物学重复的样品，先采用TMM对readcount数据进行标准化处理，之后用DEGseq进表6.2差异列-5.1608e-3.7539e-4.1443e-2.0931e-2.4452e--Gene_id:编log2FoldChange:pvalue(pval):qvalue(padj):校正后的pvalue。qvalue越小，表 差异表达筛异，我们从差异倍数和显著水平两个水平进行评估，对差异进行筛选，阈值设定一般为:|log2(FoldChange)|1qvalue0.005。对于有生物学重复的实验，由于DESeq已经进行了生物学变异的消除，我们对差异筛选的标准一般为:padj<0.05。图6.3差异火山的表表差异聚类分因成类，从而识别未知的功能或已知的未知功能；因为这些同类的可能具有相达水平，做层次聚类(hierarchicalclustering)分析，不同的颜色的区域代表不同的聚类分组信息，同组内的表达模式相近，可能具有相似的功能或参与相同的生物学过程。除了差异表达量rpkm层次聚类分析，我们对还分别用H-cluster、K-means和SOM等三种方法对差异的相对表达水平值log2(ratios)进行聚类。不同的聚类算法分别将差异分为若干cluster，同一cluster中的在不同的处理条件下具有相似的表达水平变化趋势。图6.4差异聚类图一为整体rpkm层次聚类，log10RPKM达达红到蓝，表示log10(RPKMlog2(ratioscluster中的在同实验条件相对表达量折线图，蓝线条表示这个cluster中的有在不实验条件下xy差异GO富集分GeneOntology（简称GO,）是功能国际标准分类体系。根据实验目的筛选差异后，研究差异在GeneOntology中的分布状况将阐明实验中样本差异在功能上的体现。GO富集分析方法为GOseq（Youngetal,2010）,此方法基于Walleniusnon-centralhyper-geometricdistribution。相对于普通的超几何分布(Hyper-geometricdistribution)，此分布的特点是从某个类别中抽取的概率与从某个类别之外抽取一个的概率是不同的，这种概率的不同是通过对长度的偏进行估计得到的，从而能更为准确地计算出GOterm被差异富集的概率。差异GO富集列表7.1差异GO富集列metabolicprocess2.3219e-8.96e-2.8685e-5.4238e-4.2165e-5.4238e-1.6237e-1.2532e-GO_accession：GeneOntologyDescription：GeneOntology Corrected_pValue：矫正后的P-Value，一般情况下，P-value<0.05DEG_item：与该GO相关的差异的数DEG_list：GO注释的差异数差异GO富集DAG有向无环图(DirectedAcyclicGraph，DAG)为差异GO富集分析结果的图形化展示方位作为有向无环图的主节点，并通过包含关系，将相关联的GOTerm一起展示，颜色的深浅代表富集程度。我们的项目中分别绘制生物过程(biologicalprocess)、分子功能(molecularfunction)和细胞组分(cellularcomponent)的DAG图。图7.2GOGOTOP10的GOTERMp-value差异GO富集柱状差异GO富集柱状图，直观的反映出在生物过程(biologicalprocess)、细胞组分(cellularcomponent)和分子功能(molecularfunction)富集的GOterm上差异的个数分布情况。我们挑选了富集最显著的30个GOterm在图中展示，如果不足30条，则全部展示。图7.3GO图一：纵坐标为富集的GOterm，横坐标为该term中差异个数。不同颜色用来区分生物过程、细胞组分和分子功能，带“*”为富集的GOterm图二：对图一中的GO，按生物过程、细胞组分和分子功能三大类别及差差异KEGG富集分在生物体内，不同相互协调行使其生物学功能，通过Pathway显著性富集能确定差异表达参与的最主要生化代谢途径和信号转导途径。KEGG（KyotoEncyclopediaofGenesandPathway为单位，应用超几何检验，找出与整个组背景相比，在差异表达中显著性富集差异KEGG富集列表8.1差异KEGG富集列P-CorrectedP- e- e-Glutathione e- Samplenumber：该通路下差异的个数Backgroundnumber：该通路下的个数CorrectedP-value：矫正后的统计学显著水平，一般情况下，P-value0.05差异KEGG富集散点factor、Qvalue和富集到此通的个数来衡量。其中Richfactor指差异表达的中位于该pathway条目的数目与所有有注释中位于该pathway条目的总数的比值。Rich图8.2差异KEGG富集散点纵轴表示pathway名称，横轴表示Richfactor，点的大小表示此pathway中差异表达个数多少，而点的富集KEGG标色同上，括号中数字为log2(Foldchange)。以上步骤可脱机实现，如连接互联网，点击各个节点，可以连接到KEGG数据库中各个KO的具体信息页。图8.3显著富集的KEGGpathway差异蛋白互作网络分析我们主要应用STRING蛋白质互作数据库()中的互作关系，针对数据库中包含的物种，直接从数据库中提取出目标集(比如差异list)的互作关系构建网络；针对数据库中不包含的物种，我们首先将目标集序列应用blastx比对到string数据库中包含我们提供差异蛋白互作网络数据文件，此文件可以直接导入Cytoscape软件进行可视化编辑。Cytoscape软件使用方法可参考我们提供的使用说明文档(CytoscapeQuickStart.pdf户可以针对一些网络的拓扑属性进行统计和标示作图，比如：互作网络图点(node)此节点的度(degree)成正比，即与此节点相连的边越多，它的度越大，节点也就越大，这些节点在网络中可能处于较为的位置。节点的颜色与此节点的系数(clustering相关，颜色梯度由绿到红对应系数的值由低到高；系数表示此节点的邻接点之间的连通性好坏，系数值越高表示此节点的邻接点之间的连通性越好等等。根据不同的研究目的和需求，客户还可以在网络图中进行调整节点位置和颜色、标注表达量水操作。需要注意的是，通过blast比对得到的结果不能保证较好的准确性，这部分的工作只是给客户提供参考，辅助客户发现一些可能的重要的。按我们提供的使用说明将文件导入Cytoscape软件后的效果图如下：图9Cytoscape四、参考文献Anders,S.(2010).HTSeq:ysinghigh-throughputsequencingdatawithAnders,S.,andHuber,W.(2010).Differentialexpression ysisforsequencecountdata.GenomeBiol.Anders,S.andHuber,W.(2012).DifferentialexpressionofRNA-Seqdataatthegenelevel-theDESeqKanehisa,M.,M.Araki,etal.(2008).KEGGforlinkinggenomestolifeandtheenvironment.NucleicacidsKim,D.,G.Pertea,etal.(2012).TopHat2:Parallelmapoftranscriptomestodetectindels,genefusions,andmore.(TopHat2)Langmead,B.,Trapnell,C.,Pop,M.&Salzberg,S.L.(2009).Ultrafastandmemory-efficientalignmentofshortDNAsequencestothehumangenome.GenomeBiol.(Bowtie)Langmead,B.andS.L.Salzberg(2012).Fastgapped-readalignmentwithBowtie2.Naturemethods.(Bowtie2)Mao,X.,Cai,T.,Olyarchuk,J.G.,Wei,L.(1995).AutomatedgenomeannotationandpathwayidentificationusingtheKEGGOrthology(KO)asacontrolledvocabulary.Bioinformatics.(KOBAS)Marioni,J.C.,C.E.Mason,etal.(2008).R