无参考基因的转录组分析

无参考基因的转录组分析一、实验流程提取样品总RNA后，用带有Oligo（dT）的磁珠富集真核生物mRNA（若为原核生物，则用试剂盒去除rRNA后进入下一步）。加入fragmentationbuffer将mRNA打断成短片段，以mRNA为模板，用六碱基随机引物（randomhexamers）合成第一条cDNA链，然后加入缓冲液、dNTPs、RNaseH和DNApolymeraseI合成第二条cDNA链，在经过QiaQuickPCR试剂盒纯化并加EB缓冲液洗脱之后做末端修复、加poly（A）并连接测序接头，然后用琼脂糖凝胶电泳进行片段大小选择，最后进行PCR扩增，建好的测序文库用IlluminaHiSeq™2000进行测序。二、信息分析流程1、产量统计原始序列数据测序得到的原始图像数据经basecalling转化为序列数据，我们称之为rawdata或rawreads，结果以fastq文件格式存储，fastq文件为用户得到的最原始文件，里面存储reads的序列以及reads的测序质量。在fastq格式文件中每个read由四行描述：\@FC61FL8AAXX:1:17:1012:19200#GCCAAT/1CCACTGTCATGTGAACATCACAGAGACATTTCTTGA+bbbbbbbbbbbbbbbbbbbbbbbbbaaaaaaaaa_\\每个序列共有4行，第1行和第3行是序列名称（有的fq文件为了节省存储空间会省略第三行“＋”后面的序列名称），由测序仪产生；第2行是序列；第4行是序列的测序质量，每个字符对应第2行每个碱基，第四行每个字符对应的ASCII值减去64，即为该碱基的测序质量值，比如c对应的ASCII值为99，那么其对应的碱基质量值是35。从IlluminaGAPipelinev1.3开始（目前为v1.6），碱基质量值范围为2到35。表1为测序错误率与测序质量值简明对应关系。具体地，如果测序错误率用E表示，碱基质量值用sQ表示，则有下列关系：sQ=-10lgE表1测序错误率与测序质量值简明对应关系测序错误率测序质量值对应字符5%13M1%20T0.1%30^去除杂质数据某些原始序列带有adaptor序列，或含有少量低质量序列。我们首先经过一系列数据处理以去除杂质数据，得到Cleanreads。数据处理的步骤：1.去除含adaptor的reads2.去除N的比例大于10%的reads3.去除低质量reads（质量值Q≤5的碱基数占整个read的50％以上）4.获得CleanreadsCleanReads数据原始序列数据经过去除杂质后得到的数据。产量统计和后续信息分析分析都基于CleanReads。测序产量统计表格示例SamplesTotalReadsTotalNucleotides(nt)Q20percentageNpercentageGCpercentage*Sample_A1,634,670122,600,25089.47%0.00%48.50%*TotalNucleotides=TotalReads1xRead1size+TotalReads2xRead2size;TotalReadsandTotalNucleotidesareactuallycleanreadsandcleannucleotides;Q20percentageisproportionofnucleotideswithqualityvaluelargerthan20;Npercentageisproportionofunknownnucleotidesincleanreads;GCpercentageisproportionofguanidineandcytosinenucleotidesamongtotalnucleotides.2、组装结果我们使用短reads组装软件SOAPdenovo[5]做转录组从头组装。SOAPdenovo首先将具有一定长度overlap的reads连成更长的片段，这些通过readsoverlap关系得到的不含N的组装片段我们称之称为Contig。然后，我们将reads比对回Contig，通过paired-endreads能确定来自同一转录本的不同Contig以及这些Contig之间的距离，SOAPdenovo将这些Contig连在一起，中间未知序列用N表示，这样就得到Scaffold。进一步利用paired-endreads对Scaffold做补洞处理，最后得到含N最少，两端不能再延长的序列，我们称之为Unigene。如果同一物种做了多个样品测序，则不同样品组装得到的Unigene可通过序列聚类软件做进一步序列拼接和去冗余处理，得到尽可能长的非冗余Unigene。最后，将Unigene序列与蛋白数据库nr、Swiss-Prot、KEGG和COG做blastx比对（evalue0.00001），取比对结果最好的蛋白确定Unigene的序列方向。如果不同库之间的比对结果有矛盾，则按nr、Swiss-Prot、KEGG和COG的优先级确定Unigene的序列方向，跟以上四个库皆比不上的Unigene我们用软件ESTScan[3]预测其编码区并确定序列的方向。对于能确定序列方向的Unigene我们给出其从5'到3'方向的序列，对于无法确定序列方向的Unigene我们给出组装软件得到的序列。组装出来的序列长度是组装质量的一个评估标准。我们会对组装出来的Contig、Scaffold、Unigene做一个长度分布统计。如下图所示，给出的bar图统计Contig的长度分布。横坐标是组装出来的Contig的长度，纵坐标是对应长度的Contig的数目。组装成功的Contig结果在文件夹1.Contig，Scaffold相关结果在文件夹2.Scaffold，Unigene相关结果在文件夹3.Unigene。文件的详细意义可见各个文件夹下面对应的readme。注：文件夹下面svg图可能需要安装svg插件才能打开3、Unigene功能注释功能注释信息给出Unigene的蛋白功能注释、COG功能注释。首先，通过blastx将Unigene序列比对到蛋白数据库nr、Swiss-Prott、KEGG和COG（evalue0.00001），得到跟给定Unigene具有最高序列相似性的蛋白，从而得到该Unigene的蛋白功能注释信息。COG是对基因产物进行直系同源分类的数据库，每个COG蛋白都被假定来自祖先蛋白，COG数据库是基于细菌、藻类、真核生物具有完整基因组的编码蛋白、系统进化关系进行构建的，我们将Unigene和COG数据库进行比对，预测Unigene可能的功能并对其做功能分类统计，从宏观上认识该物种的基因功能分布特征。Unigene的COG功能注释结果样式示例4、Unigene的GO分类根据nr注释信息我们能得到GO功能注释。GeneOntology（简称GO）是一个国际标准化的基因功能分类体系，提供了一套动态更新的标准词汇表（controlledvocabulary）来全面描述生物体中基因和基因产物的属性。GO总共有三个ontology，分别描述基因的分子功能（molecularfunction）、所处的细胞位置（cellularcomponent）、参与的生物过程（biologicalprocess）。我们根据nr注释信息，使用Blast2GO软件（Conesa,Gotzetal.2005）得到Unigene的GO注释信息。Blast2GO已被其它文献引用超过150次，是同行广泛认可的GO注释软件。得到每个Unigene的GO注释后，我们用WEGO软件（Ye,Fangetal.2006）对所有Unigene做GO功能分类统计，从宏观上认识该物种的基因功能分布特征。基因注释到GO条目结果文件示例GO条目与All-Unigene对应结果文件示例5、Unigene代谢通路分析KEGG是系统分析基因产物在细胞中的代谢途径以及这些基因产物的功能的数据库，利用KEGG可以进一步研究基因在生物学上的复杂行为。根据KEGG注释信息我们能进一步得到Unigene的Pathway注释。注释到代谢通路结果文件示例6、预测编码蛋白框（CDS）首先，我们按nr、Swiss-Prot、KEGG和COG的优先级顺序将Unigene序列与以上蛋白库做blastx比对（evalue0.00001），如果某个Unigene序列比对上高优先级数据库中的蛋白，则不进入下一轮比对，否则自动跟下一个库做比对，如此循环直到跟所有蛋白库比对完。我们取blast比对结果中rank最高的蛋白确定该Unigene的编码区序列，然后根据标准密码子表将编码区序列翻译成氨基酸序列，从而得到该Unigene编码区的核酸序列（序列方向5'-3'）和氨基酸序列。最后，跟以上蛋白库皆比对不上的Unigene我们用软件ESTScan（Iseli,Jongeneeletal.1999）预测其编码区，得到其编码区的核酸序列（序列方向5'-3'）和氨基酸序列。'7、Unigene表达差异分析差异表达分析找出在不同样本间存在差异表达的基因，并对差异表达基因做GO功能分析和KEGGPathway分析。差异表达基因的筛选参照AudicS.等人发表在GenomeResearch上的基于测序的差异基因检测方法[1]，我们开发了严格的算法筛选两样本间的差异表达基因。筛选差异表达基因所用的到的假设检验的零假设和备择假设如下：:0H某一个基因在两个样本中表达量相同:aH某一个基因在两个样本中表达量不同假设观测到基因A对应的reads数为x，已知在一个大文库中，每个基因的表达量只占所有基因表达量的一小部分，在这种情况下，p(x)的分布服从泊松分布[1]：!)(xexpx（为基因A的真实转录数）已知，样本一能比对到所有Unigene的总reads数为1N，样本二能比对到所有Unigene的总reads数为2N，基因A在样本一中对应的reads数为x，在样本二中对应的reads数为y，则基因A在两样本中表达量相等的概率可由以下公式计算[1]：yixip0)|(2（当5.0)|(0yixip时）或者))|(1(20yixip（当5.0)|(0yixip时）其中)1(1212)1(!!)!()()|(ixiNNyxixNNxip然后，我们对差异检验的p-value作多重假设检验校正，通过控制FDR（FalseDiscoveryRate）来决定pvalue的域值。假设挑选了R个差异表达基因，其中S个是真正有差异表达的基因，另外V个是其实没有差异表达的基因，为假阳性结果。希望错误比例Q＝V/R平均而言不能超过某个可以容忍的值（比如1％），则在统计时预先设定FDR不能超过0.01（Benjamini,Yekutieli.2001）。在得到差异检验的FDR值同时，我们也根据基因的表达量（RPKM值）计算该基因在不同样本间的差异表达倍数。FDR值越小，差异倍数越大，则表明表达差异越显著。在我们的分析中，差异表达基因定义为FDR≤0.001且倍数差异在2倍以上的基因。得到差异表达基因之后，我们对差异表达基因做GO功能分析和KEGGPathway分析。表达量注释Unigene表达量的计算使用RPKM法（ReadsPerkbperMillionreads）[6]，其计算公式为：3610/10NLCRPKM设RPKM(A)为UnigeneA的表达量，则C为唯一比对到UnigeneA的reads数，N为唯一比对到所有Unigene的总reads数，L为UnigeneA的碱基数。RPKM法能消除基因长度和测序量差异对计算基因表达的影响，计算得到的基因表达量可直接用于比较不同样品间的基因表达差异。8、差异Unigene的GO和Pathway分析GO功能分析GeneOntology（简称GO）是一个国际标准化的基因功能分类体系，提供了一套动态更新的标准词汇表来全面描述生物体中基因和基因产