1、1/38F:/结题报告+老销售培训/结题报告模板修改/真核无参转录组_Report.html真核无参转录组生物信息分析结题报告建库测序流程Total RNA样品检测文库构建上机测序生物信息分析流程结果展示及说明原始序列数据测序数据质量评估测序错误率分布检查A/T/G/C 含量分布检查测序数据过滤测序数据质量情况汇总转录本拼接转录本拼接拼接转录本长度分布基因功能注释基因功能注释GO分类KOG分类KEGG分类CDS 预测SNP分析SSR分析SSR分析SSR引物设计基因表达水平分析参考序列比对基因表达水平统计表基因RPKM 密度分布图RNA-seq整体质量评估表达水平的饱和曲线检查均一化分布样品间相
2、关性检查差异表达分析不同实验条件下，基因表达水平对比图基因差异表达分析差异基因筛选差异基因表达水平聚类分析GO 富集分析差异表达基因GO富集分析topGO有向无环图分析KEGG 富集分析差异基因KEGG富集分析差异基因KEGG富集散点图富集KEGG通路图备注参考文献2/38F:/结题报告+老销售培训/结题报告模板修改/真核无参转录组_Report.html北京诺禾致源生物信息科技有限公司一、建库测序流程从RNA样品到最终数据获得，样品检测、建库、测序每一个环节都会对数据质量和数量产生影响，而数据质量又会直接影响后续信息分析的结果。为了从源头上保证测序数据的准确性、可靠性，诺禾致源对样品检测、建
3、库、测序每一个生产步骤都严格把控，从根本上确保了高质量数据的产出。实验流程图如下：1 Total RNA样品检测诺禾致源对RNA样品的检测主要包括4种方法：(1)琼脂糖凝胶电泳分析RNA降解程度以及是否有污染(2)Nanodrop检测RNA的纯度（OD260/280比值）(3)Qubit对RNA浓度进行精确定量(4)Agilent 2100精确检测RNA的完整性2 文库构建及库检样品检测合格后，用带有Oligo（dT）的磁珠富集真核生物mRNA（若为原核生物，则通过试剂盒去除rRNA来富集mRNA）。随后加入fragmentation buffer将mRNA打断成短片段，以mRNA为模板，用六
4、碱基随机引物（random hexamers）合成一链cDNA，然后加入缓冲液、dNTPs、RNase H和DNA polymerase I合成二链cDNA，随后利用AMPure XP beads纯化双链cDNA。纯化的双链cDNA再进行末端修复、加A尾并连接测序接头，然后用AMPure XP beads进行片段大小选择，最后进行PCR富集得到最终的cDNA文库。文库构建完成后，先使用Qubit2.0进行初步定量，稀释文库至1ng/ul，随后使用Agilent 2100对文库的insert size进行检测，insert size符合预期后，使用Q-PCR方法对文库的有效浓度进行准确定量（文库
5、有效浓度 2nM），以保证文库质量。文库构建原理图如下：3/38F:/结题报告+老销售培训/结题报告模板修改/真核无参转录组_Report.html3 上机测序库检合格后，把不同文库按照有效浓度及目标下机数据量的需求pooling后进行Illumina HiSeq/MiSeq测序。北京诺禾致源生物信息科技有限公司二、生物信息分析流程对于无参考基因组的转录组分析，可先将测序所得的序列拼接成转录本，以转录本为参考序列，进行后续分析。信息分析流程图如下：4/38F:/结题报告+老销售培训/结题报告模板修改/真核无参转录组_Report.html北京诺禾致源生物信息科技有限公司三、结果展示及说明1原始
6、序列数据高通量测序(如Illumina HiSeqTM2000/MiseqTM)得到的原始图像数据文件经CASAVA碱基识别(Base Calling)分析转化为原始测序序列（Sequenced Reads），我们称之为 Raw Data或Raw Reads，结果以 FASTQ(简称为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
7、 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 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 测序标识别符(选择性部分)；第四行是对应
8、碱基的测序质量，该行中每个字符对应的 ASCII 值减去 33，即为对应第二行碱基的测序质量值。Illumina测序标识符详细信息如下：HWI-ST1276Instrument unique identifier of the sequencer71run number Run number on instrumentC1162ACXXFlowCell ID ID of flowcell1LaneNumber positive integer1101TileNumber positive integer1208X x coordinate of the spot.Integer which c
9、an be negative2458Y y coordinate of the spot.Integer which can be negative1ReadNumber-1 for single reads;1 or 2 for paired endsNwhether it is 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
10、 even numberCGATGThaoIllumina index sequences5/38F:/结题报告+老销售培训/结题报告模板修改/真核无参转录组_Report.html北京诺禾致源生物信息科技有限公司2测序数据质量评估2.1测序错误率分布检查如果测序错误率用e表示，Illunima HiSeqTM2000/MiseqTM的碱基质量值用Qphred表示，则有：Qphred=-10log10(e)。Illunima Casava 1.8版本碱基识别与Phred分值之间的简明对应关系见下表：Phred分值不正确的碱基识别碱基正确识别率Q-sorce101/1090%Q10201/100
11、99%Q20301/100099.9%Q30401/1000099.99%Q40对于RNA-seq技术，测序错误率分布具有两个特点，具体见图1：(1)测序错误率会随着测序序列(Sequenced Reads)的长度的增加而升高，这是由于测序过程中化学试剂的消耗而导致的，并且为illumina高通量测序平台都具有的特征(Erlich and Mitra,2008;Jiang et al.)。(2)前6个碱基的位置也会发生较高的测序错误率，而这个长度也正好等于在RNA-seq建库过程中反转录所需要的随机引物的长度。所以推测这部分碱基的测序错误率较高的原因为随机引物和RNA模版的不完全结合所致(Ji
12、ang et al.)。一般情况下，单个碱基位置的测序错误率应该低于1%。图1测序错误率分布图横坐标为reads的碱基位置，纵坐标为单碱基错误率前100bp 为双端测序序列的第一端测序 Reads的错误率分布情况，后100bp为另一端测序reads 的错误率分布情况。6/38F:/结题报告+老销售培训/结题报告模板修改/真核无参转录组_Report.html北京诺禾致源生物信息科技有限公司2.2A/T/G/C 含量分布检查GC含量分布检查用于检测有无AT、GC分离现象。对于RNA-seq来说，因随机性打断及G/C和A/T含量分别相等的原则，理论上GC及AT含量每个测序循环上应分别近似相等(若为
13、链特异性建库，可能会出现AT分离和/或GC分离)，且整个测序过程基本稳定不变，呈水平线。但在现有的高通量测序技术中，反转录成cDNA 时所用的6bp的随机引物会引起前几个位置的核苷酸组成存在一定的偏好性，这种波动属于正常情况。如图 2 所示：图2GC含量分布图横坐标为reads的碱基位置，纵坐标为单碱基所占的比例；不同颜色代表不同的碱基类型前100bp为双端测序序列的第一端测序Reads的GC分布情况，后100bp为另一端测序reads的GC分布情况。7/38F:/结题报告+老销售培训/结题报告模板修改/真核无参转录组_Report.html北京诺禾致源生物信息科技有限公司2.3测序数据过滤测
14、序得到的原始测序序列（Sequenced Reads）或者 raw reads，里面含有带接头的、低质量的reads，如图3所示。为了保证信息分析质量，必须对raw reads过滤，得到clean reads，后续分析都基于 clean reads。数据处理的步骤如下：(1)去除带接头(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
15、-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 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
16、 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 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
17、 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 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
18、 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 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
19、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 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
20、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 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
21、 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=后面为转录本的长度，即该转录本的碱基数，path为从 de Bruijn Graph subComponent中经历的路径。其后为该转录的碱基序列。每个转录本的id号构成都为comp_c_seq,其中comp为拼接过程形成的de Bruijn Graph Component,c为10/38F:/结题报告+老销售培训/结题报告模板修改/真核无
22、参转录组_Report.htmlsubcomponet，可以看作为广泛意义上的基因，seq 代表转录本。详细解释见Trinity官方网站 拼接转录本长度分布用 Trinity 将测序序列拼接成一个转录组，以此作为后续分析的参考序列。取每条基因中最长的转录本作为Unigene。对转录本及Unigene的长度进行统计，结果见表2，表3和图4，图5。表2拼接长度频数分布情况一览表Transcript length interval200-500bp500-1kbp1k-2kbp2kbpTotalNumber of transcripts100027443455729343241244906Numbe
23、r of unigenes7364317186113397632109800表3拼接长度分布情况一览表Min LengthMean LengthMedian LengthMax LengthN50N90Transcripts2011112712167961886457Unigenes201661340167961166258N50/N90的定义为：按照长度将拼接转录本从大到小排序，累加转录本的长度，到不小于总长50%/90%的拼接转录本的长度就是N50/N90。图4 拼接转录本长度分布图横坐标为拼接转录本的长度区间，纵坐标为对应长度的频数。图4 拼接得到的Unigene长度分布图横坐标为拼接的
24、到的Unigene的长度区间，纵坐标为对应长度的频数。12/38F:/结题报告+老销售培训/结题报告模板修改/真核无参转录组_Report.html北京诺禾致源生物信息科技有限公司4基因功能注释4.1 基因功能注释基因功能注释所用到的数据库：Nr(NCBI non-redundant protein sequences)是NCBI官方的蛋白序列数据库，它包括了GenBank基因的蛋白编码序列，PDB(Protein DataBank)蛋白数据库、SwissProt蛋白序列及来自PIR（Protein Information Resource）和PDF（Protein Research Foun
25、dation）等数据库的蛋白序列。Nt(NCBI nucleotide sequences)是NCBI官方的核酸序列数据库，包括了GenBank，EMBL和DDBJ（但不包括EST,STS,GSS,WGS,TSA,PAT,HTG序列）的核酸序列。Pfam(Protein family)是最全面的蛋白结构域注释的分类系统。蛋白质是由一个个结构域组成的，而每个特定结构域的蛋白序列具有一定保守性。PFAM将蛋白质的结构域分为不同的蛋白家族，通过蛋白序列的比对建立了每个家族的氨基酸序列的HMM统计模型。PFAM家族按注释结果可靠性分为两大类：手工注释的可靠性高的Pfam-A家族和程序自动产生Pfam-
26、B家族。我们通过HMMER3程序，搜索已建好的蛋白结构域的HMM模型，对unigene进行了蛋白家族的注释。详见http:/pfam.sanger.ac.uk/。KOG/COG:COG是Clusters of Orthologous Groups of proteins的简称，KOG为euKaryotic Ortholog Groups。这两个注释系统都是NCBI的基于基因直系同源关系，其中COG针对原核生物，KOG针对真核生物。COG/KOG结合进化关系将来自不同物种的同源基因分为不同的Ortholog簇，目前COG有4873个分类，KOG有4852个分类。来自同一ortholog的基因具有
27、相同的功能，这样就可以将功能注释直接继承给同一COG/KOG簇的其他成员。详见http:/www.ncbi.nlm.nih.gov/COG/。Swiss-Prot（A manually annotated and reviewed protein sequence database）搜集了经过有经验的生物学家整理及研究的蛋白序列。详见http:/www.ebi.ac.uk/uniprot/。KEGG是Kyoto Encyclopedia of Genes and Genomes的简称，是系统分析基因产物和化合物在细胞中的代谢途径以及这些基因产物的功能的数据库。它整合了基因组、化学分子和生化系统
28、等方面的数据,包括代谢通路（KEGG PATHWAY）、药物（KEGG DRUG）、疾病（KEGG DISEASE）、功能模型（KEGG MODULE）、基因序列（KEGG GENES）及基因组（KEGG GENOME）等等。KO（KEGGORTHOLOG）系统将各个KEGG注释系统联系在一起，KEGG已建立了一套完整KO注释的系统，可完成新测序物种的基因组或转录组的功能注释。详见 http:/www.genome.jp/kegg/。GO(Gene Ontology)是一套国际标准化的基因功能描述的分类系统。GO分为三大类ontology：生物过程（Biological Process）、分子
29、功能（Molecular Function）和细胞组分(Cellular Component)，分别用来描述基因编码的产物所参与的生物过程、所具有的分子功能及所处的细胞环境。GO的基本单元是term，每个term有一个唯一的标示符（由“GO:”加上7个数字组成，例如GO:0072669）；每类ontology的term通过它们之间的联系（is_a,part_of,regulate）构成一个有向无环的拓扑结构。详见http:/www.geneontology.org/。各数据库及功能注释所用到的软件及方法：与NR、NT、SwissProt、KOG、KEGG GENES序列数据库的比对：NCBI
30、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 Server。注释结果统计见表4。表4 基因注释成功率统计Number of UnigenesPercentage(%)Annotated in NR5226647.6Annotated in NT2461222.41Annotated in KO2419922.03Annotate
31、d in SwissProt4271438.9Annotated in PFAM4309839.25Annotated in GO5447949.61Annotated in KOG2790725.41Annotated in all Databases93018.47Annotated in at least one Database6159156.0913/38F:/结题报告+老销售培训/结题报告模板修改/真核无参转录组_Report.htmlTotal Unigenes109800100Annotated in NR:NR注释成功的Unigene数目及其占总Unigene数的比例Anno
32、tated in NT:NT比对成功的Unigene数目及其占总Unigene数目的比例Annotated in KO:KO注释成功的Unigene数目及其占总Unigene数的比例Annotated in Swissprot:Swissprot注释成功的Unigene数目及其占总Unigene数的比例Annotated in Pfam:Pfam注释成功的Unigene数目及其占总Unigene数的比例Annotated in GO:GO注释成功的Unigene数目及其占总Unigene数的比例Annotated in KOG:KOG注释成功的Unigene数目及其占总Unigene数的比例A
33、nnotated 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
