专利名称：：一种基于EST数据库和UniGene数据库的基因挖掘方法
技术领域：
：本发明涉及应用生物信息学领域，具体是涉及一种挖掘UniGene数据库中特异表达基因的方法。
背景技术：
：EST(表达序列标签)测序是高通量检测基因表达信息的方法之一，近些年来，由于EST序列长度较长，特异性较好，比基因芯片的噪音低，使得EST测序得到了广泛的应用。但是，由于EST测序的成本较高，单个研究人员在建立EST序列信息时，对EST序列的挖掘工作，主要集中在EST序列信息的发现方面，例如，EST序列的拼接，拼接结果的注释以及SSR和SNP的发现等方面，往往EST测序的EST序列数目较少，即对EST序列的定量方面的研究较少，只能基于EST序列信息定性分析生物体的性状，无法通过EST序列定量分析基因和性状间的关系。但是，随着大规模EST测序计划的进行，公共数据库(EST数据库和UniGene数据库)中的EST序列积累越来越多，研究人员开始利用EST序列来定量分析基因和性状间的关系，例如，Aouacheria.A.等人在2006年第7期的BMCGenomics上，发表了名为"BioinformaticscreeningofhumanESTsfordifferentiallyexpressedgenesinnormalandtumortissues(利用人类EST的生物信息学方法筛选正常组织和肿瘤组织中差异表达基因)"的文章，阐述了采用BLAST搜索的方法对EST表达量进行计数统计，定量分析了基因和癌症的关系，挖掘出了人类癌症响应基因。但是，由于旁系同源基因和一些高相似性的保守结构域序列的干扰，使得该方法计数会出现假阳性，导致EST表达量计数错误，从而使得后续的统计分析产生较大的误差，无法正确挖掘与性状相关的诱导基因。
发明内容为了克服现有的基因挖掘方法中EST表达量计数出现假阳性进而导致无法正确挖掘与性状相关的诱导基因的问题，本发明提供一种基于EST数据库和UniGene数据库的基因挖掘方法
技术领域：
：本发明的一种基于EST数据库和UniGene数据库的基因挖掘方法的具体过程为步骤A:下载EST数据库和UniGene数据库，并对所述的UniGene数据库中的EST序列信息按照物种类型进行分类，接下来执行步骤B;步骤B:对所述的EST文库注释信息进行信息检索，进而对EST文库注释信息进行分类，接下来执行步骤C;步骤C:根据EST文库注释信息的分类信息计算表达基因UniGene转录组的EST表达量，接下来执行步骤D;步骤D:对所获得的表达基因UniGene转录组的EST表达量进行超几何分布检验，计算获得表达基因UniGene转录组差异表达的超几何分布检验值Pialue，接下来执行步骤E;步骤E:采用FDR方法调整表达基因UniGene转录组差异表达的超几何分布检验值Pialue，接下来执行步骤F;步骤F:设置表达基因UniGene转录组差异表达的超几何分布检验值Pialue阈值为O.01，筛选异常状态响应基因，接下来执行步骤G;步骤G:利用RT-PCR技术验证所筛选出的异常状态响应基因为与性状相关的诱导基因。本发明综合EST数据库和UniGene数据库的数据，避免了拼接或是比对过程中人为引入的误差；本发明还完成了对EST数据的定量分析，而不是简单的定性的分析，从表达量水平上揭示了异常状态所响应的基因本质；本发明还采用超几何分布检验方法计算差异表达的超几何分布检验值P-value,结合FDR方法调整所计算的超几何分布检验值P-value,避免了多次检验引入的误差，使得基因挖掘的结果更加精确。本发明可以应用于人类疾病发生、动植物生长发育调控、动植物疾病调控过程以及动植物逆境胁迫等性状相关诱导基因的挖掘，本发明从公共数据库(EST和UniGene)中挖掘得到了目的性状的响应基因，高效、准确地解析了生物性状，为生物重要过程的揭示奠定了重要的基础。图l是本发明的一种基于EST数据库和UniGene数据库的基因挖掘方法的工作流程图，图2是大豆逆境胁迫诱导基因的RT-PCR半定量分析结果示意图。具体实施例方式具体实施方式一参见图l，本具体实施方式所述的一种基于EST数据库和UniGene数据库的基因挖掘方法的具体过程为步骤A:下载EST数据库和UniGene数据库，并对所述的UniGene数据库中的EST序列信息按照物种类型进行分类，接下来执行步骤B;步骤B:对所述的EST文库注释信息进行信息检索，进而对EST文库注释信息进行分类，接下来执行步骤C;步骤C:根据EST文库注释信息的分类信息计算表达基因UniGene转录组的EST表达量，接下来执行步骤D;10步骤D:对所获得的表达基因UniGene转录组的EST表达量进行超几何分布检验，计算获得表达基因UniGene转录组差异表达的超几何分布检验值Pialue，接下来执行步骤E;步骤E:采用FDR方法调整表达基因UniGene转录组差异表达的超几何分布检验值Pialue，接下来执行步骤F;步骤F:设置表达基因UniGene转录组差异表达的超几何分布检验值Pialue阈值为O.01，筛选异常状态响应基因；步骤G:利用RT-PCR技术验证所筛选出的异常状态响应基因为与性状相关的诱导基因。本具体实施方式中的EST数据库含有EST序列，本具体实施方式中的UniGene数据库是NCBI开发的，是通过蛋白质的相似性、基因表达信息、cDNA克隆和基因组位置等信息，将同一转录位点的一组转录序列拼接成的一致性序列。UniGene数据库包含两部分信息即EST序列信息及EST文库注释信息，所述信息具有较高可信度，是真实存在的转录本信息，适用于性状相关诱导基因的挖掘。本具体实施方式中所述步骤B的具体过程为首先，收集生物体异常状态的关键词，将所收集的每一个生物体异常状态的关键词在所述的EST文库注释信息中进行信息检索，筛选异常状态EST文库，并提取异常状态EST文库的ID，然后，收集生物体正常状态的关键词，将所收集的每一个生物体正常状态的关键词在所述的EST文库注释信息中进行信息检索，筛选正常状态EST文库，并提取正常状态EST文库的ID，所述的在EST文库注释信息中进行信息检索的检索项包括主题TITLE、发育状态DEVELOPMENTAL—STAGE和组织VERBAT頂—TISSUE三项。本具体实施方式中所述步骤C的具体过程为在步骤C中，根据EST文库注释信息的分类信息，从所述UniGene数据库中的UniGene转录组文件中提取表达基因UniGene转录组的正常状态EST序列信息的条数和异常状态EST序列信息的条数进行计数，同时对UniGene转录组文件中所有的UniGene转录组的正常状态EST序列信息的条数和UniGene转录组文件中所有的UniGene转录组的异常状态EST序列信息的条数进行计数，将上述所获得的EST序列信息的所有计数信息转化为EST表达量；如果某生物体的UniGene数据库中已经对每个UniGene转录组按照EST文库注释信息建立完善的分类信息，则直接提取EST序列信息的所有计数信息，将所述计数信息转化为相应的EST表达量。本具体实施方式中所述步骤D的具体过程为设表达基因UniGene转录组的正常状态EST表达量为a，表达基因UniGene转录组的异常状态EST表达量为b，且UniGene转录组文件中所有UniGene转录组的总的正常状态EST表达量为c，所述UniGene转录组文件中所有UniGene转录组的总的异常状态EST表达量为d，以a、b、c和d四个数构建超几何分布检验，利用免费开源软件R平台的phyper完成超几何分布检验，公式一Y。"'daJ利用公式一计算获得表达基因UniGene转录组的超几何分布检验值Pialue。本具体实施方式中步骤E的具体过程为将UniGene转录组文件中所有UniGene转录组差异表达的超几何分布检验值Pialue按照从小到大的顺序依次排列，将p(i)定义为排列后的第i个超几何分布检验值P-value,同时设所述UniGene转录组文件中所有UniGene转录组的总数为n，所述的i和n均为自然数，且i《n，则调整后的第i个超几何分布检验值Pialue为P-Value=p(i)*n/i，此等式中所述的P-value和p(i)是重新赋值关系。本具体实施方式中步骤F的具体过程为将调整后的超几何分布检验值Pialue阈值以下的超几何分布检验值Pialue所对应的表达基因UniGene转录组作为异常状态的响应基因。在本具体实施方式中，以表达基因UniGene转录组的EST序列信息的条数作为表达基因UniGene转录组的EST表达量，依据表达基因UniGene转录组的EST文库描述信息，对所述EST文库进行分类，通过对比正常状态EST表达量和异常状态EST表达量构建超几何分布检验，进行统计分析，采用FDR方法进行基因差异表达的误差控制，筛选异常状态的响应基因。本具体实施方式适用于生物异常过程差异表达基因的挖掘，如人类疾病发生、动植物生长发育调控、动植物疾病调控过程以及动植物逆境胁迫等过程。具体实施方式二本实施方式所述的一种基于EST数据库和UniGene数据库的基因挖掘方法，与具体实施方式一的不同之处在于它是对人类癌症响应基因的挖掘在步骤A中，下载人类EST数据库和人类UniGene数据库；在步骤B中，直接提取人类EST文库注释信息，首先，提取人类异常状态EST文库，并提取所述的人类异常状态EST文库的ID，然后，提取人类正常状态EST文库，并提取所述的人类正常状态EST文库的ID;在步骤C中，根据人类EST文库注释信息的分类信息，从所述人类UniGene数据库中的人类UniGene转录组文件Hs.profiles中提取人类表达基因UniGene转录组的健康状态HealthState下孚L房月中瘤breasttumor，白血病leukemia，卵巢月中瘤ovariantumor，中枢神经系统CNS的原始神经夕卜胚层月中瘤primitiveneuroectodermaltumor，前歹[J腺癌prostatecancer和正常状态normal的人类EST序列信息的条数进行计数，同时对人类UniGene转录组文件Hs.profiles中所有的UniGene转录组的正常状态人类EST序列信息的条数和人类UniGene转录组文件Hs.profiles中所有的UniGene转录组的异常状态人类EST序列信息的条数进行计数，将上述所获得的人类EST序列信息的所有计数信息转化为人类表达基因UniGene转录组的EST表达量；在步骤D中，对所获得的人类表达基因UniGene转录组的EST表达量进行超几何分布检验，其中，以人类UniGene转录组文件Hs.profiles中正常状态normal下的数据作为人类UniGene转录组的正常状态EST表达量，以人类UniGene转录组文件Hs.profiles中其它癌症病变状态作为异常状态，并将所述异常状态下的数据作为人类表达基因UniGene转录组的异常状态EST表达量，利用所获得的人类表达基因UniGene转录组的正常状态EST表达量、人类表达基因UniGene转录组的异常状态EST表达量、人类UniGene转录组文件Hs.profiles中所有UniGene转录组的总的正常状态EST表达量和人类UniGene转录组文件Hs.profiles中所有UniGene转录组的总的异常状态EST表达量构建超几何分布检验，利用免费开源软件R平台的phyper完成超几何分布检验，计算获得人类表达基因UniGene转录组差异表达的超几何分布检验值Pialue;在步骤E中，采用FDR方法调整人类表达基因UniGene转录组差异表达的超几何分布检验值Pialue，该调整过程为将人类UniGene转录组文件Hs.profiles中所有UniGene转录组差异表达的超几何分布检验值Pialue按照从小到大的顺序依次排列，将p(i)定义为排列后的第i个超几何分布检验值Pialue，同时设人类UniGene转录组文件Hs.profiles中所有UniGene转录组的总数为n，所述的i和n均为自然数，且i《n，则调整后的第i个超几何分布检验值P-value为P-value=p(i)*n/i，此等式中所述的P-value和p(i)是重新赋值关系；在步骤F中，设置人类癌症所响应的表达基因UniGene转录组差异表达的超几何分布检验值Pialue阈值为O.Ol，筛选人类癌症响应基因，在此筛选过程中，将调整后的超几何分布检验值Pialue阈值以下的超几何分布检验值Pialue所对应的表达基因UniGene转录组作为人类癌症的响应基因；在步骤G中，利用RT-PCR技术验证所筛选出的人类癌症响应基因为与人类癌症相关的诱导基因。本实施方式提供了一种对人类癌症响应基因的挖掘方法，本实施方式中所述的人类UniGene转录组文件Hs.profiles内容为>2IHealthStateadrenaltumor□/12841breast(mammsrytumoi:cervicaltumor0/chondros&rc0/colorectaltumor□esophagealtumoi:□gei:mcelltumor0/glioma□/107493hesd&ndnecktuinor3458182865/114863/173552647250/13740Sinsulinoma□/30305kidneytumor□/69383leukeniis0/956271ivertumor1/96641lungtumor□/103480lymphoma0/72064norL-neoplasis□/97513noi:脆l16/3374366ovariantumoi:□/77210pancreaticcsncei:□/74633primitiveneuroectodermaltumoi:oftheCNS□prostatecancer0/103951retinolnlastoina0/46517skintumoi:□/124881sins11intestineadenocsrcinorna□/12684125405本实施方式中提取的Hs.profiles中健康状态HealthState下乳房肿瘤breasttumor，白血病leukemia，卵巢肿瘤ovariantumor，中枢神经系统CNS的原始神经外胚层肿瘤primitiveneuroectodermaltumor，前列腺癌prostatecancer禾口正常状态normal的EST序列信息的计数信息为Hs.2163374366Hs.46343374366Hs.ll113374366Hs.127337436601633743668337436603873374366Hs.6688337436612033743"18033743662Hs.1021593374366133374366Hs.109463374366Hs.Ill233374366Hs.120403374366Hs.129233743"□2053374366□94573□95627077210□126405□103951394573□077210□126405□103951□945731196627□772100126405010395194573196627□7721001264050103951094573095627094573249562701294573696627□945731956270194573209662729457349562743945733966270945730966271945733094573□96627094573□9662794573□0772100126405010395177210□126405□103951□77210491264051010395177210□126405□1039512772100126405010395177210131264051103951□772102126405010395127721001264050103951□772100126405o103951□772100126405o103951□772100126405110395177210□126405□1039511264057103951从人类UniGene转录组文件Hs.profiles中提取人类EST序列信息计数信息的方法可采用下述程序实现open(human—file,〃〈$ARGV〃)；while(defined($in—liiie=〈humaii—file>))14chomp$in—line;if($in—line=/"〉/&&$in—line=/HealthState/)$ugid=(split(/\s+/，$in—line))[1];$ugid=(split(AI/，$ugid));$ugid="Hs.$ugid";if($in—line=/breast\(mammarygland\)tumor/)錢temp二split(/\s+/，$in—line);@br=@temp[$frtemp—2，$frtemp];if($in—line=厂leukemia/)錢temp二split(/\s+/，$in—line);@leu=@temp[$frtemp—2，$frtemp];if($in—line=厂normal/)錢temp二split(/\s+/，$in—line);@nm=@temp[$frtemp—2，$frtemp];if($in—line=/ovariantumor/)錢temp二split(/\s+/，$in—line);@ov=@temp[$frtemp—2，$frtemp];if($in—line=/primitiveneuroectodermaltumoroftheCNS/)錢temp二split(/\s+/，$in—line);@cns=@temp[$frtemp-2，$frtemp];if($in—line=/prostatecancer/)@temp=split(/\s+/，$in—line);@pr=@temp[$frtemp-2，$frtemp];if($in—line=/BodySites/)print〃$ugid\t$nm\t$nm[l]\t$br\t$br[1]\t$leu\t$leu[l]\tprint〃$ov\t$ov[1]\t$cns\t$cns[1]\t$pr\t$pr[1]\t〃；close(human—file);本实施方式中对人类表达基因UniGene转录组差异表达的EST表达量进行超几何分布检验，计算获得的人类表达基因UniGene转录组在各种癌症状态下与在正常状态(normal)下差异表达的超几何分布检验值Pialue，所述超几何分布检验值Pialue具体为UGID朋LEovPR0s-2..11111Bs-41111Bs.111t13E—12111Bs-12111.1Bs-3611111Bs-4113_細—31111OL40594401OL9^304031p1i17E—120-74^512125Bs-6611110l9647277740lGOO4S31111Bs-960l65^43390-206558605Ol119573633u07.E—050.91*364717Bs_10^.pL3314743021p_"11333491Bs.10411OL04331600911Bs.I*OL727236779OL15571853711111111111Bs.12011110.7117GS821Bs.129111110l1331908276l35E~05Ol506409974&"E—05Ol44326299上述获得超几何分布检验值Pialue可通过下述程序实现data—all=read.table("human—est.txt〃，header=F，sep=〃\t〃)data=data—all[(data—all[，1]%in%microarray.pvalue[，l])，];16鼎鼎BR鼎鼎鼎鼎鼎est二cbind(data[，4]-1，data[，5]，data[，3]，data[，2]+data[，4]);row.names(est)=data[，1];pvalue二l:nrow(est);for(iin1:nrow(est))pvalue[i]二phyper(est[i，l]，est[i，2]，est[i，3]，est[i，4]，lower.tail二F);BR=pvalue;鼎鼎LE鼎鼎鼎鼎鼎est二cbind(data[，6]-1，data[，7]，data[，3]，data[，2]+data[，6]);row.names(est)=data[，1];pvalue二l:nrow(est);for(iin1:nrow(est))pvalue[i]二phyper(est[i，l]，est[i，2]，est[i，3]，est[i，4]，lower.tail二F);LE=pralue;鼎鼎OV鼎鼎鼎鼎鼎est二cbind(data[，8]-1，data[，9]，data[，3]，data[，2]+data[，8]);row.names(est)=data[，1];pvalue二l:nrow(est);for(iin1:nrow(est))pvalue[i]二phyper(est[i，l]，est[i，2]，est[i，3]，est[i，4]，lower.tail二F);0V=pralue;鼎鼎CNS鼎鼎鼎鼎鼎est二cbind(data[，IO]-I，data[，ll]，data[，3]，data[，2]+data[，10]);row.names(est)=data[，1];pvalue二l:nrow(est);for(iin1:nrow(est))pvalue[i]二phyper(est[i，l]，est[i，2]，est[i，3]，est[i，4]，lower.tail二F);CNS=pvalue;鼎鼎PR鼎鼎鼎鼎鼎est二cbind(data[，12]-1，data[，13]，data[，3]，data[，2]+data[，12]);row.names(est)=data[，1];pvalue二l:nrow(est);for(iin1:nrow(est))pvalue[i]二phyper(est[i，l]，est[i，2]，est[i，3]，est[i，4]，lower.tail二F);PR=pvalue;UGID=data[，1];est.pvalue=data.frame(UGID，BR，LE，0V，CNS，PR);write,table(est.pvalue，〃Hs—est—pvalue.txt〃，col.names=T，row.names=F，sep=〃\tquote二F);他.off=0.05;temp二sort(est.pvalue[，〃BR〃]);for(iin1:length(temp))if(temp[i]氺length(temp)/i〉=fdr.off){fdr.off二temp[i];break;};unigene.est二as.character(est.pvalue[est.pvalue[，〃BR〃]〈=fdr.off，1]);他.off=0.05;temp二sort(est.pvalue[，〃LE〃]);for(iin1:length(temp))if(temp[i]氺length(temp)/i〉=fdr.off){fdr.off二temp[i];break;};unigene.est二as.character(est.pvalue[est.pvalue[，〃LE〃]〈=fdr.off，1]);他.off=0.05;temp二sort(est.pvalue[，〃0V〃]);for(iin1:length(temp))if(temp[i]氺length(temp)/i〉=fdr.off){fdr.off二temp[i];break;};unigene.est二as.character(est.pvalue[est.pvalue[，〃0V〃]〈=fdr.off，1]);他.off=0.05;temp二sort(est.pvalue[，〃CNS〃]);for(iin1:length(temp))if(temp[i]氺length(temp)/i〉=fdr.off){fdr.off二temp[i];break;};unigene.est二as.character(est.pvalue[est.pvalue[，〃CNS〃]〈=fdr.off，1]);他.off=0.05;temp二sort(est.pvalue[，〃PR〃]);for(iin1:length(temp))if(temp[i]氺length(temp)/i〉=fdr.off){fdr.off二temp[i];break;};unigene.est二as.character(est.pvalue[est.pvalue[，〃PR〃]〈=fdr.off，1]);具体实施方式三本实施方式所述的一种基于EST数据库和UniGene数据库的基因挖掘方法，与具体实施方式一的不同之处在于它是对大豆逆境胁迫响应基因的挖掘在步骤A中，下载大豆EST数据库和大豆UniGene数据库；在步骤B中，以"cold"，"salt"和"drought"作为与大豆逆境胁迫相关的异常状态的关键词，将所述关键词在大豆EST文库注释信息文件Gma.lib.info中进行信息检索，筛选大豆异常状态EST文库，并提取大豆异常状态EST文库的ID，然后，将大豆EST文库中的其他EST文库作为大豆正常状态EST文库，并提取大豆正常状态EST文库的ID，所述的在大豆EST文库注释信息文件Gma.lib.info中进行信息检索的检索项包括主题TITLE、发育状态DEVELOPMENTAL—STAGE和组织VERBAT頂—TISSUE三项；在步骤C中，根据大豆EST文库注释信息的分类信息，从所述大豆UniGene数据库中的大豆UniGene转录组文件Gma.gb—cid—lid中提取大豆表达基因UniGene转录组的正常状态EST序列信息的条数和异常状态EST序列信息的条数进行计数，同时对大豆UniGene转录组文件Gma.gb—cid—lid中所有的UniGene转录组的正常状态大豆EST序列信息的条数和大豆UniGene转录组文件Gma.gb—cid—lid中所有的UniGene转录组的异常状态大豆EST序列信息的条数进行计数，将上述所获得的大豆EST序列信息的所有计数信息转化为大豆表达基因UniGene转录组的EST表达量；在步骤D中，对所获得的大豆表达基因UniGene转录组差异表达的EST表达量进行超几何分布检验，其中，以大豆UniGene转录组文件Gma.gb—cid—lid中正常状态normal下的数据作为大豆表达基因UniGene转录组的正常状态EST表达量，以大豆UniGene转录组文件Gma.gb—cid—lid中异常状态下的数据作为大豆表达基因UniGene转录组的异常状态EST表达量，利用所获得的大豆表达基因UniGene转录组的正常状态EST表达量、大豆表达基因UniGene转录组的异常状态EST表达量、大豆UniGene转录组文件Gma.gb—cid—lid中所有UniGene转录20组的总的正常状态EST表达量和大豆UniGene转录组文件Gma.gb—cid—1id中所有UniGene转录组的总的异常状态EST表达量构建超几何分布检验，利用免费开源软件R平台的phyper完成超几何分布检验，计算获得大豆表达基因UniGene转录组差异表达的超几何分布检验值P-value;在步骤E中，采用FDR方法调整大豆表达基因UniGene转录组差异表达的超几何分布检验值Pialue，该调整过程为将大豆UniGene转录组文件Gma.gb—cid—lid中所有UniGene转录组差异表达的超几何分布检验值Pialue按照从小到大的顺序依次排列，将p(i)定义为排列后的第i个超几何分布检验值Pialue，同时设大豆UniGene转录组文件Gma.gb—cid—lid中所有UniGene转录组的总数为n，所述的i和n均为自然数，且i《n，则调整后的第i个超几何分布检验值P-value为P-value=p(i)*n/i，此等式中所述的P-value和p(i)是重新赋值关系；在步骤F中，设置大豆逆境胁迫所响应的表达基因UniGene转录组的超几何分布检验值Pialue阈值为O.Ol，筛选大豆逆境胁迫响应基因，在此筛选过程中，将调整后的超几何分布检验值Pialue阈值以下的超几何分布检验值Pialue所对应的表达基因UniGene转录组作为大豆逆境胁迫的响应基因。在步骤G中，利用RT-PCR技术验证所筛选出的大豆逆境胁迫响应基因为与大豆逆境胁迫相关的诱导基因，具体验证过程为第一挑选步骤F中所述的大豆逆境胁迫响应基因进行RT-PCR验证具体挑选过程为将步骤F中筛选出的大豆逆境胁迫的响应基因按照超几何分布检验值Pialue从小到大的顺序排列，挑选出最小的8个超几何分布检验值Pialue所对应的大豆逆境胁迫的响应基因，所述8个超几何分布检验值Pialue所对应的大豆逆境胁迫的响应基因分别为Gma.2044、Gma.33428、Gma.32516、Gma.34982、Gma.8518、Gma.38、Gma.3774和Gma.22054;第二大豆总RNA提取以正常培养的大豆21日龄幼苗进行4。C低温cold、200mmolL—1NaCl盐salt和lOOmmolL—工甘露醇干旱drought胁迫处理l小时后的大豆21日龄幼苗叶片为材料，采用TRIZOL试剂法提取大豆总RNA，所述采用TRIZOL试剂法提取总RNA的方法参照试剂盒说明书进行，具体为i.称取100mg所述幼苗叶片样品放入高温高压消毒研钵，迅速在液氮中研磨成粉末，然后加入lmlTRIZOL研磨成匀浆组织；ii.将所述匀浆组织作为标本转移到1.5ml的离心管中，在153(TC环境下放置5min，21以彻底分离核蛋白复合体；iii.向所述离心管中加入0.2ml的氯仿，加盖好后用手剧烈摇晃15s，在1530。C环境下放置23min，然后离心12000rmin—、15min，28。C，离心后离心管内物质分成三层，下层为红色的酚，中间层为氯仿相，上层为无色的水相，RNA只存在于水相中，水相占总TRIZOL的60o/o;iv.将上层水相转移到另一个干净的离心管中，加入0.5ml异丙醇，在153(TC环境下静置IOmin，然后离心12000rmin—、10min，28。C;v.去步骤iv中所述离心管中上清，加入lml75。/。乙醇洗涤RNA沉淀，采用振荡器混匀，然后离心7500r.min—、5min，28°C;vi.继续去步骤v中离心管中上清，并将去上清后的离心管置于真空或空气中510min，干燥RNA沉淀，加DEPC处理的无菌水15yl，并将加入无菌水的离心管置于-2(TC环境下保存备用，提取总RNA完成；第三cDNA的合成以所提取的总RNA为模板，在下游引物的引导下，合成cDNA第一链，反应体系如下在O.2mlPCR管中依次加入总RNA2X)卩1下游引物0O卩mol/L)1.0卩1dNTPs(10mmol/L)2.0卩1Supe「HI(反婦酶)1.0卩1RNaseInhibito「〔RNA酶抑制剂)1.O卩l5xRTBuffe「(5xRT缓冲液)4.0卩1灭菌重蒸水_9.0MlTotalVolume〔总体枳)20.0pl所述反应体系混合均匀后稍稍离心，将管壁上的液滴收集到管底，将所述反应体系按如下程序运行65°C，5min;4°C，lmin;55°C，60min;7CTC，15min;瞬时离心，放入PCR仪中进行扩增，进行样品cDNA的合成；第四RT-PCR半定量分析合成大豆逆境胁迫诱导基因Gma.2044、Gma.33428、Gma.32516、Gma.34982、Gma.8518、Gma.38、Gma.3774、Gma.22054和内参基因Actin的引物，引物如表l:表l大豆逆境胁迫诱导基因RT-PCR引物<table>tableseeoriginaldocumentpage23</column></row><table>PCR反应总体积为25yl，其中O.2ymolL—^勺前向引物和反向引物，800ymolL—MnTP，1.5mmolL—1MgCl2，1UrTaq聚合酶和2.5y1的10XPCRbuffer,其余体积的物质用去离子水补齐，轻弹离心管底部使所加各试剂充分混匀，稍稍离心将管壁液滴收集到离心管底部，并将所述液滴放置于PCR自动扩增仪进行扩增；取8.0y1PC財广增产物，用2.0%的琼脂糖凝胶进行电泳，在10Vcm—、恒压条件下电泳30min后，在凝胶成像仪下观察并照相，对基因表达进行RT-PCR半定量分析，所述凝胶成像仪为紫外灯；从图2中可以看出，内参基因Actin在4种大豆材料中表达水平相同，说明基因表达的内参基因Actin(表达量参照基因)已经调平，所述4种大豆材料分别为正常培养下大豆，低温、盐和干旱逆境胁迫下大豆，由于在低温，盐和干旱逆境胁迫下的基因表达量与正常培养下的基因表达量相比上升，故所挑选出的大豆逆境胁迫响应基因为与大豆逆境胁迫相关的诱导基因。本实施方式提供了一种对大豆逆境胁迫响应基因的挖掘方法，本实施方式中所述的R平台来源于http:〃www.r-project.org/，本实施方式中所述的大豆转录组文件Gma.gb—cid—lid内容为本实施方式中提取的大豆UniGene转录组文件Gma.gb—cid—1id中大豆逆境胁迫状态下和大豆正常状态下EST序列信息的计数信息为.1000339753351551Qua.10001339753351552Qna.1000433397533515521&na.100060339753351554Qna.10010339753351554Qna.1001373397533515531Qua.1001533397533515512(^na.63397533515524Qna.10020339753351559Qua.100240339753351552从大豆UniGene转录组文件Gma.gb—cid—lid中提取大豆EST序列信息计数信息可利用下述程序实现open(lib—file,〃〈$ARGV〃)；while(defined($in—line=〈lib—file》)chomp$in—line;$library{$in—line}++;close(libfile);■08678墜^卿亟涵i17241，Z.〖涵JL涵雄^i涵,io5S2:2:2iCN-:2;2;oio:o:oio一o;o！o;o:OIojo;lili1:li1;1.i1:2i24open(cid—file,〃〈$ARGV[1]〃)；while(defined($in—line=〈cid—file》)chomp$in—line;if($in—line!八-/)@data=split(/\s+/，line);$gb=$data;$ug="Gma.$data[l]";$library=$data[2];$all—contig{$ug}++;if(!exists($ist—contig{$ug})){$ist—contig{$ug}=0;};{$ist—contig{$ug}++;};close(cid—file);foreach(sort(keys%all—contig))$all—num+=$all—contig{$—};$ist—num+=$ist—contig{$—};foreach(sort(keys%all—contig))print"$—\t$ist—contig{$—}\t$all—num\t$ist—num\t$all—contig{$—}\n差异表达的EST表达量进行超几何分布检验计算，获得的大豆表达基因UniGene转录组在逆境胁迫状态下与在正常状态下差异表达的超几何分布检验值Pialue，所述超几何分布检验值Pialue具体为Gna.20441223515530459^1314.01E—110Gna.33428106351563045981207.19E~90Gna.3251659351563045961056.50E—32Gna.349622035156304596224.04E—20Gna.85182935155304596497.02E—18Gna.381935156304596232.56E—17Gna.37742435155304596383.2犯—1642351553045961312.81E—12Gna.175171435155304596214.90E—11pna一15941235155304596166.43E—11Gna.78SO2235155304596498.舰—11G^ia.508281351553045963801.15E—10Gna.195311035155304596139.29E—10Gna.215841135155304598161.83E~09G叫.268^71635155045984.03E""Q9Gna.142722435155304598674.88E""09Gna.3006120304^98506.81E~09Gna.166441335155304598241.15E""08Gna.220588351553Q4598101.Gna.3507173515530459881.31E""08上述数据的获得方法可以采用实施方式二中的获得人类超几何分布检验值Pialue的方法。具体实施方式四本实施方式所述的一种基于EST数据库和UniGene数据库的基因挖掘方法，与具体实施方式三的不同之处在于步骤G中所述内参基因Actin的调平通过调节RT-PCR模板cDNA的上样量完成，具体为如果PCR仪中扩增产物条带过亮，则按比例减少模板cDNA，如果PCR仪中扩增产物条带较暗，则增加模板cDNA，直到内参基因Actin在4种大豆材料下的扩增条带亮度相同为止，此时内参基因Actin已经调平。权利要求1.一种基于EST数据库和UniGene数据库的基因挖掘方法，其特征在于它的挖掘过程为步骤A下载EST数据库和UniGene数据库，并对所述的UniGene数据库中的EST序列信息按照物种类型进行分类，接下来执行步骤B；步骤B对所述的EST文库注释信息进行信息检索，进而对EST文库注释信息进行分类，接下来执行步骤C；步骤C根据EST文库注释信息的分类信息计算表达基因UniGene转录组的EST表达量，接下来执行步骤D；步骤D对所获得的表达基因UniGene转录组的EST表达量进行超几何分布检验，计算获得表达基因UniGene转录组差异表达的超几何分布检验值P-value，接下来执行步骤E；步骤E采用FDR方法调整表达基因UniGene转录组差异表达的超几何分布检验值P-value，接下来执行步骤F；步骤F设置表达基因UniGene转录组差异表达的超几何分布检验值P-value阈值为0.01，筛选异常状态响应基因，接下来执行步骤G；步骤G利用RT-PCR技术验证所筛选出的异常状态响应基因为与性状相关的诱导基因。2.根据权利要求l所述的一种基于EST数据库和UniGene数据库的基因挖掘方法，其特征在于所述步骤B的具体过程为首先，收集生物体异常状态的关键词，用所收集的每一个生物体异常状态的关键词在EST文库注释信息中进行信息检索，筛选异常状态EST文库，并提取异常状态EST文库的ID，然后，收集生物体正常状态的关键词，用所收集的每一个生物体正常状态的关键词在所述的EST文库注释信息中进行信息检索，筛选正常状态EST文库，并提取正常状态EST文库的ID，所述的在EST文库注释信息中进行信息检索的检索项包括主题TITLE、发育状态DEVELOPMENTAL—STAGE和组织VERBAT頂—TISSUE三项。3.根据权利要求l所述的一种基于EST数据库和UniGene数据库的基因挖掘方法，其特征在于所述步骤C的具体过程为根据EST文库注释信息的分类信息，从所述UniGene数据库中的UniGene转录组文件中提取表达基因UniGene转录组的正常状态EST序列信息的条数和异常状态EST序列信息的条数进行计数，同时对所述UniGene转录组文件中的所有的UniGene转录组的正常状态EST序列信息的条数和所述UniGene转录组文件中的所有的UniGene转录组的异常状态EST序列信息的条数进行计数，将上述所获得的EST序列信息的所有计数信息转化为EST表达量。4.根据权利要求l所述的一种基于EST数据库和UniGene数据库的基因挖掘方法，其特征在于所述步骤C的具体过程为根据UniGene数据库中按照EST文库注释信息的分类信息对UniGene转录组文件中所有的UniGene转录组所完善建立的EST序列信息，直接提取EST序列信息的所有计数信息，将所述计数信息转化为相应的EST表达量。5.根据权利要求l所述的一种基于EST数据库和UniGene数据库的基因挖掘方法，其特征在于所述步骤D的具体过程为设表达基因UniGene转录组的正常状态EST表达量为a，表达基因UniGene转录组的异常状态EST表达量为b，且UniGene转录组文件中所有UniGene转录组的总的正常状态EST表达量为c，所述UniGene转录组文件中所有UniGene转录组的总的异常状态EST表达量为d，以a、b、c和d四个数构建超几何分布检验，利用免费开源软件R平台的phyper完成超几何分布检验，公式一<formula>formulaseeoriginaldocumentpage3</formula>利用公式一计算获得表达基因UniGene转录组差异表达的超几何分布检验值Pialue。6.根据权利要求l所述的一种基于EST数据库和UniGene数据库的基因挖掘方法，其特征在于所述步骤E的具体调整过程为将UniGene转录组文件中所有UniGene转录组差异表达的超几何分布检验值Pialue按照从小到大的顺序依次排列，将p(i)定义为排列后的第i个超几何分布检验值Pialue，同时设所述UniGene转录组文件中所有UniGene转录组的总数为n，所述的i和n均为自然数，且i《n，则调整后的第i个超几何分布检验值P-value为P-value=p(i)*n/i，此等式中所述的P-value和p(i)是重新赋值关系。7.根据权利要求l所述的一种基于EST数据库和UniGene数据库的基因挖掘方法，其特征在于所述步骤F的具体筛选过程为将调整后的超几何分布检验值Pialue阈值以下的超几何分布检验值Pialue所对应的表达基因UniGene转录组作为异常状态的响应基因。8.根据权利要求l所述的一种基于EST数据库和UniGene数据库的基因挖掘方法，其特征在于所述挖掘方法是对人类癌症响应基因的挖掘在步骤A中，下载人类EST数据库和人类UniGene数据库；在步骤B中，直接提取人类EST文库注释信息，首先，提取人类异常状态EST文库，并提取所述的人类异常状态EST文库的ID，然后，提取人类正常状态EST文库，并提取所述的人类正常状态EST文库的ID;在步骤C中，根据人类EST文库注释信息的分类信息，从所述人类UniGene数据库中的人类UniGene转录组文件Hs.profiles中提取人类表达基因UniGene转录组的健康状态HealthState下乳房月中瘤breasttumor，白血病leukemia，卵巢月中瘤ovariantumor，中枢神经系统CNS的原始神经夕卜胚层月中瘤primitiveneuroectodermaltumor,前歹[J腺癌prostatecancer和正常状态normal的人类EST序列信息的条数进行计数，同时对人类UniGene转录组文件Hs.profi1es中所有的UniGene转录组的正常状态人类EST序列信息的条数和人类UniGene转录组文件Hs.profiles中所有的UniGene转录组的异常状态人类EST序列信息的条数进行计数，将上述所获得的人类EST序列信息的所有计数信息转化为人类表达基因UniGene转录组的EST表达量；在步骤D中，对所获得的人类表达基因UniGene转录组的EST表达量进行超几何分布检验，其中，以人类UniGene转录组文件Hs.profiles中正常状态normal下的数据作为人类表达基因UniGene转录组的正常状态EST表达量，以人类UniGene转录组文件Hs.profiles中其它癌症病变状态作为异常状态，并将所述异常状态下的数据作为人类表达基因UniGene转录组的异常状态EST表达量，利用所获得的人类表达基因UniGene转录组的正常状态EST表达量、人类表达基因UniGene转录组的异常状态EST表达量、人类UniGene转录组文件Hs.profiles中所有UniGene转录组的总的正常状态EST表达量和人类UniGene转录组文件Hs.profiles中所有UniGene转录组的总的异常状态EST表达量构建超几何分布检验，利用免费开源软件R平台的phyper完成超几何分布检验，计算获得人类表达基因UniGene转录组差异表达的超几何分布检验值Pialue;在步骤E中，采用FDR方法调整人类表达基因UniGene转录组的超几何分布检验值Pialue，该具体调整过程为将人类UniGene转录组文件Hs.profiles中所有UniGene转录组差异表达的超几何分布检验值Pialue按照从小到大的顺序依次排列，将p(i)定义为排列后的第i个超几何分布检验值Pialue，同时设人类UniGene转录组文件Hs.profiles中所有UniGene转录组的总数为n，所述的i和n均为自然数，且i《n，则调整后的第i个超几何分布检验值P-value为P-value=p(i)*n/i，此等式中所述的P-value和p(i)是重新赋值关系；在步骤F中，设置人类癌症所响应的表达基因UniGene转录组差异表达的超几何分布检验值Pialue阈值为O.Ol，筛选人类癌症响应基因，在此筛选过程中，将调整后的超几何分布检验值Pialue阈值以下的超几何分布检验值Pialue所对应的表达基因UniGene转录组作为人类癌症的响应基因；在步骤G中，利用RT-PCR技术验证所筛选出的人类癌症响应基因为与人类癌症相关的诱导基因。9.根据权利要求l所述的一种基于EST数据库和UniGene数据库的基因挖掘方法，其特征在于它是对大豆逆境胁迫响应基因的挖掘在步骤A中，下载大豆EST数据库和大豆UniGene数据库；在步骤B中，以"cold"，"salt"和"drought"作为与大豆逆境胁迫相关的异常状态的关键词，将所述关键词在大豆EST文库注释信息文件Gma.lib.info中进行信息检索，筛选大豆异常状态EST文库，并提取大豆异常状态EST文库的ID，然后，将大豆EST文库中的其他EST文库作为大豆正常状态EST文库，并提取大豆正常状态EST文库的ID，所述的在大豆EST文库注释信息文件Gma.lib.info中进行信息的检索，检索项包括所述EST文库注释信息文件Gma.lib.info中的主题TITLE、发育状态DEVELOPMENTAL—STAGE和组织VERBAT頂—TISSUE三项；在步骤C中，根据大豆EST文库注释信息的分类信息，从所述大豆UniGene数据库中的大豆UniGene转录组文件Gma.gb—cid—lid中提取大豆表达基因UniGene转录组的正常状态EST序列信息的条数和异常状态EST序列信息的条数进行计数，同时对大豆UniGene转录组文件Gma.gb—cid—lid中所有的UniGene转录组的正常状态大豆EST序列信息的条数和大豆UniGene转录组文件Gma.gb—cid—lid中所有的UniGene转录组的异常状态大豆EST序列信息的条数进行计数，将上述所获得的大豆EST序列信息的所有计数信息转化为大豆表达基因UniGene转录组的EST表达量；在步骤D中，对所获得的大豆表达基因UniGene转录组的EST表达量进行超几何分布检验，其中，以大豆UniGene转录组文件Gma.gb—cid—lid中正常状态normal下的数据作为大豆表达基因UniGene转录组的正常状态EST表达量，以大豆UniGene转录组文件Gma.gb—cid—lid中异常状态下的数据作为大豆表达基因UniGene转录组的异常状态EST表达量，利用所获得的大豆表达基因UniGene转录组的正常状态EST表达量、大豆表达基因UniGene转录组的异常状态EST表达量、大豆UniGene转录组文件Gma.gb—cid—1id中所有UniGene转录组的总的正常状态EST表达量和大豆UniGene转录组文件Gma.gb—cid—lid中所有UniGene转录组的总的异常状态EST表达量构建超几何分布检验，利用免费开源软件R平台的phyper完成超几何分布检验，计算获得大豆表达基因UniGene转录组差异表达的超几何分布检验值Pialue;在步骤E中，采用FDR方法调整大豆表达基因UniGene转录组差异表达的超几何分布检验值Pialue，该调整过程为将大豆UniGene转录组文件Gma.gb—cid—lid中所有UniGene转录组差异表达的超几何分布检验值Pialue按照从小到大的顺序依次排列，将p(i)定义为排列后的第i个超几何分布检验值Pialue，同时设大豆所有UniGene转录组的总数为n，所述的i和n均为自然数，且i《n，则调整后的第i个超几何分布检验值Pialue为Pialue=p(i)*n/i，此等式中所述的Pialue和p(i)是重新赋值关系；在步骤F中，设置大豆逆境胁迫所响应的表达基因UniGene转录组差异表达的超几何分布检验值Pialue阈值为O.Ol，筛选大豆逆境胁迫响应基因，在此筛选过程中，将调整后的超几何分布检验值Pialue阈值以下的超几何分布检验值Pialue所对应的表达基因UniGene转录组作为大豆逆境胁迫的响应基因；在步骤G中，利用RT-PCR技术验证所筛选出的大豆逆境胁迫响应基因为与大豆逆境胁迫相关的诱导基因，具体验证过程为第一挑选步骤F中所述的大豆逆境胁迫响应基因进行RT-PCR验证具体挑选过程为将步骤F中筛选出的大豆逆境胁迫的响应基因按照超几何分布检验值Pialue从小到大的顺序排列，挑选出最小的8个超几何分布检验值Pialue所对应的大豆逆境胁迫的响应基因，所述8个超几何分布检验值Pialue所对应的大豆逆境胁迫的响应基因分别为Gma.2044、Gma.33428、Gma.32516、Gma.34982、Gma.8518、Gma.38、Gma.3774和Gma.22054;第二大豆总RNA提取以正常培养的大豆21日龄幼苗进行4。C低温cold、200mmolL-1NaCl盐salt和100mmolL-l甘露醇干旱drought胁迫处理l小时后的大豆21日龄幼苗叶片为材料，采用TRIZOL试剂法提取大豆总RNA，所述采用TRIZOL试剂法提取总RNA的方法参照试剂盒说明书进行，具体为i.称取100mg所述幼苗叶片样品放入高温高压消毒研钵，迅速在液氮中研磨成粉末，然后加入lmlTRIZOL研磨成匀浆组织；ii.将所述匀浆组织作为标本转移到1.5ml的离心管中，在153(TC环境下放置5min，以彻底分离核蛋白复合体；iii.向所述离心管中加入0.2ml的氯仿，加盖好后用手剧烈摇晃15s，在1530。C环境下放置23min，然后离心12000rmin-1，15min，28。C，离心后离心管内物质分成三层，下层为红色的酚，中间层为氯仿相，上层为无色的水相，RNA只存在于水相中，水相占总TRIZOL的60o/o;iv.将上层水相转移到另一个干净的离心管中，加入0.5ml异丙醇，在153(TC环境下静置IOmin，然后离心12000rmin-1，10min，28。C;v.去步骤iv中所述离心管中上清，加入lml75。/。乙醇洗涤RNA沉淀，采用振荡器混匀，然后离心7500r*min-1，5min，28°C;vi.继续去步骤v中离心管中上清，并将去上清后的离心管置于真空或空气中510min，干燥RNA沉淀，加DEPC处理的无菌水15y1，并将加入无菌水的离心管置于-2(TC环境下保存备用，提取总RNA完成；第三cDNA的合成以所提取的总RNA为模板，在下游引物的引导下，合成cDNA第一链，反应体系如下在O.2mlPCR管中依次加入总RNA2X)卩1下游引物0O卩mol/L)1.0卩1dNTPs(10mmol/L)2.0卩1Supe「HI(反婦酶)1.0卩1RNaseInhibito「〔RNA酶抑制剂)1.O卩l5xRTBuffe「(5xRT缓冲液)4.0卩1灭菌重蒸水_9.0MlTotalVolume〔总体枳)20.0pl所述反应体系混合均匀后稍稍离心，将管壁上的液滴收集到管底，将所述反应体系按如下程序运行65°C，5min;4°C，lmin;55°C，60min;7CTC，15min;瞬时离心，放入PCR仪中进行扩增，进行样品cDNA的合成；第四RT-PCR半定量分析合成大豆逆境胁迫诱导基因Gma.2044、Gma.33428、Gma.32516、Gma.34982、Gma.8518、Gma.38、Gma.3774、Gma.22054和内参基因Actin的引物，引物如表l:表lt豆.細办迫i紹細RT-PCR引物_站E"I时向引物反向3f物Gma-2044MACTTTGACTGGCAAfiACCATTATCTGAACTCTO:ACCTCCAAGG恥GAAGTGAAC下CAGACMGACCCfflXAAGCTGAGAGAGGAAACcOTCTGcrGAAAGAGA肌虹ACCACTACCTTCACAACGma.33ACATGCCTTCTACAACACCOTCKTGC(TCAtCCCT(TATAnTG亂3774CTGGTTCTATGCCACCrTCTT匚TCTCCTCTGTATTTTCTCCTCGGTG亂22054ActknCKAGCACTG恥TCATCACAACTACTGCTKAGCAGTGKAAATGTPCR反应总体积为25y1，其中0.2ymo1L-l的前向引物和反向引物，800ymolL-ldNTP，1.5mmolL-lMgC12，1UrTaq聚合酶和2.5y1的10XPCRbuffer，其余体积的物质用去离子水补齐，轻弹离心管底部使所加各试剂充分混匀，稍稍离心将管壁液滴收集到离心管底部，并将所述液滴放置于PCR自动扩增仪进行扩增；取8.0y1PC財广增产物，用2.0%的琼脂糖凝胶进行电泳，在10Vcm-l恒压条件下电泳30min后，在凝胶成像仪下观察并照相，对基因表达进行RT-PCR半定量分析，所述凝胶成像仪为紫外灯；内参基因Actin在4种大豆材料中表达水平相同，说明基因表达的内参基因Actin已经调平，所述4种大豆材料分别为正常培养下大豆，低温、盐和干旱逆境胁迫下大豆，由于在低温，盐和干旱逆境胁迫下的基因表达量与正常培养下的基因表达量相比上升，故所挑选出的大豆逆境胁迫响应基因为与大豆逆境胁迫相关的诱导基因。10根据权利要求9所述的一种基于EST数据库和UniGene数据库的基因挖掘方法，其特征在于步骤G中所述内参基因Actin的调平通过调节RT-PCR模板cDNA的上样量完成，具体为如果PCR仪中扩增产物条带过亮，则按比例减少模板cDNA，如果PCR仪中扩增产物条带较暗，则增加模板cDNA，直到内参基因Actin在4种大豆材料下的扩增条带亮度相同为止，此时内参基因Actin已经调平。全文摘要一种基于EST数据库和UniGene数据库的基因挖掘方法，它涉及应用生物信息学领域。它克服了现有的基因挖掘方法中无法正确挖掘与性状相关的诱导基因的问题。本发明的方法利用EST数据库中的EST序列将UniGene数据库中的表达基因UniGene转录组的EST表达量数字化，构建超几何分布检验，并结合FDR方法调整表达基因UniGene转录组的差异表达的超几何分布检验值P-value，筛选异常状态响应基因，最后利用RT-PCR技术验证所述响应基因为与性状相关的诱导基因。本方法可以用于人类疾病发生、动植物生长发育调控、动植物疾病调控过程以及动植物逆境胁迫等性状相关诱导基因的挖掘。文档编号G06F19/00GK101661536SQ20091030848公开日2010年3月3日申请日期2009年10月20日优先权日2009年10月20日发明者季佐军,华才,朱延明,勇李,束永俊,锡柏,巍纪申请人:东北农业大学