全长转录组有参报账用结题报告

合同关键指标

2+3全长有参转录组服务合完成1个样品的全长转录组，样品产出不少于8GbCleanData。完成可变剪接分析。完成长链非编码RNA预测和长链非编码RNA靶预测。完成转录本功能注释分析。实验流全长转录组实验流程包括样品检测、文库构建和上机。实验流程见下图全长转录组实验流程样品检高质量的RNA是整个项目成功的基础，为保证数据准确性，我们使用以下方法对样品进行检测，检测结果达到要求后方可进行库Nanodrop检测RNA的纯度（OD260/280）、浓度、核酸吸收峰是否正常Agilent2100精确检测RNA的完整性，检测指标包括:RIN值、28S/18S、图谱基线有无上抬、5S峰电泳检测RNA样品是否有组DNA的污文库构样品检测合格后，进行文库构建，主要流程如下使用SMARTer™PCRcDNASynthesisKit合成mRNA的cDNABluePippin筛选全cDNA片段，构建不同大小cDNA文库:。(3)再次PCR扩增放大筛选的全长cDNA。对全长cDNA进行末端修复连接SMRT哑铃型接头进行核酸外切酶消化使用BluePippin进行二次筛选，获得文库文库质文库构建完成后，对文库质量进行检测，检测结果达到要求后方可进行上机，检测方法如下使用Qubit2.0进行准确定量使用Agilent2100对文库大小进行检测，文库大小符合预期后才可进行上机上库检合格后，按照目标下机数据量，使用PacBio仪器进行全长转录组生物信息学转录组研究是理解生命过程必不可少的工具之一，然而基于第二代高通量平台的R-e2.0技术往往不能准确得到或组装出完整转录本，无法识别Isoform、同源、超、等位表达的转录本，使人们难以理解这一生命活动更次的含义。基于acioRT单分子实时技术的全长转录组无须打断A片段，使用AE技术反转录得到的全长cDA。该平台的超长读取（edian10kb）包含了单条完整转录本序列信息，后期分析无需组装，所测即所得[1][2][3]。获取全长转录组的分析过程主要包括3个阶段[4]，全长序列识别、isoform水平聚类得到一致性序列和一致性序列polishin，详细步骤如下：从原始下机序列中提取ROI（ReadsOfInsert）序列，过滤序列中的cDNA引物和polyA，并根据序列中是否存在3'引物、5'引物和PolyA（可选）将序列分成全长序列与非全长序列、嵌合序列与非嵌合序列。使用ICE（Iterativeisoform-clustering）算法将来自同一isoform的全长序列进行聚类，将具有相似序列的全长序列聚成一簇（cluster），每个cluster得到一条一致性序列使用Quiver算法聚类非全长序列，对得到的一致性序列进行校正（polishing），筛选高质量的序列进行后续分析。考虑到cNA建库的局限性[4]，我们筛选得到的高质量序列由于建库时5'端序列的缺失可能是非全长序列，因此，仅5'末端外显子存在差异，其余外显子均一致的序列进行合并，取其中最长的序列作为最终的转录本序列最终得到的转录本序列可直接用于后续的Isoform、同源、、等位、SS、可变剪接、lncNA等分析。引导人们更次的理解位于中心法则中心地位的这一生命活动，另外还可用于对所在组的注释升级，完善组数据库。转录组生物信息分析流程见下图全长转录组生物信息分析流程2.2.1数基于单分子（SequencingBySynthesis，SMRT）技术，Pacbio高通量平台对cDNA文库进序，产出大量的数据，称为原始数据（RawData）。结果文件列(1)pacbio文件命名规则pacbio产出文件或序列命名方式示例如下pacbio命名示"m"=TimeofRunStart(yymmdd_hhmmss)上机时InstrumentSerialNumber仪器序列编SMRTCellBarcodeSMRTcellbarcode编SetNumber(a.k.a."LookNumber".Deprecatedfield,usedinearlierversionofRS)6PartNumber(usually"p0","X0"whenusingexpiredreagents)ZMWholenumberSubreadRegion(start_stopusingpolymerasereadcoordinates)†subread在原始序列中的起始和终止(2)主要文件介绍主要生成*.bas.h5,*.1.bax.h5,*.2.bax.h5,和*.3.bax.h5文件。bax.h5文件包括每个循环中碱基识别信息。bas.h5文件bax.h5文件的描述（pointer）信息。bas.h5文件说明请考 /software/instrument/2.0.0/bas.h5+Reference+Guide.pdfetadat.xl文件主要包括数据的基本信息，包括所使用的酶和化学试剂信息、样品名称和其他实验信息。详细说明请参考:https: /software/instrument/2.0.0/Metadata+Output+Guide.pdf。sts.xml文件主要包括每轮（movie）中的统计信息说明，详细说明请考 /software/instrument/1.3.1/Statistics+Output+Guide.pdfsubread序列文件：原始序列（Polymerasereads）从接头处打断并过滤掉接头序列后得到subread，该序列FASTQFASTA格式保存，其中FASTQ主要包括序列和质量信息，FASTQ格式文件示意图如下:FASTQ格式文件示意注:FASTQ文件中通常每4行对应一个序列单元:第一行以@开头，后面接着序列标识（I）以及其它可选的描述信息;第二行为碱基序列，即es;第三行以开头，后面接着可选的描述信息;第四行为as每个碱基对应的质量打分编码，长度必须和es的序列长度相同，该行中每个字符对应的ASII值减去，即为对应第二行碱基的质量值。中每个ZMW(zero- waveguide)中得到的原始序列如下图所示序列介Polymerase聚合酶序列（polymeraseread）:DNA聚合酶以SMRTbell™环状模板链合成的核酸序列，可用于过程中每轮（run）的质控。olyerasereds经过滤后仅剩余高质量片段，包含接头序列和通过环状模板链合成的含多个序列的拷贝，如图6中'ra'所示。每个聚合酶序列（polymeraseread）可以分割成一个或多个子序列（Subread），subread是聚合酶以SMRTbell™一条模板链经过一轮（passes）合成的，不包括接头序列。每个subreads包含质量值和相关酶活参数。图6共有2个半的subread。其中fullpasses（完整的subread）数为2。numberoffull指原始序列中存在两端均含有SMRTbell™接头（adapter，图6中黑域）的子序列（接头间的序列）个数，图6中fullpasses2CircularConsensus(CCS)CCS序列是通过每个ZMW孔中的（subreads）子序列得到的一致性序列，无需进行参考序列的比对。不同于ROI序列，序列要求每个插入序列中至少含2个完整（full-pass）的subreadsReadsofInsert根据需求筛选原始序列中满足最小numberoffullpasses（<低序列准确性（0.70~1.00）的原始序列，通过子序列获得一致序列即为每个单分子反应器ZMW的ROI序列。ROI是每个单分子反应器ZMW里插入序列的最高质量序列。ROI序列可以从一定程度上评估建库质量和SMRT®Cell上样时序列的长度。为充分利用数据量，本项目不考虑完整的subread个数（即设置最小fullpasses=0），最低序列准确度为0.75。（详细内容请参见：/acificiocince/c_prier/iki/ndersanding-acio-tansripoe-daa#radxplined）。例如，如果根据模板链最终合成了一个半的subread，最终将合并成一条插入序列（ReadofInsert，ROI），而CCS序列要求最少有两个完整的subreads，所以它是ROI序列的特殊情况。ROI序列可以很好的评估SMRTcell上样时的插入序列的长度。对于较长的模板链，最终产生的ROI序列可能和polymeraseread序列等长。Full-Length(FL)ReadversusNon-Full-Length(nFL)SMRTLink软件中IsoSeq2定义两端同时含有3'引物和5'3'引物前polyA尾（可选）的序列称为全长序列（Full-Length(FL)Read）。5'或3'primer可以是Clontech或其他全长cDNA建库引物，或特异性的RT-PCR引物。反之，则为非全（non-full-lengthread）Full-Lengthnon-chimericRead(FLNC)：建库过程中因接头浓度或SMRTbell浓度过低造成两个cDNA模板链直接相连而生的嵌合序列称为人工嵌合序列，如下图所示全长序列中的人工嵌合序SMRT所使用的PacBio仪器每个cell含有ZMWs，reads进入ZMW孔中被，一个ZMW中含一条的reads（P1）为有据本项目F01下机数据过滤PolymeraseRead片段长度小于50bp、序列准确性小于0.75的序列，将剩余序列从接头处打断并过滤掉接头序列后得到subreads，过滤长度小于50bp的sbred，剩余subred即为cleandata。 数据统计SampleSamplecDNASMRTData1-22-13-1根据条件fullpasses>=0且序列准确性大于0.75从原始序列提取ROI序列。统计各建库中ROI（ReadsofInsert）序列数、ROI的碱基数和插入序列的平均长度对下机数据进行评估。该项目各建库数据产出统计见下表:ROI数据统计部分结果展cDNAReadsofReadBasesofMeanReadLengthofMeanReadQualityof1-2-3-(1)ROI序列长度分布:ROI序列长度与建库时的cDNA长度选择有关，ROI序列随cDNA长度的增加而增加cDNA长度建库ROI序列长度分布如下图所示cDNA建库ROI长度分布示意注:横坐标:ROI长度分布;左侧纵坐标：ROI序列长度频数分布直方图；右侧纵坐标：ROI序列长度累积频率曲线(3)ROIfullpasses数分布ROI的各孔fullpassescDNA的长度有关，一般cDNA长度的增加而减少。ROI序列的准确性fullpasses的影响，fullpasses越高，序列准确性越高，各cDNA建库数据的fullpasses数分布如下图所示:生成ROI序列的fullpasses数分布示注:横坐标:生成ROI序列的fullpasses数;纵坐标:相应fullpasses数的ROI序列数(3)ROI质量分布 质量值反映序列的准确性，各建库数据的ROI序列质量分布如下图所示ROI质量值分布示意注:横坐标:ROI质量值分布;左侧纵坐标:ROI序列质量频数分布直方图;右侧纵坐标:大于相应质量值的碱基通过检测ROI序列中是否包含5'引物，3'引物及polyA尾，可将序列分成全长序列（包含5'引物，3'引物及polyA尾）和非全长序列。classify过程去除ROI序列中cDNA引物序列及polyA序列获得建库时的插入序列，同时根据建库时两端引物的差别确定链合成方向，并将序列分为全长序列和非全长序列、嵌合序列和非嵌合序列。全长序列数据统计cDNAReadsofNumberoffiveprimeNumberofthreeprimeNumberofpoly-A1-2-3-全长序列长度反映了建库时cDA序列长度，可通过统计全长序列的长度评估建库质量，本次数据中的全长序列长度与建库大小一致。不同cDNA建库时FLNC序列的长度分布如下图所示：全长非嵌合序列长度分布示意注:横坐标:FNC序列长度分布;左侧纵坐标：FNC序列长度频数分布直方图；右侧纵坐标：FLNC序列长度累积频率曲线。LC长度除受建库时cDA的长度影响外，还受试剂和时间等因素的影响，故全长序列的百分比（L%）在不同cell中会有很大变化。因建库过程中SMRTbell或adaptor浓度过低，会产生嵌合（ArtificialConcatemers）序列，本次全长序列中ArtificialConcatemers比例0.38%，可认为该数SMRTbell浓度适中。各库中所有ROI序列分类结果如下图所示:序列分类（classify）示意isoformSMRTysis(v2.3.0)软件使用ICE（IterativeClusteringforErrorCorrection）算法对序列进行迭代聚类，将相似的序列（即同一转录本的多个拷贝）聚类到一簇cluster，每个cluster得到一个一致性序列（consensusisoform）。结合非全长序列，使用quiver程序对各cluster中的一致序列进行校正,最终得到准确度大于99%的高质量转录（HQ，high- isoforms）。不同建库长度得到的HQ和LQ(低质量转录本，low-quality)转录本统计结果见下表ICE聚类结果统计NumberofconsensusAverageconsensusisoformsreadNumberofpolishedhigh-quality本项目F01得到consensusisoform65,257个，其中HQ转录本数为24,889，LQ转录本40,345个，得到的consensus序列长度分布如下所示consensusisoform长度分布示注:横坐标:consensusisoform序列长度分布;左侧纵坐标：consensusisoform序列长度频数分布直方图；右侧纵坐consensusisoform序列长度累积频率曲线利用proovread软件[5]通过RNA-Seq数据对各个样品得到的低质量一致性序列进行校正，提高序列的准确性。转录本去冗在全长转录本cluster过程中参数设置较严格，为得到质量较高的一致性序列，同一转录本的多拷贝序列分到不同cluster的可能性比随机将两个不属于同一转录本拷贝的序列分到同一cluster中的可能性大，不可避免的产生了冗余序列。同时，全长转录本过程中，3'端因存在poly-A结构，可以确定3'端比较完整，而5'端序列可能存在降解，导致同一转录本的不同拷贝分到不同的cluster中，如下图绿色圈中所示，5'端差异造成不同转录本，导致冗余序列的产生。5'序列差异造成转录本序列通过GMAP（GenomicMapandAlignmentProgram）[6]将得到的校正后的一致序列与参考组进行序列比对（设置参数--cross-species--allow-close-indels0），使用cDNA_Cupcake（ /Magdoll/cDNA_Cupcake/wiki）软件对比对结果去冗余，过滤identity小于0.9，coverage小于0.85的序列，合并仅5’端外显子有差异的比对。过滤的转录本列1.转录本序列1.转录本分类结1.转录本注释结果1.转录组eads与参考组序列比对结果文件（通常为BAM格式）、物种参考组序列和注释文件，推荐使用整合组浏览器（IGV，IntegrativeGenomicsViewer）进行可视化浏览。IGV具有以下特点:能在不同尺度下显示单个或多个eads在参考组上的位置，包括eads在各个上的分布情况和在注释的外显子、内含子、剪接接合区、间区的分布情况等;能在不同尺度下显示不同区域的Reads丰度，以反映不同区域的转录水平能显示及其剪接异构体的注释信息能显示其他注释信息既可以从服务器端各种注释信息，又可以从本地加载注释信息IGV浏览器界利用BUSCO[7]对去冗余后的转录组进行完整性评估，评估结果如图所示转录组完整性评估结果示意转录本融合发因为去冗余分析过程中，通过设置最小coverage为0.85过滤了嵌合序列，故融合转录本分析中需要使用去冗余前的一致序列，按下条件进行融合转录本筛选比对到2个或多个位点每个位点必须比对上至少5%的转录本长度,最小比对长度为所有位点比对到的总长度必须占转录本总长度的95%以上两位点间距离必须达10kb以上。gff注释文件:融合转录本GFF文件1.提供组注释补充信息的同时，也提供以FATA格式的新转录本序列。融合转录本序列的FSTA文件:融合转录本FASTA文件1.转录本结构分转录生成的前体mN（pre-RA），有多种剪接方式，选择不同的外显子，产生不同的成熟mN，从而翻译为不同的蛋白质，构成生物性状的多样性。这种转录后的mNA加工过程称为可变剪接或选择性剪接（lternativespliing）。通过Astalavista软件[8]获取每个样品存在的可变剪接类型，主要的可变剪接类型如下图所示可变剪接类注：(A外显子跳跃；(B)可变转录终止位点；(C)可变外显子；(D)可变转录起始位点；(E内含子保留我们从Astalavista软件分析结果中,对转录本发生上述5种可变剪接事件情况进行统计,各样品中预测的可变剪接事件数量统计见下可变剪接事件数量统注：Alternative3'splicesite:可变转录终止位点；Alternative5'splicesite:可变转录起始位点；Exonskip外显子跳跃；Intronretention内含子保留；Mutuallyexclusiveexon:可变外显子可变剪接事件统计1.由于使用的软件或数据本身的局限性，导致所选参考组的注释往往不够精确，这样就有必要对原有注释的结构进行优化。如果在原有边界之外的区域有MappedReads支持，将的非翻译区（UntranslatedRegion，UTR）向上下游延伸，修正的边界。结构优化结果见下面文件：结构优化结注：GeneID：ID；Locus：座，格式为“编号:起点坐标-终点坐标”；Strand：正负链；Site：优化的位置，3'或5'UTR；OriginalSite：原来注多聚腺苷酸化是指多聚腺苷酸与信使RNA（mRNA）分子的共价链结。在蛋白质生物合成的过程中，这是产生准备作翻译的成mRNA的方式的一部份。在真核生物中，多聚腺苷酸化是一种机制，令mRNA分子于它们的3'端中断。多聚腺苷酸尾（A尾）保护mRNA，免受核酸外切酶，并且对转录终结、将mRNA从细胞核输出及进行翻译都十分重要。在原核生物中，前体mRNA的可变多聚腺苷酸化（alternativepolyadenylation,APA）可能贡献于转录组多样性，组的编码能力以及的调控机制。我们采用TAPISpipeline[9]来识APA。各个样品所识别出的APA如下：可变多聚腺苷酸化结1.多聚腺苷酸化位点个数分布示意注:横坐标:多聚腺苷酸化位点个数;纵坐标：个数利用MEME对所有转录本polyA位点上游50bp的序列进行分析，鉴定得到的motif如下图所示polyA位点上游motifSSR分MISA(MIcroSAliteidentificationtool)是一款鉴定简单重复序列的软件，其参考见附表。它可以通过对转录本序列的分析，鉴定出7种类型的SSR:Mono-nucleotide（单碱基）、Di-nucleotide（双碱基）、Tri-nucleotide（三碱基）、Tetra-nucleotide（四碱基）、enta-nuceotde（五碱基）、exa-nuclotie（六碱基）、copoundR（混合微，两个SR距离小于100bp）从新转录本中筛选500bp以上的转录本，利用MISA软件做SSR分析，结果见下表SSR分析结果统SearchingTotalnumberofsequencesTotalsizeofexaminedsequencesTotalnumberofidentifiedNumberofSSRcontainingNumberofsequencescontainingmorethan1NumberofSSRspresentincompoundMonoDiTriTetraPentaHexaSSR分析结注：每一行为一个SSR标记及其引物设计结果，每个SSR标记最多设计三对引物;eI：转录本编号；SSn：同一转录本上的SSR序号；SSty：SSR类型，包括完美单碱基重复（）、完美双碱基重复（p）、完美三碱基重复（p）、完美四碱基重复（p）、完美五碱基重复（）、完美六碱基重复（6）和混合SS（c，即包含至少两个完美SS，且之间距离小于1b）；SSR：SSR序列，括号内为重复单元，括号外数字表示重复次数；Siz：SSR的长度；SSStat：SSR在转录本上的开始位置；SSREd：SSR在转录本上的结束位置;FP1'-：第一条正向引物序列；Tm：第一条正向引物序列的退火温度，单位为C；Siz：第一条正向引物序列的长度；P1'-')：第一条反向引物序列；Tm：第一条反向引物序列的退火温度，单位为C；Siz：第一条反向引物序列的长度；Psiz：产物的长度；PStt：产物在上的开始位置；PEd：产物在上的结束位置。对不同SSR类型的密度分布进行统计，结果见下图SSR类型分布示意新编码区序列预TransDecoder[10](v3.0.0)软件基于开放阅读框（OpenReadingFrame，ORF）长度、对数似然函数值（Log-likelihoodcore）、氨基酸序列与fam数据库蛋白质结构域序列的比对等信息，能够从转录本序列中识别可靠的潜在编码区序列（odingequence，CDS），是常用的CDS预测软件对可变剪接分析中得到的新转录本使用TransDecoder软件对其编码区序列及其对应氨基酸序列的预测。本次共获得orf个，其中完整orf20,456条，预测的CDS文件结果如下图CDS分析结果文件示意注:文件为标准的FASTA格式，每个序列单元以开始到下一个“<”之前结束。“<”后面紧接编码区序列编号（由转录本编号和开放阅读框序号构成），转录本编号，之后是序列描述信息，包括序列类型（OFty）、序列长度（l）、序列在转录本上的位置（lc）和链类型。预测的编码区序列类型有完整（cmlt，即同时预测到起始子和终止子）、’端部分（im_til，即仅预测到起始子）、’端部分（pim_atil，即仅预测到终止子）和内部区段（itrl，即起始子和终止子都没有预测到）。从第二行开始直到出现之前为止为编码区序列。预测得到的完整ORF区编码蛋白序列长度分布如下预测的CDS编码蛋白长度分布示意LncRNA预因lncA不编码蛋白，因此，通过对转录本进行编码潜能筛选，判断其是否具有编码潜能，从而可以判定该转录本是否为lncRNA。百迈客综合目前应用最广泛的编码潜能分析方法对新发现的转录本进行lncRNA的预测，主要包括:CPC[11]分析、CNCI[12]分析、pfam蛋白结构域分析、CPAT[13]分析四种方法。CPC(CodingPotentialCalculator)是一种基于序列比对的蛋白质编码潜能计算工具。通过将转录本与已知蛋白数据库比对根据转录本各个编码框的生物学序列特征评估其编码潜能。Score<0时，为noncodingRNA。CPC分析结果CPC分析结果统注:第1列:transcript_id:转录本ID;第2列 length:ORF长度;第3列:type:转录本类型;第4列:score:转录本得分，当score＜0时，为NoncodingCNCICNCI(Coding-Non-CodingIndex)分析是一种通过相邻核苷酸三联体特征区分编码-非编码转录本的方法。该工具不依赖于已知的注释文件，可以有效对不完整的转录本和反义转录本进行预测。NI工具提供两种比对模式：ve（脊椎物种）；pl（植物物种），本项目选择pl参数。当score<0时，为 CI分析结果如下:CNCI分析结果统注:第1列:transcript_id:转录本ID;第2列:type:转录本类型;第3列:score:转录本得分，当score<0时，为Noncoding;第4列:start:转录本起始位置;第5列:end:转录本终止位置。CPATCPAT(CodingPotentialAssessmentTool)分析是一种通过构建逻辑回归模型，基于ORF长度、ORF覆盖度，计算Fickett得分和Hexamer得分来判断转录本编码和非编码能力的分析方法。CPAT分析结果如下:CPAT分析结果统fam数据库是最全面的蛋白结构域注释的分类系统。蛋白质是由一个或多个结构域组成的，而每个特定结构域的蛋白序列具有一定保守性。fam将蛋白质的结构域分为不同的蛋白，通过蛋白序列的比对建立了每个的氨基酸序列的HM统计模型。百迈客将转录本各个编码框上的蛋白序列与pfam数据库做hmscan同源搜索，能比对上的转录本即为具有某个蛋白结构域的转录本，即认为具有编码能力，而无比对结果的转录本被认为是潜在的lncRNA。结果如下:pfam分析结果统注:第1列:transcript_id:转录本ID;第2列hmmacc:比对到pfam结构域ID;:hmmstart:比对到结构域的起始位置;第3列:hmmname:pfam结构称;第4列:hmmstart:比对到结构域的起始位置;5列:hmmend:比对到结构域的终止位置;6列:hmmlength:pfam结构域的长度;7列bitscore比对打分值;第八列:E-value:比对的E值，pfam结构域筛选的条件E-value<0.001。为直观展示分析结果，将以上4种分析软件鉴定得到的noncdingtranscripts进上4种分析结果取交集，用于后续lncRNA4种筛选方法维恩根据lncRNA在参考组注释信息（gff）上的位置，对lncRNA进行分类绘图，结果如下lncRNA位置分类示lncRNA靶预对预测得到的lncRNA序列进行靶预测。基于lncRNA与其靶的作用方式，我们采用2种预测方法第一种，lncRNA调控其邻近的表达，主要根据lncRNA与mRNA的位置关系预测，定义中每100kbp范围内存在差异表达lncRNA与差异表达的mRNA;第二种，lncA与mNA由于碱基互补配对而产生作用，主要利用Lncar[14]靶预测工具对我们的LncRA进行靶预测。两种方法共预测得到的靶结果请见文件：靶预测结果注:第一列:lncRNA的id号;第二列:对应的编号lncRNA；靶id号;基于位置关系的LncRNA靶预基于物理位置的靶预测结注:第一列:lncRNA的id号;第二列:lncRNA对应的靶id号;编号基于互补序列的靶预注:#LncRNA_ID:lncRNA的id号;TargetGene_id:lncRNA对应的靶ID转录因子分转录因子（Transcriptionfactor）是指能够结合在某上游特异核苷酸序列上的蛋白质，这些蛋白质可以调控RNA聚合酶与DNA模板的结合，从而调控的转录。植物转录因子预测使用iTAK[15]软件，动物转录因子鉴定使用动物转录因子数据库—— 2.0[16]，共预测得到转录因子4,005个，预测结果如下:转录因子序列文使用BLAST[17]软件(version2.2.26)将得到的新转录本序列与NR[18]、Swissprot[19]，GO[20]、COG[21]、KOG[22]、Pfam[23]、NR数据库是NCBI中的非冗余蛋白质数据库，包含了Swissprot、PIR(ProteinInformationResource)、PRF(ProteinResearchFoundation)、PDB(ProteinDataBank)蛋白质数据库及从GenBank和RefSeq的CDS数据翻译过来的蛋白质数据信息。该数据库网址见附表。issprot数据库是由I（欧洲生物信息学）负责的数据库，包含了有相关参考文献且经过校对的蛋白质注释信息数据库，可信度很高。COG(ClustersofOrthologousGroups)数据库是对产物进行同源分类的数据库，是一个较早的识别直系同源的数据库，通过对多种生物的蛋白质序列大量比较而来的。KOG（euKaryoticOrthologGroups）数据库是针对真核生物，基于直系同源关系，结合进化关系将来自不同物种的同源分为不同的Ortholog簇。来自同一Ortholog的具有相同的功能，这样就可以将功能注释直接继承给同一KOG簇的其他成员。Pfam(Proteinfamily)是最全面的蛋白结构域注释的分类系统。蛋白质是由一个个结构域组成的，而每个特定结构域的蛋白序列具有一定保守性。Pfam将蛋白质的结构域分为不同的蛋白，通过蛋白序列的比对建立了每个的氨基酸序列的HMM统计模型。GO(GeneOntology)数据库是一个国际标准化的功能分类体系，提供了一套动态更新的标准词汇表来全面描述生物体中和产物的功能属性。该数据库总共有三大类，分别是分子功能(molecularfunction)，细胞组分(cellularcomponent)和生物学过程(biologicalprocess)，各自描述了产物可能行使的分子功能，以及所处的细胞环境和参与的生物学过程。GO数据库中最基本的概念是Term，每个条目都有一个Term名，比如“cell”、“fibroblastgrowthfactorreceptorbinding”或者“signaltransduction”，同时有一个唯一的编号，形如GO:nnnnnnn。KEGG(KyotoEncyclopediaofGenesandGenomes)数据库是系统分析产物在细胞中的代谢途径以及这些产物功能的数据库。它整合了组、化学分子和生化系统等方面的数据，包括代谢通路(AHWA)、药物(DU)、疾病(IES)、序列(GENES)及组(GENOME)等。利用该数据库有助于把及表达信息作为一个整体的网络进行研究。将可变剪接分析中得到的新转录本进行功能注释，各数据库注释到的转录本数量统计见下表注释的转录本数量统计AnnotatedNewIsoformSwiss-(genefamily)是来源于同一个祖先，由一个通过重复而产生两个或的拷贝而构成的一组，它们在结构和功能上具有明显的相似性，编码相似的蛋白质产物。通过Pfam注释，转录本的 分类结果如下:分类结注:第一列:名称;第二列:对应转录本的idSNP/InDel分SNP（SingleNucleotidePolymorphisms）是指在组上由单个核苷酸变异形成的遗传标记，其数量很多，多态性丰富。百迈客基于各样品reads与参考组序列的opat2比对结果，使用AK软件[25]识别样品与参考组间的单碱基错配，识别潜在的SNP位点。进而可以分析这些SNP位点是否影响了的表达水平或者蛋白产物的种类Inel(inertin-eleton)是指相对于参考组，样本中发生的小片段的插入缺失，该插入缺失可能含一个或多个碱基。GAK也能够检测样品的插入缺失（Inel）。Inel变异一般比SP变异少，同样反映了样品与参考组之间的差异，并且编码区的Inel会引起移码突变，导致功能上的变化。GATK识别标准如下：35bp范围内连续出现的单碱基错配不超过3个经过序列深度标准化的SNP质量值大于2.0。各样品分别按照以上条件筛选，最终获得可靠的SNP位点SnpEff[26]是一款用于注释变异（SNP、InDel）和预测变异影响的软件。根据变异位点在参考组上的位置以及参考组上的位置信息，可以得到变异位点在组发生的区域（间区、区或DS区等），以及变异产生的影响（同义非同义突变等）。由于转录完成之后，RA除了需要加帽、加loy()和可变剪接之外，较少mNA会经历RA编辑（RAeditin），从而会产生单碱基的替换、插入、缺失。NA编辑能使同一产生序列多样的mRA，但是这种多态性不是组固有的多态性。从比对结果来看，SNP和单碱基替换的RNA编辑结果是一样的。因此，通过转录组数据识别出SNP不免会含有RNA编辑的产物。SNP位点信InDel位点信注：Chr：SNP/InDel位点所在 编号；Pos：SNP/InDel位点在 上的位置；Gene_id：SNP/InDel位点所在的或原来未注释的区（表中用Intergenic表示）；Ref：所选参考 组中的SNP/InDel等位；Alt： 样品中识别到的其他的SNP/InDel等位；T*：样品T*该SNP/InDel位点的分型；Depth：样品T*该SNP/InDel位点的 深度；AlleDp：样品T*该SNP/InDel位点的各等位 深度；Effect：SNP/InDel所在区域或类型；Codon_change：编码改变方式，未改变用点表示。Effect具体说明详见： 根据SP位点碱基替换的不同方式，可以将SP位点分为转换（ransition）和颠换（ransvesion）两种类型。根据SP位点的等位（Allele）数目，可以将SNP位点分为纯合型SNP位点（只有一个等位）和杂合型SNP位点（两个或多个等位）。不同物种杂合型SNP所占的比例存在差异。对各样品筛选出的SNP位点数目、转换类型比例、颠换类型比例以及杂合型SNP位点比例进行统计，如下表SNP位点统计部分结果展BMK-SNPGenicIntergenicSNP突变类型统计分布如下图所示SNP突变类型分布示意注：横轴为SNP突变类型，纵轴为相应的SNP将每个的SP位点数目除以的长度，得到每个的SP位点密度值，统计所有的NP位点密度值并做密度分布图。部分的SNP位点密度分布图如下：SNP密度分布注：横轴为上平均每1000bp序列中分布的SNP数目，纵轴为数采用SNPEff分别对SNP，InDel注释，SNP，InDel的注释结果统计如下所示注：纵轴为SNP所在区域或类型，横轴为分类数目。InDel注释分类示注：纵轴为InDel所在区域或类型，横轴为分类数目本项目使用指定的组作为参考进行序列比对及后续分析，利用HISAT2将CleanReads与参考组进行序列比对，获取在参考组或上的位置信息。将比对到不同上的Reads进行位置分布统计，绘制MappedReads在所选参考组上的覆盖深度分布图二代数据与参考组比对部分结果展TotalMappedUniquelymappedMultiplemappedMappedReads在参考组上的位置及覆盖深度分布示意注：横坐标为位置；纵坐标为覆盖深度以2为底的对数值，以10kb作为区间单位长度，划分成多个小窗口（Window），统计落在各个窗口内的Mapped Reads作为其覆盖深度。蓝色为正链，绿色为负链。统计MappedReads在指定的参考组不同区域（外显子、内含子和间区）的数目，绘制组不同区域上各样品Mapped的分布图，如下组不同区域Reads分布示意注：图中将组分为外显子区、间区、内含子区，区域大小按Map到相应区域的Reads在所有MappedReads中所占理论上，来自成熟RA的eads应比对到外显子区。eads比对到内含子是由于RA前体和发生可变剪切的内含子保留；Reads比对到间区是由于组注释不完善。合格的转录组文库是转录组的必要条件，为确保文库的质量，从以下3个不同角度对转录组文库进行质量评估。(1)片段化随机性检验mNA片段化后的插入片段大小选择，是从RA序列中独立随机地抽取子序列，mNA数目越大、打断方式和时间控制得越合适，目的RA每个部分被抽取到的可能性就越接近，RA片段化随机性越高，mRA上覆盖的eads越均匀。通过apped eds在各mRNA转录本上的位置分布，模拟mRNA片段化结果，检验mRNA片段化的随机程度。如果mRNA存在严重降解，被降解的碱基序列不能被，即无eads比对上。因此，通过查看appedeads在NA转录本上的位置分布可了解RA的降解情况。样品MappedReads在mRNA转录本上的位置分布如下图：MappedReadsmRNA上的位置分布示注：横坐标为标准化后的mRNA位置，纵坐标为对应位置区间内Reads在总MappedReads中所占百分比。由于参考的mRNA长度不同，作图时对把每个mRNA按照长度划分成100个区间，进而统计每一区间内的MappedReads数目及所占的比例，图中反映的是所有mRNA各个区间内的MappedReads比例的汇总。插入片段长度检验插入片段长度检验插入片段长度的离散程度能直接反映出文库过程中磁珠纯化的效果。通过插入片段两端的eads在参考组上的比对起止点之间的距离计算插入片段长度。大部分的真核生物为断裂，外显子被内含子隔断，而转录组得到的是无内含子的成熟RN。当RA中跨内含子的片段两端的eads比对到组上时，比对起止点之间的距离要大于插入片段长度。因此，在插入片段长度模拟分布图中，主峰右侧形成1个或多个杂峰。各样品的插入片段长度模拟分布图如下：插入片段长度模拟分布示意注：横坐标为双端Reads在参考组上的比对起止点之间的距离，范围为0到800bp；纵坐标为比对起止点之间不同距离的双端或插入片段数量转录组数据饱和度检验为了评估数据是否充足并满足后续分析，对得到的转录本数进行饱和度检测。由于一个物种的数目是有限的，且转录具有时间和空间特异性，因此随着量的增加，检测到的转录本数目会趋于饱和。对于表达量越高的转录本，越容易被检测定量。因此，对于表达量越低的转录本，需要更大的数据量才能被准确定量。使用各样品的appedata对检测到的不同表达情况的转录本数目饱和情况进行模拟，绘制曲线图如下：转录组数据饱和度模拟示意注：本图为随机抽取10%、20%、30%……90%的总体数据单独进行定量分析的结果；横坐标代表抽取数据定位到组上Reads数占总定位的reads数的百分比，纵坐标代表所有抽样结果中表达量差距小于15%的Gene在各个FPKM范围的百分比转录组可以模拟成一个随机抽样的过程，即从一个样品转录组的任意一段核酸序列上独立地随机抽取序列片段。抽取自某一（或转录本）的片段数目服从负二项分布（BetaNegativeBinomialDistribution）。基于该数学模型，使用Cufflinks软件[27]的Cuffquant和Cuffnorm组件，通过MappedReads在上的位置信息，对转录本的表达水平进行定量。抽取自一个转录本的片段数目与数据（或apped ata）量、转录本长度、转录本表达水平都有关，为了让片段数目能反映转录本表达水平，需要对样品中的appedeads的数目和转录本长度进行归一化。uffuant和uffnorm采用FPKM（FragmentsPerKilobaseoftranscriptperMillionfragmentsmapped）作为衡量转录本或表达水平的指标，FPKM计算公式如下：FPKM计算公注：公式中，cDNAFragments表示比对到某一转录本上的片段数目，即双端Reads数目；MappedFragmentsMillions)表示比对到转录本上的片段总数，以10^6为单位；TranscriptLength(kb)：转录本长度10^3个碱基为单位。所有转录本的注：第一列表示转录本ID,其余各列表示转录本在不同样品中的表达量样品转录本表达量总体分各样品FPKM密度分布对比注：图中不同颜色的曲线代表不同的样品，曲线上点的横坐标表示对应样品FPKM的对数值，点的纵坐标表示概率密为了进一步查看单个样品转录本表达水平分布的离散程度，以及直观的比较不同样品的整体转录本表达水平。采用箱线图展示FPKM分各样品FPKM箱线注：图中横坐标代表不同的样品；纵坐标表示样品表达量FPKM的对数值。该图从表达量的总体离散角度来衡量各样品表达水平。对于设立生物学重复的项目，评估生物学重复的相关性对于分析转录组数据非常重要。生物学重复的相关性不仅可以检验生物学实验操作的可重复性；还可以评估差异表达的可靠性和辅助异常样品的筛查。将皮尔逊相关系数r（PearsonCorrelationCoefficient）作为生物学重复相关性的评估指标。r2越接近1，说明两个重复样品相关性越强。百迈客保证对同一条件的所有生物学重复样品进行同人同批样品提取、建库，同Run同Lane。对异常样品进行详细分析，并根据分析结果与沟通共识决定重新进行实验，还是剔除异常样品进行后续分析。为研究在不同实验处理下的表达模式，首先对选取样品的FPKM进行标准化归一处理，然后对标准化的FPKM值进行K-聚类分析，同一类的在不同的实验处理下具有相似的变化趋势，具有相似变化趋势的往往具有相似的功能转录本表达聚类折线示意注：x轴表示实验分组，y轴表示标准化的FPKM。黑色线条表示这个cluster中的所有在不同实验条件下相表达量的平均值的折线图差异表达分表达具有时间和空间特异性，在两个不同条件下，表达水平存在显著差异的，差异表达（EG）。差异表达分析得到的基因集合叫做差异表达集，使用“_vs_”的方式命名。根据两（组）样品之间表达水平的相对高低，差异表达可以划分为上调（Up-regulatedTranscript）和下调（Down-regulatedTranscript）。上调在样品（组）B中的表达水平高于样（组）A中的表达水平；反之为下调。上调和下调是相对的，由所给A和B的顺序决定对于有生物学重复的样本，Eeq适用于进行样品组间的差异表达分析，获得两个生物学条件之间的差异表达集；对于没有生物学重复的样本，使用EBSeq进行差异分析。在差异表达过程中，将FoldChange≥2且FDR<0.01作为筛选标准。差异倍数（FoldChange）表示两样品（组）间表达量的比值。错误发现率（alseiscoeryate，R）是通过对差异显著性p值（p-value）进行校正得到的。由于转录组的差异表达分析是对大量的表达值进行独立的统计假设检验，会存在假阳性问题，因此在进行差异表达分析过程中，采用了公认的Benjamini-Hochberg校正方法对原有假设检验得到的显著性p值（p-value）进行校正，并最终采用FDR作为差异表达筛选的关键差异表达集数目统计如下表差异表达数目统计部分结果展DEGDEGup-down-通过火山图（VolcanoPlot）可以快速地查看在两个（组）样品中表达水平的差异，以及差异的统计学显著性。差异表达火山图下图差异表达火山示意注：差异表达火山图中的每一个点表示一个，横坐标表示某一个在两样品中表达量差异倍数的对数值；纵坐标表示表达量变化的统计学显著性的负对数值。横坐标绝对值越大，说明表达量在两样品间的表达量倍数差异越大；纵坐标值越大，表明差异表达越显著，筛选得到的差异表达越可靠。图中绿色的点代表下调差异表达，红色的点代表上调差通过MA图可以直观地查看的两个（组）样品的表达水平和差异倍数的整体分布。差异表达MA图见下图差异表达MA示意注：差异表达MA图中每一个点代表一个。横坐标为A值：log2(FPKM)，即两样品中表达量均值的对数值；纵坐标为M值：log2(FC)，即两样品间表达量差异倍数的对数值，用于衡量表达量差异的大小。图中绿色的点代表下调差异表达，红色的点代表上调差异表达，黑色的点代表非差异表达。对筛选出的差异表达做层次聚类分析，将具有相同或相似表达模式的进行聚类，部分差异表达聚类结果如下图差异差异表 聚类部分结果展注：横坐标代表样品名称及样品的聚类结果，纵坐标代表的差异及的聚类结果。图中不同的列代表不同的样品，不同的行代表不同的。颜色代表了在样品中的表达量水平log2（FPKM+1）。对差异表达进行数据库的功能注释，各差异表达集注释到的数量统计见下表注释的差异表达数量部分结果展DEGSwiss-差异表达GO注释及富GO数据库是GO组织（GeneOntologyConsortium）于2000年构建的一个结构化的标准生物学注释系统，旨在建立及其产物知识的标准词汇体系，适用于各个物种。GO注释系统是一个有向无环图，包含三个主要分支，即：生物学过程（BiologicalProcess），分子功能（MolecularFunction）和细胞组分（CellularComponent）。差异表达GO分类统计结果见下图差异表达GO注释分类统计示意注：横坐标为GO分类，纵坐标左边为数目所占百分比，右边为数目。此图展示的是在差异表达背景和全部背景下GO各二级功能的富集情况，体现两个背景下各二级功能的地位，具有明显比例差异的二级功能说明差异差异表达topGO有向无环示意注：对每个GO节点进行富集，最显著的10个节点在图中用方框表示，图中还包含其各层对应关系。每个方（或椭圆）内给出了该GO节点的内容描述和富集显著性值。不同颜色代表不同的富集显著性，颜色越深，显著性越差异表达KEGG注释及富在生物体内，不同的产物相互协调来行使生物学功能，对差异表达的通路（Pathway）注释分析有助于进一步解读的功能。KEGG（KyotoEncyclopediaofGenesandGenomes）是系统分析功能、组信息数据库，它有助于研究者把及表达信息作为一个整体网络进行研究。作为是有关Pathway的主要公共数据库(Kanehisa,2008），KEGG提供的整合代谢途径(pathay)查询，包括碳水化合物、核苷、氨基酸等的代谢及有机物的生物降解，不仅提供了所有可能的代谢途径，而且对催化各步反应的酶进行了全面的注解，包含有氨基酸序列、PDB库的等等，是进行生物体内代谢分析、代谢网络研究的强有力工具。差异表 的通路注释结果见下图差异表达KEGG通路注释示注：相对于对照组来说，红色框标记的酶与上调因有关，绿色框标记的酶与下调有关。蓝色框标记的酶与上调和下调均有关，框内的数字代表酶的编号（ECnumber），而整个通路由多种酶催化的复杂生化反应构成，此通路图中与差异表达相关的酶均用不同的颜色标出，根据研究对象间的差异，重点研究某些代谢通路相关的对差异表达KEGG的注释结果按照KEGG中通路类型进行分类，分类图如下图所示差异表达KEGG分类示意注：纵坐标为KEGG代谢通路的名称，横坐标为注释到该通路下的个数及其个数占被注释上的总数的比例分析差异表达在某一通是否发生显著差异（over-presentation）即为差异表达的通路富集分析。Pathway显著性富集分析以KEGG数据库中Pathway为单位，应用超几何检验，找出与整个背景相比，在差异表达中显著性富集的Pathway。差异表达KEGG通路富集分析结果见下图，图中呈现了显著性Q值最小的前20个通路。差异表达KEGG通路富集散点示意注：图中每一个圆表示一个KEGG通路，纵坐标表示通路名称，横坐标为富集因子（EnrientFactor），表示差异中注释到某通路的比例与所有中注释到该通路的比例的比值。富集因子越大，表示差异表达在该通路中的富集水平越显著。圆圈的颜色代表qvalue，qvalue为多重假设检验校正之后的Pvalue,qvalue越小，表示差异表达在该通路中的富集显著性越可靠；圆圈的大小表示通路中富集的数目，圆圈越大，表示越多。差异表达COG分COG（ClusterofOrthologousGroupsofproteins）数据库是基于细菌、藻类、真核生物的系统进化关系构建得到的，利用COG数据库可以对产物进行直系同源分类。差异表达COG分类统计结果见下图差异表达COG注释分类统计示意注：横坐标为COG各分类内容，纵坐标为数目。在不同的功能类中，所占多少反映对应时期和环境下代谢或者生理偏向等内容，差异表达eggNOG分eggNOG（evolutionarygenealogyofgenes:non-supervisedorthologousgroups）数据库对直系同源类群进行了功能描述和功能分类的注释,包含1133个物种的直系同源类群。eggNOG分类统计结果见下图:差异表达eggNOG注释分类统计示意图注：横坐标为eggNOG各分类内容，纵坐标为数目。差异表达蛋白互作网TIG是收录多个物种预测的和实验验证的蛋白质-蛋白质互作的数据库，包括直接的物理互作和间接的功能相关。结合差异表达分析结果和数据库收录的互作关系对，构建差异表达互作网络。对于数据库中包含的物种，可直接从数据库中提取出目标集的互作关系对构建互作网络；对于数据库中未收录信息的物种，使用LAT软件，将目的与数据库中的蛋白质进行序列比对，寻找同源蛋白，根据同源蛋白的互作关系对构建互作网络。构建完成的蛋白质互作网络可导入Cytoscape软件进行可视化。Cytoscape可视化的差异表达蛋白质互作网络如下图差异表达蛋白质互作网络示意注：图中的节点为蛋白质，边为互作关系。互作网络点(node)的大小与此节点的度(degree)成正比，即与此节点相连的边越多，它的度越大，节点也就越大。节点的颜色与此节点的系数(clusteringcoefficient)相关，颜色梯度由绿到红对应系数的值由低到高；系数表示此节点的邻接点之间的连通性好坏，系数值越高表示此节点的邻接点之间的连通性越好。边(edge)的宽度表示此边连接的两个节点间的互相作用的关系强弱，互相作用的关系越参考文AuKF,SebastianoV,AfsharPT,DurruthyJD,LeeL,WilliamsBA,etal.CharacterizationofthehumanESCtranscriptomebyhybridsequencing.ProcNatlAcadSciUSA.2013;110(50):E4821–30. :10.1073/pnas. SharonD,TilgnerH,GrubertF,SnyderM.Asingle-moleculelong-readsurveyofthehumantranscriptome.NatBiotechnol.2013;31(11):1009–14. :10.1038/nbt.2705.TreutleinB,GokceO,QuakeSR,SudhofTC.Cartographyofneurexinalternativesplicingmappedbysingle-moleculelong-readmRNAsequencing.ProcNatlAcadSciUSA.2014;111(13):E1291–9. :10.1073/pnas. GordonS, gE,SalamovA,ZhangJ,MengX,etal.(2014)Widespreadpolycistronictranscriptsinmushroom-formingfungirevealedbysingle-moleculelong-readmRNAsequencing.BioRxiv, Hackl,T.,Hedrich,R.,Schultz,J.,Forster,F.(2014).proovread:large-scalehigh-accuracyPacBiocorrectionthroughiterativeshortreadconsensus.Bioinformatics,30(21),3004-3011.WuTD,WatanabeCK(2005)GMAP:agenomicmapandalignmentprogramformRNAandESTsequences.Bioinformatics21:1859–1875.PMID: FelipeA.Simão,RobertM.Waterhouse,PanagiotisIoannidis,EvgeniaV.Kriventseva,andEvgenyM.Zdobnov.BUSCO:assessinggenomeassemblyandannotationcompletenesswithsingle-copyorthologs.Bioinformatics31.19(2015):3210-FoissacS,SammethM.ASTALAVISTA:dynamicandflexibleysisofalternativesplicingeventsincustomgenedatasets.NucleicAcidsResearch2007,35(WebServerissue):W297-9.AbdelghanySE,HamiltonM,JacobiJL,etal.Asurveyofthesorghumtranscriptomeusingsingle-moleculelongreads[J].NatureCommunications,2016,7:11706.HaasBJ,PapanicolaouA,YassourM,GrabherrM,BloodPD,BowdenJ,etal.DenovotranscriptsequencereconstructionfromRNA-sequsingtheTrinityplatformforreferencegenerationandysis.NatProtoc.2013;8:1494–1512.:L.Kong,Y.Zhang,Z.Q.Ye,X.Q.Liu,S.Q.Zhao,L.Wei,andG.Gao.2007.CPC:assesstheprotein-codingpotentialoftranscriptsusingsequencefeaturesandsupportvectormachine.NucleicAcidsRes36:W345-349.LiangSun,HaitaoLuo,DechaoBu,GuoguangZhao,KuntaoYu,ChanghaiZhang,YuanningLiu,RunShengChenandYiZhao*Utilizingsequenceintrinsiccompositiontoclassifyprotein-codingandlongnon-codingtranscripts.NucleicAcidsResearch(2013),:10.1093/nar/gkt646.WangL,ParkHJ,DasariS,WangS,KocherJP,LiW.CPAT:Coding-PotentialAssessmentToolusinganalignment-freelogisticregressionmodel.NucleicAcidsRes.2013Apr1;41(6):e74.:10.1093/nar/gkt006.LiJ,MaW,ZengP,etal.LncTar:atoolforpredictingtheRNAtargetsoflongnoncodingRNAs[J].BriefingsinBioinformatics,2015,16(5):806.ZhengY,JiaoC,SunH,RosliHG,PomboMA,ZhangP,BanfM,DaiX,MartinGB,GiovannoniJJ,ZhaoPX,RheeSY,FeiZ(2016)iTAK:aprogramforgenome-widepredictionandclassificationofplanttranscriptionfactors,transcriptionalregulators,andproteinkinases.MolecularPlant9:1667-1670.Hong-MeiZhang,TengLiu,Chun-JieLiu,ShuangyangSong,XiantongZhang,WeiLiu,HaiboJia,YuXue,andAn-YuanGuo.AnimalTFDB2.0:aresourceforexpression,predictionandfunctionalstudyofanimaltranscriptionfactors.Nucl.AcidsRes.(28January2015)43(D1):D76-D81.AltschulSF,MaddenTL,Sch?fferAA,ZhangJ,etal.GappedBLASTandPSIBLAST:ANewGenerationofProteinDatabaseSearchPrograms.NucleicAcidsResearch1997,25(17):3389-3402.DengYY,LiJQ,WuSF,ZhuYP,etal.IntegratedNRDatabaseinProteinAnnotationSystemandItsLocalization.ComputerEngineering2006.,32(5):71-74.ApR,BairochA,WuCH,BarkerWC,etal.UniProt:theUniversalProteinknowledgebase.NucleicAcidsResearch2004Jan1;32(Databaseissue):D115-9.AshburnerM,BallCA,BlakeJA,BotsteinD,etal.Geneontology:toolfortheunificationofbiology.NatureGenetics2000,25(1):25-29.TatusovRL,GalperinMY,NataleDA.TheCOGdatabase:atoolforgenomescaleysisofproteinfunctionsandevolution.NucleicAcidsResearch2000,28(1):33-36.KooninEV,FedorovaND,JacksonJD,etal.Acomprehensiveevolutionaryclassificationofproteinsencodedcompleteeukaryoticgenomes.Genomebiology,2004,5(2):FinnRD,BatemanA,ClementsJ,etal.Pfam:theproteinfamiliesdatabase.NucleicAcidsResearch,2013:KanehisaM,GotoS,KawashimaS,OkunoY,etal.TheKEGGresourcefordecipheringthegenome.NucleicAcidsResearch2004,32(Databaseissue):D277-D280.MckennaA,HannaM,BanksE,etal.TheGenomeysisToolkit:aMapReduceframeworkforyzingnext-generationDNAsequencingdata[J].GenomeResearch,2010,20(9):1297-303.CingolaniP,PlattsA,LeLW,etal.Aprogramforannotatingandpredictingtheeffectsofsinglenucleotidepolymorphisms,SnpEff[J].Fly,2012,6(2):80.TrapnellC,WilliamsBA,PerteaG,etal.TranscriptassemblyandficationbyRNA-Seqrevealsunannotatedtranscriptsandisoformswitchingduringcelldifferentiation.[J].NatureBiotechnology,2010,28(5):511-5.AndersS,HuberW.Differentialexpressionysisforsequencecountdata[J].GenomeBiology,2010,LengN,DawsonJA,ThomsonJA,etal.EBSeq:anempiricalBayeshierarchicalmodelforinferenceinRNA-seqexperiments[J].Bioinformatics,2013,29(8):1035.FranceschiniA,SzklarczykD,FrankildS,et