kj-dbanbj--anbj华北林业实验中心6个无参转录组数据分析结题报告1

1、华北林业实验中心6个文冠果无参转录组数据分析结题报告2017/03/16目录1项目信息基本思想信息分析流程样本信息2组装与评估2.4组装序列说明组装结果统计组装长度分布组装后评估...6注释结果3.1比对率分析 均一性分析 准确性评估 注释比例评估Rea ds利用率评估蛋白比对率评估O R F 功能注释G O 分析C O G 注释4表达量分析4 .1 表达量估计.24.1.3表达量分布统计数据饱和度分析实验样品聚类4 .2 差异表达分析.

2、..6功能分析差异表达差异表达差异表达差异表达差异表达差异表达结果统计火山图 韦恩图 聚类图 注释蛋白互作网络55.1 G O 功能分析.25.1.3差异表达差异GO 统计GO 富集分析GO 富集DAG图5.2 K E G G 通路分析变异分析SSR分析附录6788 .18 .28 .38 .4参考文献数据库介绍软件与方法说明结果目录1项目信息1.1 基本思想安诺转录组是基于Illumina平台，通过研究某个物种在特定组织或者特定时期下的所有mRNA，针对实际样品情况采用灵活的差异分析策略，可以找到生物体不同时期、不同组织或不同间mRNA的表达差异，再通过软件

3、进行功能注释，最终可以得到mRNA参与生命活动的一个清晰的生物信息图谱。1.2 信息分析流程Illumina HiSeq所得原始下机序列，通过去除低质量序列、去接头等过程完成数据处理，得到高质量序列。对质控后的高质量序列进行组装，再基于组装结果，进行ORF分析及功能等后续信息分析，同时对不同处理组进行差异表达分析。流程图如下：图1 信息分析流程图如项目仅有一个样品，无法进行虚线所示的分析内容。1.3 样本信息项目共有6个样本用于分析，样品信息示例如下：表1 样品信息（1）Sample：结果中使用的样本名称；（2） Sample Description：客户提供的原始样本描述信息；（3） Gro

4、up：客户提供的原始样本分组信息。2组装与评估2.1 组装序列说明转录组无参分析中采用的组装软件为Trinity（20140717版）(Manfred,etal., 2011)。Trinity是由Broad和Hebrew University of Jerusalem共同研发而成，能够高效的构建de novo转录组。Trinity利用de Bruijn图论的基本原理，针对转录本具有可变剪接的特点，组装出全长转录本，其原理如下图：图2 Trinity组装原理S am p l eF l o we r bud1Sample Description.GroupFlowerbudTrinity组装原理概

5、述：（a）Inchworm基于k-mer图（中）对Read集合（上，短粗线）进行组 装，形成线性Contigs集合（下，彩线）；（b）Chrysalis对Inchworm形成的Contigs进行聚类， 若Contig（彩线）之间存在至少k-1-mer的重复，则形成一个Components，对每个Components构建一个de Bruijn图；（c）Butterfly对Chrysalis构建的de Bruijn图进行修剪、压缩、提取，得到同源的和（或）可变剪接形式的转录本。组装出来的转录本示例如下图：图3 FASTA格式说明图FASTA格式首先以大于号“>”开头，接着是序列的标识符，然后

6、是序列的描述信息。换行后是序列信息，序列中允许空格、换行和空行，直到下一个大于号，表示该序列的结束。基于过滤后的Clean Data，用Trinity组装出全长转录本序列，并基于转录本序列，取每个中最长的转录本序列作为Unigene。2.2 组装结果统计对Trinity组装出来的PotentialTranscript、Unigene进行各项指标的统计。结果如下表：表2 组装结果统计表（1）Count：总的组装出来的潜在转录本、潜在的个数；（2） Percent GC (%)：潜在转录本、潜在（3） Total Bases：潜在转录本、潜在的GC含量；的碱基数。He ade rT r i ni

7、tyU ni g e neCount127,325106,550Percent GC (% )41.010040.9900Total Bases84,780,02363,045,7302.3 组装长度分布对组装出来的Potential Transcript、Unigene统计其长度分布。其中N50代表序列从大到小排列，当其累计长度刚刚超过全部序列总长度50% 时，最后一个序列的大小即为N50的大小；N90代表序列从大到小排列，当其累计长度刚刚超过全部序列总长度90% 时，最后一个序列的大小即为N90的大小。结果如下表：表3 序列长度统计表（1）N50：Potential Transcript、

8、Unigene N50的结果；（2） N90：Potential Transcript、Unigene N90的结果；（3） Min：最短Potential Transcript、Unigene的长度；（4） Max：最长Potential Transcript、Unigene的长度；（5） Count：所有Potential Transcript、Unigene的个数；（6） Mean：Potential Transcript、Unigene的平均长度。B asi c S tatT r i ni tyU ni g e neN501,2561,046N90245230Min201201Max2

9、1,29621,296Count127,325106,550Mean665.8552591.7008表4序列长度分布表（1）200:400)：长度大于等于200，小于400的PotentialTranscript的个数，以下相同；（2）4000:-)：长度大于等于4000的Potential Transcript的个数。对Potential Transcript、Unigene统计其长度分布，分布图如下：图4 序列长度分布图Le ng th Hi st D i str i buti o nT r i ni tyU ni g e ne200:400)75,38369,004400:600)15,

10、59012,617600:800)7,7255,611800:1000)4,9263,3691000:1200)3,8502,5521200:1400)3,3412,2421400:1600)2,9171,9601600:1800)2,6171,7611800:2000)2,2111,5102000:2200)1,6801,1382200:2400)1,3909422400:2600)1,0607322600:2800)8835952800:3000)7004763000:3200)5493673200:3400)4803253400:3600)3842523600:3800)28920038

11、00:4000)2681834000:-)1,082714图中横坐标为序列长度分布，纵坐标为序列频数。对组装出来的转录本、Unigene统计每条序列的GC含量和长度，并且做成如下的图：图5 序列长度对GC含量分布图横坐标为序列长度，纵坐标为序列GC含量，图中每个点代表一条序列。一般而言，随着序列长度的增加，GC含量的波动逐渐减小。2.4 组装后评估2.4.1 . 1 比对率分析Bowtie2是一个超级快速的，较为节省内存的短序列拼接至模板组的工具。利用Bowtie2（2.2.3版）(Langmead, et al., 2009)将用于组装的序列与组装后的转录本序列进行比对，使其到组装转录本，然

12、后统计比对上序列的Reads比例。比对结果见下表：表5 比对结果统计表（1）Total Reads：过滤后Clean Reads总数；（2）Aligned Reads：能比对上转录本序列的Reads条数和百分比；Re ad_typ eC o untP e r ce nt( % )Total Reads318,119,840100Aligned Reads299,549,94594.1600Pair Mapping290,749,88891.4000（3）Pair Mapping：Read1和Read2均能比对上转录本序列的Reads条数和百分比。根据比对结果，可作成如下饼图：图6 比对结果统计

13、图2 . 4 . 2 均一性分析均一性是指的核酸序列的随机程度。若序列不偏向于的特定区域，则称其均一性好。若结果均一性很差，将直接影响转录组的各项分析结果。一般而言，的起始或终止区域，覆盖度会偏低，但在中间区域应该为一条总体趋势平稳的曲线。我们通过统计每个潜在转录本不同位置的深度来衡量测序结果均一性。由于不同转录本长度不同，将潜在转录本平均划分为100个窗口，统计每个窗口里面的碱基的深度，并计算所有潜在转录本在该对应窗口的深度平均值，依此作出如下序列均一性分布图：图7 均一性分析图将潜在转录本平均划分为100个窗口，统计每个窗口里面的碱基的深度。以潜在转录本的划分窗口为横坐标，以落在该窗口内的

14、碱基深度为纵坐标，做折线图。2 . 4 . 3 准确性评估对于一条组装序列来说，若其序列中含有较多的终止子，则这条序列的编码能力较弱。基于组装后的序列，随机选取1000条转录本，采用BLAST将其与近缘物种的蛋白序列进行比对，评估组装序列中由于框移错误导致的缺口以及过早终止（含终止子）的比例；若无近缘物种，则对于植物，默认比对到模式生物拟南芥的参考序列；对于动物，默认参考物种为人类。其组装准确性统计表如下：表6 准确性统计结果（1）Stop Codon (%)：由于框移错误导致的提前终止的转录本、Unigene序列的比例；（2）Gaps (%)：由于框移错误导致的序列缺口的转录本、Unigen

15、e序列的比例。2 . 4 . 4 注释比例评估Asse m bl i e sS to p C o do n( % )G ap s( % )Trinity2.06600.0000Unigene3.33300.0000对于组装结果，将其与近缘物种的cDNA序列进行比对。通常近缘物种cDNA序列的完整性以及物种之间的近缘程度会对注释结果有一定的影响。基于组装后的序列，随机选取1000条转录本，采用BLAT将其与近缘物种的cDNA序列进行比对，统计转录本序列的可注释比例。若无近缘物种，对于植物，默认比对到模式生物拟南芥的参考序列；对于动物，则默认参考物种为人类。其注释比例统计表如下：表7 注释比例统计

16、表（1） Total Matches (%)：比对上参考序列的转录本、Unigene的比例；（2） Matched Bases (%)：比对上参考序列的转录本、Unigene的碱基比对率的平均值；（3） Avg Accuracy (%)：比对上参考序列的转录本、Unigene的比对一致性的平均值。2 . 4 . 5 Rea ds利用率评估嵌合体指一个来源于1个或多个编码序列的一部分，组成新的。嵌合体的形成可能是反转录时转录本位置移动或异常重组等一系列原因导致的。随机选取100000条过滤后的Clean Reads，采用Bowtie将其比对到组装后序列，统计比对信息。其比对结果见下表：表8 Re

17、ads利用率及嵌合体比例统计表（1） NumReadsMatched (%)：100000条Clean Reads中比对到组装序列的Reads所占的比例；（2） AvgDist (bp)：双端均比对到组装序列的Read1与Read2之间的平均距离；（3）Chimera (%)：Read1与Read2比对到不同转录本的比例，即嵌合体比例。Asse m bl i e sT r i ni tyNumReadsMatched(% )90.0400AvgDist(Bp)59.5000Chimera(% )0.0400Asse m bl i e sT o talM atche dAvgM atche s(

18、 % )B ase s( % )Accur acy( % )Trinity24.400039.299372.5057Unigene20.600039.620373.11172 . 4 . 6蛋白比对率评估真核生物中存在一些高度保守的区域，利用这些蛋白的高度保守性，可以评估组装结果的完整性与可注释性。采用BLAST将组装后序列与2748个蛋白（真核生物的保守区域）进行比对，并统计比对率，其比对结果如下： 表9 蛋白统计表（1）Cegma Hit：比对到的蛋白的转录本、Unigene数目。3注释结果3.1开放阅读框（Open Reading Frame，ORF）是结构的正常核苷酸序列，从起始子到终

19、止子的阅读框可编码完整的多肽链，其间不存在使翻译中断的终止子。在无条件限制的情况下，DNA序列可以按六种框架阅读和翻译（每条链三种，对应三种不同的起始位点）。ORF识别包括检测这六个阅读框架并决定哪一个包含以启动子和终止子为界限的DNA序列而其内部不包含启动子或终止子，符合这些条件的序列有可能对应一个真正的单一的产物。ORF的识别是证明一个新的DNA序列为特定的蛋白质编码的部分或全部的先决条件。TransDecoder是一款专门针对RNASequencing的组装序列开放阅读框（ORF）的工具。对Unigene利用TransDecoder（20140717）鉴定编码区域，其中鉴定的标准有：1)

20、ORF的长度需要大于一定的长度；2)该序列的似然函数的log值需要大于0；3)选取6个ORF中分数最大的一个；Asse m bl i e sC e g m a Hi tTrinity2,667Unigene2,6634) 如果一个ORF完全包含另一个，那么输出最长的一个。Unigene.fasta.transdecoder.cds为出来的ORF序列，序列格式为FASTA格式。Unigene.fasta.transdecoder.pep为格式为FASTA格式。出来的蛋白序列，序列对出来的ORF序列统计其长度分布，结果如下图：图8 O RF长度分布图3.2 O RF 功能注释Trinotate是一

21、款全面的针对转录组组装序列的注释工具。它采用多种不同的方法，针对7大数据库来进行功能注释，其中包括同源性搜索、蛋白质结构域鉴定（HMMER/PFAM）、蛋白质信号（SingalP/TmHMM）以及与现存的注释数据库的比对（Uniprot/eggNOG/GO/KEGG等）。同源性搜索是根据已知的序列通过序列比对来的同源性；蛋白质结构域是指生物大中具有特异结构和功能的区域；信号肽指新多肽链中用于指导蛋白质的跨膜转移（定位）的N-末端的氨基酸序列（有时不一定在N端）。同源性搜索、蛋白质结构域鉴定与蛋白质信号肽的等均有助于理解潜在转录本的结构和功能。针对出来的ORF，采用Trinotate（20140

22、717）对其进行功能注释。注释结果示例如下：表10注释结果示例表（1）Gene_ID：组装出来的名；（2） NT_Seq-ID：（3） NT_Score：同NT数据库的最优比对结果；同NT数据库的比对得分；G e ne _I Dc52 9 4 4 _g 1NT_Seq-IDgi|568822394|ref|XM_006465553.1|NT_Score732NT_Evalue0.0000NT_DescriptionPREDICTED: Citrus sinensis protein disulfideis.NR_Seq-IDgi|568822395|ref|XP_006465616.1|NR_

23、Score2,047NR_Evalue0.0000NR_DescriptionPREDICTED: protein disulfide isomerase-like 1-5.BLASTX_IDsp|Q66GQ3|PDI16_ARATHBLASTX_Score677BLASTX_EvalueE:0BLASTX_DescriptionRecName: Full=Protein disulfide isomerase-like.BLASTP_IDsp|Q66GQ3|PDI16_ARATHBLASTP_Score678BLASTP_EvalueE:0BLASTP_DescriptionRecName:

24、 Full=Protein disulfide isomerase-like.PFAM_IDPF00085.15PFAM_NameThioredoxinPFAM_DescriptionThioredoxineggNOG_IDCOG0526eggNOG_DescriptionThiol-disulfide isomerase and thioredoxinsRNAMMER.Prot_IDc52944_g1|m.22258Prot_Coords274-1872+SignalPsigP:1270.819YESTmHMMExpAA=17.78PredHel=1Topology=i5-27oGO_bio

25、logical_processGO:0045454|cell redox homeostasisGO_molecular_functionGO:0003756|protein disulfide isomerase activityGO_cellular_componentGO:0005783|endoplasmic reticulum;GO:0005788|end.KOK09580|PDIA1, P4HB; protein disulfide- isomerase.Mapmap04141|Protein processing in endoplasmic reti.（4）NT_Evalue：

26、同NT数据库的比对Evalue值；（5）NT_Description：NT数据库中该的功能描述；（6） NR_Seq-ID：（7） NR_Score：同NR数据库的最优比对结果；同NR数据库的比对得分；（8）NR_Evalue：同NR数据库的比对Evalue值；（9）NR_Description：NR数据库中该的功能描述；（10） BLASTX_ID：（11） BLASTX_Score：同Uniprot数据库的最优比对结果；同Uniprot数据库的比对得分；（12）BLASTX_Evalue：同Uniprot数据库的比对Evalue值；（13）BLASTX_Description：Unipro

27、t数据库中该的功能描述；（14）BLASTP_ID：ORF同Uniprot数据库的最优比对结果；（15） BLASTP_Score：ORF同Uniprot数据库的比对得分；（16） BLASTP_Evalue：ORF同Uniprot数据库的比对Evalue值；（17） BLASTP_Description：Uniprot数据库中该ORF的功能描述；（18） PFAM_ID：PFAM数据库的最优比对结果；（19） PFAM_Name：PFAM数据库的蛋白名称；（20） PFAM_Description：PFAM数据库的功能描述；（21） eggNOG_ID：eggNOG数据库的最优比对结果；（2

28、2）eggNOG_Description：eggNOG数据库的功能描述；（23）RNAMMER：RNAMMERrRNA的结果；（24）Prot_ID：的ORF的ID；（25） Prot_Coords：（26） SignalP：SignalP的ORF在转录本上的位置；信号肽的结果；（27）TmHMM：TmHMM跨膜区的结果；（28）GO_biological_process：注释到的描述生物进程的GO Term；（29）GO_molecular_function：注释到的描述功能的GO Term；（30）GOponent：注释到的描述细胞组分的GO Term；（31） KO：注释到的KEGG中的

29、ID；（32） Map：注释到的通路的功能描述。为了更直观的了解注释结果，对其中注释到Uniprot、NR、NT的进行统计，作出如下韦恩图：图9 注释结果比较韦恩图图中NT、NR、BLASTX、BLASTP分别对应注释结果示例表中的4列，数字代表注释到此数据库的基因个数。3.3 G O 分析Gene Ontology（ 简称 GO,）是功能国际标准分类体系。作为本体会(Gene Onotology Consortium)所建立的数据库，它旨在建立一个适用于各种物种的，对和蛋白质功能进行限定和描述的，并能随着研究不断深入而更新的语言词汇标准。GO分为功能（MolecularFunction）、生

30、物过程（biologicalprocess）、和细胞组成（cellularcomponent）三个部分。或蛋白质可以通过ID对应或者序列注释的方法找到与之对应的GO编号，而GO编号可用于对应到Term，即功能类别或者细胞。利用Blast2GO，计算每个Term的数目。并根据GO统计结果绘制柱状图，结果如下：图10 GO 统计柱状图3.4 CO G 注释针对未知，采用Tatusov R L等提出的COG（Cluster of OrthologousGroup直系同源簇）方法(Tatusov, et al., 2000)，即用不同种族的成对相似聚类法把它们划分成各种直系同源簇，从而可以用同一簇中的

31、已知注释未知的功能。该方法通过对未知进行组水平上找寻直系同源体，从而借助相关直系同源体中已知未知“开放阅读框架”的生物学功能。1) 用BLAST软件默认参数将Unigene与COG数据库比对；2 ) 对数据库比对结果过滤，去掉e_value值大于1e-5的比对序列，同时选取同一序列比对多次中比对最好的一条；3) 对2）的结果进行注释。COG注释结果示例如下表：表11 CO G注释结果示例表（1） Gene_ID：组装出来的名；（2） Identity：一致性分值；（3） Eval：e_value值；（4） COG_gene：比对上的COG数据库中的（5） COG_num：比对上的COG数据库中

32、的名；ID；（6） COG_anno：比对上的COG数据库中的（7） Class：COG数据库对此功能的分类；（8） Class_anno：所属类的功能描述。描述；根据COG的注释结果，统计每类下的COG的个数，并做出如下的柱状图：图11 CO G统计图横坐标为COG数据库中的类别，纵坐标为匹配到的数。G e ne _I Dc59 6 77_g 4Identity32.6400Eval4e-26COG_geneSPAC167.08COG_numCOG2801COG_annoTransposase and inactivated derivativesClassLClass_annoReplic

33、ation, recombination and repair ;4表达量分析4 .1 表达量估计表达水平一般是通过该转录的mRNA的多少来衡量的。每个转录产生的mRNA的量，是受到时空等多种因素调控的，在不同的生长发育阶段，或者不同的组织水平，转录出mRNA的量都是不一样的。RPKM(Mortazavi, et al., 2008)是利用RNA-Seq技术用来定量估计表达值的有效工具，能够消除长度和量差异对计算表达的影响，计算得到的表达量可直接用于比较不同样品间的表达差异。RPKM即ReadsPer Kilobase Millon Mapped Reads，意为“一百万条Reads比对到长度

34、为1K的长度上的Reads数”。其计算公式为：图12 RPKM公式设RPKM(A)为A的表达量，则R为唯一比对到A的Reads数，N为唯一比对到的总Reads数，L为A的长度。4 . 1 . 1 表达量分布统计根据所有样品的表达量，得到该样品的表达量密度图。一般而言，差异表达的数量只占整体的小部分，因此所有样品应该具有类似的表达量分布情况。图13 表达量分布图对每组样品的表达量作密度分布图。横坐标为log2(RPKM+1)，纵坐标为的密度。不同颜色代表不同样品。根据每个样品的表达量，对每个样品进行绘制箱子图，查看样品的表达量整体分布趋势，得到所有样品的表达量的分布箱式图如下：图14 表达量箱式

35、分布图4 . 1 . 2 数据饱和度分析定量饱和曲线检查反映了表达水平定量对数据量的要求。表达量越高的，就越容易被准确定量；反之，表达量低的，需要较大的数据量才能被准确定量。根据每条饱和曲线达到平台期的数据量深度，可以推断本次实验的数据量是否满足数据分析需求。将序列按照10% ，20% .100% 随机取样（Subsample），并统计在每个数据量下获得的的表达量值，同用全部序列得到的表达量值进行比较，如果两个值的差异小于15% ，就认为该在该数据量下定量准确。根据的表达量不同，将分为四个等级，统计不同数据量下定量的准确性。图15 饱和度分析图横坐标为不同深度的数据量，纵坐标为在该深度下准确估

36、计表达量的百分比。4 . 1 . 3 实验样品聚类根据样品全部的表达量信息对样品进行系统聚类，同一实验条件下的样品通常会聚到一类，表明实验条件为影响聚类的主要因素。图16 全部样品的cluster聚类图根据每个样品的表达量，剔除掉在所有样品中表达量都小于1的进行聚类。计算两两之间的皮尔逊相关系数（Pearson Correlation Efficiency），来表示样品两两间的相似度。再利用系统聚类法（Hierarchical Cluster）将相似度高的样品归为一类，以此类推，最终得到样品的整体聚类结果。通过对所有表达进行聚类，可以找到某一类相似表达模式的，同时能展示一类在不同处理中的变化情

37、况。采用R中的kmean方法对所有进行kmean聚类分析，下图为所有的kmean聚类图：图17kmeans聚类图根据在每个样品里的表达量，取以2为底的对数后，再利用系k-均值聚类（k-meanclustering），最终得到的聚类结果。图中x轴为不同的样品，y轴为表达量的log2值，灰色的线表示在样品间表达量的变化，的线为该类的均值。4.2 2 差异表达分析4.2.1 . 1 差异表达结果统计针对不同的样品分组情况，采用不同的差异表达分析工具(Anders,etal., 2010)。无生物学重复，利用DEGseq软件(Wang, et al., 2010)，进行处理组与参考组的比较，并选取|l

38、og Ratio|1和q<0.05的2作为差异表达基因；针对存在生物学重复的情况，采用DEseq2进行处理组与参考组的比较，并选取|log Ratio|1和q<0.05的2作为差异表达。结果示例如下表： 表12 组间比较得到的差异表达示例表（1）Gene：组装得到的名；（2） *_count：样品*的Reads Count数；（3） *_normalize：每个组的标准化之后的结果；（4） FoldChange：组的标准化后的数值的倍数值；（5） Log2FoldChange：两组的标准化后的数值倍数的比例的log2值；（6）pval：计算的p值；（7） padj：校正之后的p值；

39、（8） Up/Down：上调还是下调表达，up上调，down为下调；G e nec8 32 0 5_g 1Flower_1_count8Flower_2_count0Flower_normalize1.2293Flowerbud1_count0Flowerbud2_count0Flowerbud_normalize0.1428FoldChange8.6086Log2FoldChange3.1057pval0.1822padjNAUp/DownupSignificantno（9）Significant：是否为显著性差异。根据组间比较得到的表达量差异倍数的分布散点图，结果如下图：图18 差异表达差

40、异倍数分布图横坐标表示，纵坐标表示log2foldchange值，黄色表示大于1的，表示小于-1的基因，灰色表示中间的，两条虚线分别表示倍数的对数值为1和-1。根据组间比较得到的差异，统计上下调个数，结果如下表： 表13 组间差异表达数目统计表（1） Up：表示在第一组中上调表达的（2） Down：表示在第二组中上调表达的（3） Total：表示两组间的差异表达； 个数。根据上表作差异表达个数统计图，如下图：N am eF l o we r _F l o we r budUp9,945Down5,871Total15,816图19 差异表达个数统计图4 . 2 . 2 差异表达火山图根据表达量

41、及其统计学显著性，作差异火山图，如下：图20 差异表达火山图横坐标为不同实验组中/不同样品中表达倍数变化，纵坐标为表达量变化的统计学显著程度，不同颜色表示不同的分类。4 . 2 . 3 差异表达韦恩图通过比较处理组和参考组，对不同组别之间的差异表达作韦恩图，如下：图21 差异韦恩图4 . 2 . 4 差异表达聚类图通过比较处理组和参考组，对差异表达进行聚类分析，可以很直观反映出差异表达的变化情况。我们利用R软件（3.1.1版），对差异表达和不同样本/实验条件同时进行分层聚类分析。下图为两组样本的差异表达聚类示例。图22 差异聚类图根据差异表达在每个样品里的表达量，取以2为底的对数后，计算欧氏距

42、离，再利用系统聚类法（Hierarchical Cluster），最终得到样品的整体聚类结果。在图中，表达量的变化用颜色的变化表示，表示表达量较低，黄色表示表达量较高。4 . 2 . 5 差异表达注释利用Uniprot、Pfam、GO和KEGG等数据库对差异表达进行注释，获得差异表达详细描述信息，结果示例如下表： 表14 差异表达的结果注释G e nec2 6 8 4 1 _g 1Flower_1_count90Flower_2_count42Flower_normalize108.2796Flowerbud1_count186Flowerbud2_count119Flowerbud_norm

43、alize266.7118FoldChange0.4059Log2FoldChange-1.3005pval0.0259padj0.0793Up/DowndownSignificantnoNT_Seq-IDgi|567900973|ref|XM_006442912.1|NT_Score1,372NT_Evalue0.0000NT_DescriptionCitrus clementina hypothetical protein (CICLE_v.NR_Seq-IDgi|568850053|ref|XP_006478743.1|NR_Score4,176NR_Evalue0.0000NR_Des

44、criptionPREDICTED: LRR receptor-like serine/threonine-p.BLASTX_IDsp|Q9FL28|FLS2_ARATHBLASTX_Score1,214BLASTX_EvalueE:0BLASTX_DescriptionRecName: Full=LRR receptor-like serine/threonin.BLASTP_IDsp|Q9FL28|FLS2_ARATHBLASTP_Score1,259BLASTP_EvalueE:0BLASTP_DescriptionRecName: Full=LRR receptor-like seri

45、ne/threonin.PFAM_IDPF08263.7PFAM_NameLRRNT_2PFAM_DescriptionLeucine rich repeaerminal domaineggNOG_IDCOG0515eggNOG_DescriptionSerine/threonine protein kinaseRNAMMER.Prot_IDc26841_g1|m.2451Prot_Coords191-3769-SignalPsigP:1400.504YES（1）Gene：组装得到的名；（2） *_count：样品*的Reads Count数；（3） *_normalize：各组表达量标准化后

46、的结果；（4） FoldChange：两组的标准化后的数值倍数的比例值；（5）Log2FoldChange：两组的标准化后的数值倍数的比例的log2值；（6） pval：计算的p值；（7） padj：校正之后的p值；（8） Up/Down：上调还是下调表达，Up上调，Down为下调；（9） Significant：是否为显著性差异；（10） NT_Seq-ID：（11） NT_Score：同NT数据库的最优比对结果；同NT数据库的比对得分；（12）NT_Evalue：同NT数据库的比对Evalue值；（13）NT_Description：NT数据库中该的功能描述；（14） NR_Seq-ID：

47、（15） NR_Score：（16） NR_Evalue：同NR数据库的最优比对结果； 同NR数据库的比对得分；同NR数据库的比对Evalue值；（17）NR_Description：NR数据库中该的功能描述；（18）BLASTX_ID：同Uniprot数据库的最优比对结果；（19） BLASTX_Score：（20） BLASTX_Evalue：同Uniprot数据库的比对得分；同Uniprot数据库的比对Evalue值；（21）BLASTX_Description：Uniprot数据库中该的功能描述；（22） BLASTP_ID：ORF同Uniprot数据库的最优比对结果；（23） BLA

48、STP_Score：ORF同Uniprot数据库的比对得分；（24） BLASTP_Evalue：ORF同Uniprot数据库的比对Evalue值；（25） BLASTP_Description：Uniprot数据库中该ORF的功能描述；（26） PFAM_ID：PFAM数据库的最优比对结果；（27） PFAM_Name：PFAM数据库的蛋白名称；（28） PFAM_Description：PFAM数据库的功能描述；（29） eggNOG_ID：eggNOG数据库的最优比对结果；（30） eggNOG_Description：eggNOG数据库的功能描述；（31）RNAMMER：RNAMMER

49、rRNA的结果；（32）Prot_ID：的ORF的ID；（33）Prot_Coords：的ORF在转录本上的位置；TmHMMExpAA=42.86PredHel=2Topology=i13-35o822-844iGO_biological_processGO:0052544|defense response by callosedepositi.GO_molecular_functionGO:0005524|ATP binding;GO:0004675|transmembrane.GO_cellular_componentGO:0005768|endosome;GO:0010008|endos

50、omemembran.KOK13420|FLS2; LRR receptor-like serine/threonine.Mapmap04626|Plant-pathogen interaction（34）SignalP：SignalP信号肽的结果；（35）TmHMM：TmHMM跨膜区的结果；（36）GO_biological_process：注释到的描述生物进程的GO Term；（37）GO_molecular_function：注释到的描述功能的GO Term；（38）GOponent：注释到的描述细胞组分的GO Term；（39） KO：注释到的KEGG中的ID；（40） Map：注释到

51、的通路的功能描述。4 . 2 . 6 差异表达蛋白互作网络应用STRING蛋白质互作数据库（）中的互作关系，针对数据库中包含的该物种的近缘物种，将组装后序列与近缘物种的蛋白序列进行比对。通过比对结果以及近缘物种的蛋白互作关系，将差异表达集（比如差异List）到近缘物种的蛋白互作网络。我们提供参考物种蛋白互作网络数据文件（蛋白互作网络参考数据集）和目标集，构建蛋白互作网络的网络数据文件可以直接导入Cytoscape软件(Shannon, et al. , 2003)，并根据目标集中的属性对网络进行可视化编辑。Cytoscape软件使用方法可参考我们提供的使用说明文档。您可以针对一些网络的拓扑属性

52、进行统计和标示作图，比如互作网络图中节点（Node）的颜色与此的上下调特性相关。若上调，则其节点为红色；若下调，则其节点为青色。根据不同的研究目的和需求，您还可以在网络图中进行调整节点位置和颜色、标注表达量水平等操作。按我们提供的使用说明将文件导入Cytoscape软件后的效果图如下：图23 差异蛋白互作网络图5功能分析5.1 G O 功能分析5.1.1 . 1 差异表达GO 统计如果研究物种具有相关GO注释数据库，直接以该数据库进行GO的分析；无相应GO注释数据库的可以采用Blast2GO得到每个对应的GO条目。针对GO数据库中第三层的条目，统计差异表达（区分上调表达和下调表达）在该条目里的个数（百分比），得到的结果如下：表15 GO 统计示例表（1） GO（2） GO（3） UpTerm：GO大类； Subterm：GO子类名称；Count：位于该子类的上调差异表达数目；（4）Up的比例；Percent：位于该子类的上调差异表达占注释到所有GOTerm的上调差异表达（5） Down Count：位于该子类的下调差异表达（6） Down Percent：位于该子类的下调差异表达因的比例。数目；占注释到所有GO Term的下调差异表达基为了直观的展示差异表达集合的GO统