诺禾致源真核无参转录组生物信息分析结题报告2013年8月资料下载.pdf

1、北京诺禾致源生物信息科技有限公司二、生物信息分析流程对于无参考基因组的转录组分析，可先将测序所得的序列拼接成转录本，以转录本为参考序列，进行后续分析。信息分析流程图如下：4/38F:/结题报告+老销售培训/结题报告模板修改/真核无参转录组_Report.html北京诺禾致源生物信息科技有限公司三、结果展示及说明1原始序列数据高通量测序（如Illumina HiSeqTM2000/MiseqTM）得到的原始图像数据文件经CASAVA碱基识别（Base Calling）分析转化为原始测序序列（Sequenced Reads），我们称之为 Raw Data或Raw Reads，结果以 FASTQ（简

2、称为fq）文件格式存储，其中包含测序序列（reads）的序列信息以及其对应的测序质量信息。FASTQ格式文件中每个read由四行描述，如下：H W I-S T 1 2 7 6:7 1:C 1 1 6 2 A C X X:1:1 1 0 1:1 2 0 8:2 4 5 8 1:N:0:C G A T G T N A A G A A C A C G T T C G G T C A C C T C A G C A C A C T T G T G A A T G T C A T G G G A T C C A T+#5 5?B B B B B?B A D E E F F C F F H H F F C

3、 F F H H H H H H H F A E 0 E C F F D/A E H H其中第一行以“”开头，随后为Illumina 测序标识别符（Sequence Identifiers）和描述文字（选择性部分）；第二行是碱基序列；第三行以“+”开头，随后为Illumina 测序标识别符（选择性部分）；第四行是对应碱基的测序质量，该行中每个字符对应的 ASCII 值减去 33，即为对应第二行碱基的测序质量值。Illumina测序标识符详细信息如下：HWI-ST1276Instrument unique identifier of the sequencer71run number Run n

4、umber on instrumentC1162ACXXFlowCell ID ID of flowcell1LaneNumber positive integer1101TileNumber positive integer1208X x coordinate of the spot.Integer which can be negative2458Y y coordinate of the spot.Integer which can be negative1ReadNumber-1 for single reads;1 or 2 for paired endsNwhether it is

5、 filtered-NB：Y if the read is filtered out,not in the delivered fastq file,N otherwise0control number-0 when none of the control bits are on,otherwise it is an even numberCGATGThaoIllumina index sequences5/38F:/结题报告+老销售培训/结题报告模板修改/真核无参转录组_Report.html北京诺禾致源生物信息科技有限公司2测序数据质量评估2.1测序错误率分布检查如果测序错误率用e表示，I

6、llunima HiSeqTM2000/MiseqTM的碱基质量值用Qphred表示，则有：Qphred=-10log10（e）。Illunima Casava 1.8版本碱基识别与Phred分值之间的简明对应关系见下表：Phred分值不正确的碱基识别碱基正确识别率Q-sorce101/1090%Q10201/10099%Q20301/100099.9%Q30401/1000099.99%Q40对于RNA-seq技术，测序错误率分布具有两个特点，具体见图1：（1）测序错误率会随着测序序列（Sequenced Reads）的长度的增加而升高，这是由于测序过程中化学试剂的消耗而导致的，并且为ill

7、umina高通量测序平台都具有的特征（Erlich and Mitra,2008;Jiang et al.）。（2）前6个碱基的位置也会发生较高的测序错误率，而这个长度也正好等于在RNA-seq建库过程中反转录所需要的随机引物的长度。所以推测这部分碱基的测序错误率较高的原因为随机引物和RNA模版的不完全结合所致（Jiang et al.）。一般情况下，单个碱基位置的测序错误率应该低于1%。图1测序错误率分布图横坐标为reads的碱基位置，纵坐标为单碱基错误率前100bp 为双端测序序列的第一端测序 Reads的错误率分布情况，后100bp为另一端测序reads 的错误率分布情况。6/38F:/

8、结题报告+老销售培训/结题报告模板修改/真核无参转录组_Report.html北京诺禾致源生物信息科技有限公司2.2A/T/G/C 含量分布检查GC含量分布检查用于检测有无AT、GC分离现象。对于RNA-seq来说，因随机性打断及G/C和A/T含量分别相等的原则，理论上GC及AT含量每个测序循环上应分别近似相等（若为链特异性建库，可能会出现AT分离和/或GC分离），且整个测序过程基本稳定不变，呈水平线。但在现有的高通量测序技术中，反转录成cDNA 时所用的6bp的随机引物会引起前几个位置的核苷酸组成存在一定的偏好性，这种波动属于正常情况。如图 2 所示：图2GC含量分布图横坐标为reads的碱

9、基位置，纵坐标为单碱基所占的比例；不同颜色代表不同的碱基类型前100bp为双端测序序列的第一端测序Reads的GC分布情况，后100bp为另一端测序reads的GC分布情况。7/38F:/结题报告+老销售培训/结题报告模板修改/真核无参转录组_Report.html北京诺禾致源生物信息科技有限公司2.3测序数据过滤测序得到的原始测序序列（Sequenced Reads）或者 raw reads，里面含有带接头的、低质量的reads，如图3所示。为了保证信息分析质量，必须对raw reads过滤，得到clean reads，后续分析都基于 clean reads。数据处理的步骤如下：（1）去除带

10、接头（adapter）的reads；（2）去除N（N表示无法确定碱基信息）的比例大于10%的reads；（3）去除低质量reads（质量值sQ c o m p 0 _ c 0 _ s e q 1 l e n=3 0 5 p a t h=3 0 6:0-2 2 4 6 8:2 2 5-2 6 4 5 3 1:2 6 5-3 0 4 C A G G A C A C A A C A T A T A T C C G A G T T G G T G T G G C T A T C G A G C A G G G C T C G G A C A C A C T G A G GA T G C A T G C

11、 T C T A T G C T G G A G C C T C A T G G A T C G G G A C C G C C C A G C C T C A T T A C A A G G G T T TT A C T C T G T C G G A A A T G T G G C T T T G A T T T A G A G A A G G C A A C A A C T C C A G G C C A G G A C T G T AG G A G A G A G C A C A T C T A C A A C C A C T G C A G T C C C G G A C T G T

12、 A G G A G A G A C C A C A C C T A C A A CC A C T G C A G G C C A G G G C T G T A G G A G G G A C C A C A C C T A C A A C C A C T G C A G G C C A G G G C T G T AG G A G G c o m p 1 _ c 0 _ s e q 1 l e n=6 2 3 p a t h=1 5 6 5:0-6 2 2 A A A A A A A A G G G A A A C G T G A A A T C A A C G A C T A T T A

13、 T A A C C T T A A A T C A A A C T T A T C A C C A GG T G A A T T A A G C C C A T C T C T G T G G G C C A G A G C A T G T G T A T A A T T A C T T A A A T A C C A A T G T AG T G G G T T T T T A A C A A A C A T G A C A A C C T G T A G G A A A T G A A T G G T A G A T T A T G A A A T T A C T CC T C A T

14、T C A A C A T A C A T T T C T C A A C A G A A G C T A C A T T T G A T T A T G T G T T T G C T A G T T A C A TA T T A T G A C G C T T T G T A T T T T G A C A T T A A A G G G G T T T A A T A A C T T A C A G C A A G A C A G A T G TT A C A T A C C C A A T T A T A G A C C G G T T A C A A T T A C A G C T

15、G T T G G G T A G G A T T T G T G T T C C T T TC A A G A C A G T G A C A A C C T T C A C A G T T G G T A G T A T T G A G A G T A C A C A C T T G C A C A T A C A T G CT C T T C C C A T C T A A G C A C A T A A A T A A A C C A A A C A A A A A A A A A C T G G T T G G C A A G G T G T G T GT G T A T A C A

16、 C A T G A A A G C A G A T T A T A T C T T A A C C T T A G A T C A A A C T T G A A T T A C T T G A T C CA A G A A G A G G A A C T A C A G T C C A G C A G T A A A G G G A C A T A C A T T T G T T T A C A T T T A C T T T G G GA C C T G C T C C A T G T T T A T A C A A T T其中大于号后紧跟转录本的id号，len=后面为转录本的长度，即该

17、转录本的碱基数，path为从 de Bruijn Graph subComponent中经历的路径。其后为该转录的碱基序列。每个转录本的id号构成都为comp_c_seq,其中comp为拼接过程形成的de Bruijn Graph Component,c为10/38F:/结题报告+老销售培训/结题报告模板修改/真核无参转录组_Report.htmlsubcomponet，可以看作为广泛意义上的基因，seq 代表转录本。详细解释见Trinity官方网站 拼接转录本长度分布用 Trinity 将测序序列拼接成一个转录组，以此作为后续分析的参考序列。取每条基因中最长的转录本作为Unigene。对转录

18、本及Unigene的长度进行统计，结果见表2，表3和图4，图5。表2拼接长度频数分布情况一览表Transcript length interval200-500bp500-1kbp1k-2kbp2kbpTotalNumber of transcripts100027443455729343241244906Number of unigenes7364317186113397632109800表3拼接长度分布情况一览表Min LengthMean LengthMedian LengthMax LengthN50N90Transcripts2011112712167961886457Unigene

19、s201661340167961166258N50/N90的定义为：按照长度将拼接转录本从大到小排序，累加转录本的长度，到不小于总长50%/90%的拼接转录本的长度就是N50/N90。图4 拼接转录本长度分布图横坐标为拼接转录本的长度区间，纵坐标为对应长度的频数。图4 拼接得到的Unigene长度分布图横坐标为拼接的到的Unigene的长度区间，纵坐标为对应长度的频数。12/38F:/结题报告+老销售培训/结题报告模板修改/真核无参转录组_Report.html北京诺禾致源生物信息科技有限公司4基因功能注释4.1 基因功能注释基因功能注释所用到的数据库：Nr（NCBI non-redundan

20、t protein sequences）是NCBI官方的蛋白序列数据库，它包括了GenBank基因的蛋白编码序列，PDB（Protein DataBank）蛋白数据库、SwissProt蛋白序列及来自PIR（Protein Information Resource）和PDF（Protein Research Foundation）等数据库的蛋白序列。Nt（NCBI nucleotide sequences）是NCBI官方的核酸序列数据库，包括了GenBank，EMBL和DDBJ（但不包括EST,STS,GSS,WGS,TSA,PAT,HTG序列）的核酸序列。Pfam（Protein famil

21、y）是最全面的蛋白结构域注释的分类系统。蛋白质是由一个个结构域组成的，而每个特定结构域的蛋白序列具有一定保守性。PFAM将蛋白质的结构域分为不同的蛋白家族，通过蛋白序列的比对建立了每个家族的氨基酸序列的HMM统计模型。PFAM家族按注释结果可靠性分为两大类：手工注释的可靠性高的Pfam-A家族和程序自动产生Pfam-B家族。我们通过HMMER3程序，搜索已建好的蛋白结构域的HMM模型，对unigene进行了蛋白家族的注释。详见http:/pfam.sanger.ac.uk/。KOG/COG:COG是Clusters of Orthologous Groups of proteins的简称，KO

22、G为euKaryotic Ortholog Groups。这两个注释系统都是NCBI的基于基因直系同源关系，其中COG针对原核生物，KOG针对真核生物。COG/KOG结合进化关系将来自不同物种的同源基因分为不同的Ortholog簇，目前COG有4873个分类，KOG有4852个分类。来自同一ortholog的基因具有相同的功能，这样就可以将功能注释直接继承给同一COG/KOG簇的其他成员。/www.ncbi.nlm.nih.gov/COG/。Swiss-Prot（A manually annotated and reviewed protein sequence database）搜集了经过有

23、经验的生物学家整理及研究的蛋白序列。/www.ebi.ac.uk/uniprot/。KEGG是Kyoto Encyclopedia of Genes and Genomes的简称，是系统分析基因产物和化合物在细胞中的代谢途径以及这些基因产物的功能的数据库。它整合了基因组、化学分子和生化系统等方面的数据,包括代谢通路（KEGG PATHWAY）、药物（KEGG DRUG）、疾病（KEGG DISEASE）、功能模型（KEGG MODULE）、基因序列（KEGG GENES）及基因组（KEGG GENOME）等等。KO（KEGGORTHOLOG）系统将各个KEGG注释系统联系在一起，KEGG已建立

24、了一套完整KO注释的系统，可完成新测序物种的基因组或转录组的功能注释。详见 http:/www.genome.jp/kegg/。GO（Gene Ontology）是一套国际标准化的基因功能描述的分类系统。GO分为三大类ontology：生物过程（Biological Process）、分子功能（Molecular Function）和细胞组分（Cellular Component），分别用来描述基因编码的产物所参与的生物过程、所具有的分子功能及所处的细胞环境。GO的基本单元是term，每个term有一个唯一的标示符（由“GO:”加上7个数字组成，例如GO:0072669）；每类ontology

25、的term通过它们之间的联系（is_a,part_of,regulate）构成一个有向无环的拓扑结构。/www.geneontology.org/。各数据库及功能注释所用到的软件及方法：与NR、NT、SwissProt、KOG、KEGG GENES序列数据库的比对：NCBI blast 2.2.27+;PFAM蛋白结构域预测：HMMER 3.0 package，hmmscan;GO功能注释：基于NR和Pfam两部分的蛋白注释结果:软件为Blast2GO v2.5（Gtz et al.,2008）和自写脚本;KEGG相关注释：KAAS，KEGG Automatic Annotation Serv

26、er。注释结果统计见表4。表4 基因注释成功率统计Number of UnigenesPercentage（%）Annotated in NR5226647.6Annotated in NT2461222.41Annotated in KO2419922.03Annotated in SwissProt4271438.9Annotated in PFAM4309839.25Annotated in GO5447949.61Annotated in KOG2790725.41Annotated in all Databases93018.47Annotated in at least one D

27、atabase6159156.0913/38F:/结题报告+老销售培训/结题报告模板修改/真核无参转录组_Report.htmlTotal Unigenes109800100Annotated in NR:NR注释成功的Unigene数目及其占总Unigene数的比例Annotated in NT:NT比对成功的Unigene数目及其占总Unigene数目的比例Annotated in KO:KO注释成功的Unigene数目及其占总Unigene数的比例Annotated in Swissprot:Swissprot注释成功的Unigene数目及其占总Unigene数的比例Annotated

28、in Pfam:Pfam注释成功的Unigene数目及其占总Unigene数的比例Annotated in GO:GO注释成功的Unigene数目及其占总Unigene数的比例Annotated in KOG:KOG注释成功的Unigene数目及其占总Unigene数的比例Annotated in all Databases:在以上7个数据库中都注释成功的Unigene数目及其占总Unigene数的比例Annotated in at least one Database:在以上7个数据库中至少1个数据库注释成功的Unigene数目及其占总Unigene数的比例Total Unigenes:总的Unigene条数,占总Unigene比例为100%注释结果示意部分见表5。表5 基因注释结果（此处只展示 NR的部分比对结果，其它数据库比对注释结果详见结果文件夹中表格）Gene IDGene LengthNR GINR IDNR ScoreNR Evaluecomp943620_c022547087029NP_998528.13352.20751e-35comp116203_c02571255543897XP_002513011.125640comp453168_c0