畜牧兽医学报  2019, Vol. 50 Issue (3): 474-484. DOI: 10.11843/j.issn.0366-6964.2019.03.002    PDF    
转录组数据分析与功能基因挖掘
李欣1, 李小俊1, 陈晓丽1, 赵毅强2, 王栋1     
1. 中国农业科学院北京畜牧兽医研究所, 北京 100193;
2. 中国农业大学生物学院, 北京 100193
摘要:高通量测序技术的不断发展和应用,为挖掘重要功能基因提供了转录组分析方法,但如何利用海量测序数据准确、高效地挖掘功能基因,仍是转录组学分析方法研究的重要瓶颈。本文综述了RNA-seq数据质量控制与读段定位、基因组注释、转录本拼接、表达水平评估、差异表达分析等环节分析方法,比较了数据分析常用软件、算法和数据库等的性能和适用范围;同时,又综述了蛋白调控互作网络和加权基因共表达网络等差异表达基因的功能分析方法。转录组分析正在从只利用物种内信息挖掘差异基因,向引入其他物种参考系进行目标物种功能基因挖掘分析方向发展。结合同源基因预测候选基因法、选择信号法、极端数据法、GO注释和KEGG富集分析法及BSR-Seq(bulked segregant RNA-Seq)法等鉴定方法,使分析结果更加科学可靠。随着测序技术和数据分析方法不断进步、数据库资源不断完善,测序数据中隐含的基因表达调控和生命规律将会逐渐得到准确、深入揭示。
关键词转录组    数据分析    基因挖掘    
Transcriptomics Analysis and Functional Genes Mining
LI Xin1, LI Xiaojun1, CHEN Xiaoli1, ZHAO Yiqiang2, WANG Dong1     
1. Institute of Animal Science, Chinese Academy of Agricultural Sciences, Beijing 100193, China;
2. College of Biological Sciences, China Agricultural University, Beijing 100193, China
Abstract: With the continuous development and application of high-throughput sequencing technology, transcriptome analysis method is developed for mining genes with important function. However, a lot work needs to be done for efficient and accurate transcriptome analysis based on massive sequencing data. Here, we reviewed methods for reads quality control, reads mapping, genome annotation, transcripts assembling, expression quantification, differential expression analysis for RNA-seq data. We summarized the performance and scope of application of the common softwares, algorithms and databases used. We also reviewed analysis methods such as protein regulatory interaction networks as well as weighted gene co-expression networks. Transcriptome analysis has been evolved from identifying differentially expressed gene within-species to utilizing related species as reference to mine the functional genes in target species. By combining with various methods, such as the homologous gene prediction, select signal detection, extreme data analysis, GO annotation and KEGG enrichment and bulked segregant RNA-Seq (BSR-Seq) methods, the results from RNA-seq analysis are more scientific and reliable. With the development of sequencing technology and data analysis methods as well as continuous improvement of database resources, the underline gene regulation and the law of life implied in the sequencing data will be uncovered accurately and deeply in future.
Key words: transcriptome     data analysis     genes mining    

1997年,Velculescu等[1]在对酵母的研究中首次提出转录组(transcriptome)概念,转录组是指特定物种特定组织或细胞转录的所有RNA集合,可揭示基因表达的时空动态性,反映生物个体在某一特定生长发育阶段的特定细胞、组织或器官所有基因的转录表达水平[2-4],同时,也可用来比较某一器官、组织或细胞在不同环境条件下的基因表达差异[5-6]。转录组研究开启了畜禽重要性状遗传解析与功能基因挖掘及鉴定研究的新纪元。

然而,测序数据的不断涌现并没有使畜禽功能基因和调控通路等挖掘效率有较大改进和提高,对测序数据进行表达量验证后,再对众多差异表达基因进行GO(gene ontology, GO)[7]和KEGG(Kyoto encyclopedia of genes and genomes, KEGG)[8]富集分析。许多科研人员都在各自的研究中进行着类似的重复。虽然这些研究的确为后续功能基因挖掘提供了重要参考依据,但是,转录组数据耗费了大量人力、物力和财力,这些分析远未使其得到深度挖掘。如何从海量序列信息中深入挖掘畜禽重要经济性状的调控通路、理解生命现象的本质及疾病发生发展的过程,成为生物信息学研究工作中亟待解决的问题。为深入挖掘测序数据信息,除进行物种内转录组数据分析外,研究人员又引入其他物种测序数据,希望借助物种间数据信息深度挖掘转录组数据。比如Li等[9]通过比较果蝇和秀丽线虫转录组相似性,探索了不同发育阶段转录组数据研究方法,发现从早期胚胎到晚期幼虫两个物种基因表达呈强共线性,以某一物种转录组数据为参考,预测另一物种同一发育过程有共线性特征的目标表达基因,这一研究方法开启了利用比较转录组研究方法深入挖掘转录组数据信息的新思路。陈巍[10]通过比较极端数据集同源基因表达的方法,寻找人脑中与小鼠大脑转录基因表达与功能相似及相异的基因,为构建人神经退行性疾病小鼠模型,并评价其有效性和准确性提供理论支持。本文综述了利用转录组数据挖掘鉴定畜禽重要功能基因的研究方法,旨在为构建科学高效的转录组测序研究方法,实现更深层次的转录组数据挖掘提供理论指导。

1 RNA-seq数据分析流程及方法

要进行RNA-seq数据分析,首先要制备目标组织总RNA,对其RNA样品片段化处理,并反转录成为cDNA片段,构建测序文库,序列测定后,采用类Unix(Linux或者MacOSX)64位操作系统进行数据处理。其分析流程如图 1所示。

图 1 RNA-seq数据分析流程图 Fig. 1 The flow chart of RNA-seq data analysis
1.1 原始数据质控

序列测定后,以FASTQ格式输出原始数据(raw data),使用Trimmomatic软件[11]过滤掉raw data的接头(adapter)、低质量、重复及未测出序列,得到clean data。再用FastQC软件[12]检测clean data的碱基质量值(quality score, Q-score)和碱基分布情况,确定测序及过滤效果,并用碱基测序错误率(Q20与Q30)表示质量的高低。Q20和Q30分别表示测序错误率为1%和0.1%的碱基比例。GC含量及样品间相关性也用于原始测序质量鉴定。序列质量控制软件主要有FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/)、Trimmomatic等。通过质量控制才能使基因挖掘和功能分析结果科学可靠。

1.2 Clean Reads的比对与拼接

为序列组装而进行的序列比对方法有以下两种:1)有参考基因组序列的测序数据组装时,先将所有reads通过序列映射定位(mapping)到参考基因组上;再把比对到相应位置的reads聚类,形成一张表示所有可能可变剪切形式的图谱(graph);最后转换graph信息为转录本信息。该方法先从UCSC数据库(http://genome.ucsc.edu/)下载参考基因组序列,并利用bowtie[13]软件将其构建为基因组索引(index)文件。然后,利用TopHat和RNASEQR软件在构建好的index文件、全基因组序列信息及参考基因注释信息指导下,将质控测序读段定位到基因组上。常用的序列比对软件主要有Bowtie2[14]、Tophat2[15]、STAR[16]、GSNAP[17]、SOAP2[18]、S-MART[19]等。2)无参考基因组序列的序列组装也称为从头测序组装(de novo sequence assembly),是将各测序读长按顺序拼接成连叠群(contig),再组装成支架(scaffle),最后填充支架中间空缺的部分gap,组装成连续的长序列,并通过与模式动植物序列比对(basic local alignment search tool, BLAST),确定基因序列。以Trintiy[20]软件为代表的de novo组装平台,为转录组从头组装提供了有效工具[21]。该软件包含Inchworm、Chrysalis和Butterfly 3个独立模块,3个软件模块依次运行,将reads组装成完整的转录本,其拼接策略是将clean data分割再拼接,得到多个独立的de Bruijn图(基因转录产物),随后通过reads溯源将这些de Bruijn图分类,最终得到全长转录本,并根据图分类分辨旁系同源基因。从头组装常用的软件有SOAPdenovo (http://soap.genomics.org.cn/soapdenovo.html)、Velvet (https://www.ebi.ac.uk/zerbino/velvet/)、Edena (http://www.genomic.ch/edena.php)、Abyss (http://www.bcgsc.ca/platform/bioinfo/software/abyss)、ALLPATH (http://genome.cshlp.org/content/18/5/810.full.html)等。

由于真核细胞中每个基因可产生多个转录本,导致同一个测序片段可能会同时出现在多个转录本中,并使RNA拼接出多个互不连通的图结构,每个图对应该基因相应的转录产物。进行序列拼接时,两种方法都希望将算法问题简化为对每个基因的拼接问题,并提出了剪接图、重叠图以及de bruijn graph等很多图模型,以实现每个图与相应基因的一一对应。两种方法采用了不同的拼接策略,各有各的使用范围,互为补充。一般来说,由于没有参考信息,又受到测序错误、覆盖度不均匀等因素影响,从头拼接算法准确率明显低于基于参考基因组的拼接算法。但有些物种没有完整的基因组序列,从头拼接方法对这些物种具开创意义。另外,尽管某些物种具有完整基因组序列,由于某些疾病等因素,使得其基因组发生严重的变异及缺失,从头拼接法将显示出明显优势[22]

虽然经历了几十年研究,序列拼接仍未得到令人满意的解决方案[23-24]。由于repeat、序列拼接的时间和空间挑战性、序列片段中的错误(杂质序列和碱基读出错误)、contig间的位置及距离不能准确确定等问题[25],常导致GenBank数据库中难以避免的拼接错误[26]。柳军涛[27]针对转录组拼接精度和效率问题,设计出基于参考基因组的拼接算法TransComb和从头拼接算法BinPacker,灵敏度和准确度均有所提高,并提高了计算速度,减少了内存消耗。王春宇[28]针对高通量测序导致的数据读长短、数量大和错误率高等问题,基于MapReduce和序列聚类提出了SeedsGraph方法,对全基因组de novo拼接方法进行了改进,对较大规模基因组拼接问题进行了深入研究。

1.3 转录本表达分析

已构建的转录组分析方法,可根据比对到基因组上的序列数量,对各基因表达进行均一化处理、定量估计表达情况和差异表达分析。常用软件有Cufflinks/Cuffdiff2[29]、SAMtools[30]、RSEM[31]、DESeq[32]、edgeR[33]等,利用这些软件统计表达基因定位的reads数量、基因长度及结构,得到该基因的RPKM(reads per kilobase per million mapped reads)或FPKM(fragments per kilobase of exon model per million mapped reads)值,进行表达差异分析。

通过TopHat软件分析,得到定位在基因组上来自同一转录本的片段,再利用Cufflinks(http://cole-trapnell-lab.github.io/cufflinks/)软件将这些片段组装为全长转录本[34]。进而,使用Cuffcompare软件将全长转录本与已有基因组注释文件进行比较,评估转录本构建情况,并根据已知Ensembl数据库中转录本信息,定义所构建转录本的内含子、外显子、基因间等区域。每组样品都可拼接得到各自的转录本信息,可利用Cuffmerge软件将它们合并为一个转录本集合,作为下一步差异表达信息分析的依据。Cufflinks中包含的Cuffdiff软件可用来计算两个或更多样品的基因表达,并且针对基因表达量,对表达丰度进行统计分析及统计学检验,获得不同样本间差异表达RNA分子。也可利用DEGSeq R package v1.10.1[35]进行不同处理间差异表达分析,通过计算同一基因在两个处理中表达量相等的P值,再采用多重假设检验对P值进行校正,校正后的P[36]用于差异表达基因显著性检验。

进行可变剪接分析时,首先利用TopHat、SOAPsplice(http://soap.genomics.org.cn/soapsplice.html)、SpliceR (http://www.bioconductor.org/packages/release/bioc/html/spliceR.html)、SpliceGrapher (http://splicegrapher.sourceforge.net/)、Asprofile(https://omictools.com/asprofile-tool)和Astalavista v2.2[37]等软件将reads比对到基因组序列后,根据位置、长度及结构信息预测剪接体的类型,然后用Cufflinks软件的isoform center程序包进行剪切体差异表达分析。

目前,主要有两种转录组序列剪接位点识别算法,分别为以基因组序列潜在位点为中心的剪接识别,如TopHat算法,以及以reads分割匹配为中心的剪接位点识别,如MapSplice算法。因为TopHat能够同时利用多机多核CPU资源,并行化运行,使得分析效率大幅提升,但其剪接模式搜索依赖于对基因组潜在位点的预识别,仅能实现对特定距离内保守剪接方式的探测。MapSplice(http://www.netlab.uky.edu/p/bioinfo/MapSplice)完全以reads为中心,进行剪接位点搜索,不受保守剪接方式限制,但由于基因组本身重复序列较多,并且测序过程中形成了过多短片段,影响了计算效率和匹配精度,降低了剪接位点的可靠性。常用的可变剪切数据库有H-InvDB(http://www.h-invitational.jp/)、MAASE(http://maase.genomics.purdue.edu/)、SpliceInfo(http://spliceinfo.mbc.nctu.edu.tw/)、TassDB(https://www.tassdb.info/)、SpliceAid(http://www.introni.it/splicing.html)、MiasDB(http://47.88.84.236/Miasdb/index.php)和EuSplice(https://www.genome.com/)等。

1.4 基因功能富集分析

基因功能指的是众多代表一定功能特征和生物过程的基因功能集[38-39]。通过基因功能富集分析,可将成百上千个不同功能的基因、蛋白或其它分子聚类到不同生物学通路,其中,生物体内实现某种功能性状的一组基因或蛋白等,常被富集到某一条通路中。基因功能富集分析减少了后续分析的复杂度,还可发现生物学过程中起关键作用的生物学通路,有助于揭示和理解生物学过程基本分子机制。包含这些基因功能集的常用基因功能数据库有GO、KEGG、Reactome(https://reactome.org/)[40-41]、BioCarta(https://cgap.nci.nih.gov/Pathways/BioCarta_Pathways)[42]、MsigDB[43-44]、wikipathways(https://www.wikipathways.org/index.php/WikiPathways)[45]、geneDB(http://www.genedb.org/Homepage)、ERGO(https://ergo.integratedgenomics.com/login.ergo)、GenMAPP(http://www.genmapp.org)[46]等。

KEGG数据库向公众开放了最著名的生物学通路资源,可系统分析基因产物的细胞代谢途径及功能,通过KEGG分析实现表达基因的pathway富集,并获得相应注释。该分析可通过KOBAS[47-48]或DAVID等软件实现,KEGG网站为每一种生物学通路提供了专门的图示说明。BioCarta数据库提供了生物学通路绘制模板,研究者可根据BioCarta数据库要求和规定的标准上传相应的生物学通路分析结果,但该数据库不负责对用户上传的生物学通路分析结果进行质量检验,导致该数据库资源质量参差不齐,且有许多重复。尽管如此,该数据库数据资源量巨大,且不同于KEGG数据库,其包含了大量代谢通路之外的生物学通路,对于有较强甄别能力的科研人员来说,有重要的参考价值。另外,基因功能富集分析还采用GenMAPP数据库,该数据库可免费使用,使用者可利用数据库资源自己绘制生物学通路,并以MAPP格式保存,这个文件很小,易于在网络上传播,有利于研究者间随时交流。

1.5 蛋白质互作网络(protein-protein interaction network, PPI network)分析

分析差异基因表达蛋白参与的生物信号传递、基因表达调控、能量和物质代谢及细胞周期调控等生命过程,不仅可了解蛋白间的功能联系,也可了解疾病等特殊生理状态下生物信号和能量物质代谢的反应机制[49-50]。常用数据库主要有SWISS-PROT(https://web.expasy.org/docs/swiss-prot_guideline.html)、TrEMBL(https://www.ebi.ac.uk/uniprot)、PDB(https://www.rcsb.org/)、SDSPB(http://lifecenter.sgst.cn/)和String(https://string-db.org/)等。将差异表达基因数据导入Cytoscape[51-52]软件,可实现对蛋白质互作网络的可视化分析。对于数据库中没有注释信息的物种,需将候选基因序列与数据库收录的亲缘较近物种的基因进行blastx比对,进而借助该物种在数据库中已经存在的蛋白质相互作用关系,预测候选靶基因间的PPI。这种方法首先需要根据文献信息确定候选分子间是否或可能存在相互关系,再利用Cytoscape软件构建互作网络并分析关键节点,最终确定PPI。Li等[53]在鸡上应用PPI分析鉴定出POMC基因是网络hub基因,且找到了该基因与摄食pathway相关的基因集。然而,目前报道的生物医学研究大多基于还原论,研究方法仅关心个别蛋白质或基因,不足以从整体水平系统地解释生物行为。网络系统生物学可利用网络分析方法,从整体水平研究生命体内各种相关生命活动间的相互关系及其动态变化,为全面揭示生命奥秘提供了重要的研究手段。网络系统生物学整合分析文献挖掘得到的信息弥补了还原论分析方法的不足,更有利于深入理解生命现象的科学本质。

1.6 加权基因共表达网络分析(weighted gene co-expression network analysis, WGCNA)

WGCNA可找出协同表达的基因、基因网络与表型间的相关关系,以及基因网络中的核心基因。适用于不同器官或组织发育调控、同一组织不同时间发育调控等复杂数据模式的转录组数据挖掘[54]。WGCNA分为表达量聚类分析和表型关联两个环节,含基因间相关系数计算、基因模块确定、共表达网络构建、模块与性状关联4个步骤。通过上述的网络分析找到处于调控中心的hub基因,再深入分析该基因并挖掘其功能。由于传统表达调控网络分析方法以基因间表达量相关系数为基础,样本数过低时,相关系数可靠性差,得到的调控网络价值不大。然而,WGCNA摒弃了通过设定相关系数阈值判断基因间“相关”或“不相关”的做法,对基因间表达量的相关系数进行适当加权,避免了信息损耗,使相应基因表达调控网近似服从无尺度网络分布。WGCNA的有效性已经得到了实践验证,利用该技术成功鉴定了猪卵母细胞与脂肪组织的脂质沉积共表达模式[55],并在比较牛和鸡转录组后找出了鸡SYFs细胞表达的调控因子ESR2[56],目前WGCNA还在不断完善,其功能也将越来越强大。

2 转录组测序数据的基因发掘

通过RNA测序和对测序数据的生物信息学分析,可鉴定重要畜禽器官、组织或细胞中各种生物学过程的差异表达基因及其调控关系,并发现新转录本,鉴别新基因。近年来,随着测序技术不断取得新突破,转录组测序数据如井喷式增长,针对当前转录组测序数据发展的迅猛势头,如何高效、精准地从海量测序数据中筛选目标候选基因和发掘新基因,已成为本领域研究的焦点。除利用质量性状表型差异寻找转录基因的常规方法[57-58]外,也研发出了一些转录组特有的挖掘方法。在挖掘策略上,从仅利用物种内数据进行目标基因挖掘,发展到了引入其他物种测序数据参考系指导主要畜禽重要基因挖掘与功能鉴定。

2.1 利用同源基因预测候选基因功能

对转录组数据分析得到的差异基因进行富集分析,找出富集在相关通路上的候选基因,然后在NCBI上找到该基因或蛋白序列,再利用序列同源性比较预测新基因功能。序列同源性比较的基本假设为:如果基因A与基因B同源,基因A可能就具有类似基因B的功能。利用同源比对算法,在DNA或蛋白质序列数据库中检索待测序列同源基因,得到一系列与该基因同源性较高的基因或片段,这些基因或片段的已知功能信息为进一步研究该基因功能提供了导向。主要有同源检索和多序列比对两种方法:1)同源检索:可通过OrthoMCL (http://orthomcl.org/orthomcl/)软件进行同源基因检索,序列检索时,根据序列的结构和可能功能,通过直系同源基因和旁系同源基因两种途径,采用软件对待检索序列进行搜索,通常认为,同一功能域的基因序列即使来源于不同物种也可定义为直系同源基因,而分属于不同功能域的基因序列,即使来源于同一物种也被定义为旁系同源。同源检索可为基因组进化分析提供证据,并预测候选基因功能。但是,基于该软件的直系同源基因检索必须要在比较完整的基因组之间进行;另外,也可基于隐马可夫模型,使用HMMER3.0软件检测序列比对结果中的保守区域,识别出序列中已知的核酸或蛋白质结构域,阐明序列间的超家族、家族、亚家族和种属特异性等关系,其优点是可检测到远源同源序列。同源检索为发掘进化或物种分化过程中的关键功能基因提供了重要工具。2)多序列比对:是同时对多个序列进行同源比较,以发现其共同结构特征,该方法为寻找基因家族或蛋白质家族保守区域提供了重要技术支撑。保守区域与家族成员的功能密切相关,通过这些方法建立的蛋白质家族数据库,有助于新基因所属蛋白家族及其保守区域的搜寻,并提供这个家族其他成员的结构和功能信息。以对燕雀的分析为例,Wang等[59]通过比较3类燕雀转录本、蛋白序列,提出了研究线粒体复合酶I进化的新方法。多序列比对软件主要有ClustalW/X(http//www.ebi.acuk/clustalw)、PRALINE(http://www.ibi.vu.nl/programs/pralinewww)、Multalin(http://www-archbac.u-psud.fr/genomics/multalin.html)、SAGA(http://www.saga-gis.org/en/index.html)、MEGA7 (http://www.megasoftware.net/)[60]、MSA[61]、MAFFT(http://mafft.cbrc.jp/alignment/software/)[62]。这些多序列比对软件是应用各种程序开发的,其中,使用最广泛的比对程序为基于渐进比对算法的ClustalW程序,但其精度不高;因此,又相继研发了其他比对程序,其中,比较有代表性和影响力的有基于迭代细化策略的PRRN/PRRP、基于傅立叶变换的MAFFrJ和基于多重迭代的MUSCLE和T-COFFEEJ[63]等,这些新程序提升了多序列比对的速度和精度。

2.2 利用选择信号法挖掘基因

物种形成经历了自然选择或自然和人工双重选择,长期处于定向选择压力作用下生物的某些性状表型将发生定向性改变,并使基因组某些区域基因组成发生质和量的改变,进而导致不同等位基因频率的变化,甚至产生新突变,基因组中这种由于选择作用留下的印迹被称为选择信号(selection signature)[64]。利用选择信号法对转录组数据进行功能基因挖掘时,首先利KisSplice、Gatk或samtools等软件,将转录组测序得到的SNP(sigle nucleotide polymorphism)数据筛选出来,再结合XP-EHH(cross-population extend haplotype homozygosity)法对SNP数据进行检测,获得群体选择信号,然后通过location and extension法[65-67]对群体选择信号进行扫描,最终获得相关性状的基因组选择区域。因为每个选择区域有许多基因,首先计算每个基因中各SNP位点的XP-EHH score,然后进行XP-EHH score由低到高排序,将该基因分值最高SNP位点的得分定为其受到的选择强度,最后,将该区域中XP-EHH score最高的SNP位点得分定为该选择区域受到的选择强度,该区域中最高XP-EHH score SNP位点所在的基因即为该选择区域候选基因。构建系统进化树分析筛选出的候选基因,挖掘出与物种类群和目标性状密切相关的受选择基因,同时,将其与自然选择和人工选择联系起来,进一步说明受选择基因在当时环境和进化中的作用。通过构建中国荷斯坦和西门塔尔牛群体选择信号图谱、欧洲猪与中国猪群体分化与混合[68]等实例说明,选择信号法可用于研究家畜驯化和品种形成机制,解释不同群体经历的育种过程和进化历史,并进一步揭示与表型相关的基因变异。

2.3 利用极端数据挖掘基因

健康组织中高表达基因对于生物体充分发挥生理功能具有潜在重要意义[69-70],由此,构建了通过比较极端数据集挖掘畜禽同源奢侈基因的方法。同源基因在不同畜禽中可能存在序列差异,表达量也会存在差异。利用百分位数法将同源基因集按照表达量高低分为(表达量位于前10%的基因集为高表达组,位于10%~50%的基因集为中表达组,位于50%~90%的基因集为低表达组,位于后10%的基因集为极低表达组)4组,定义极低表达组与高表达组为极端数据集。将A畜种高表达组与B畜种对应组织极低表达组共有的同源基因进行比较,两者交集为A畜种表达的同源奢侈基因,即该基因在A畜种中高表达,而在B畜种对应组织极低表达或不表达。同样,可筛选出在B畜种组织特异高表达,而在A畜种对应组织极低表达或不表达的同源基因。目前,百分位数法已被广泛用于挖掘相对高表达或物种特异表达同源基因[69-70],但该方法还因测序平台、比对方法和软件等差异易产生一定系统误差,导致假阳性结果比例较高,影响基因挖掘的准确性,划分高表达基因的标准不合适也可影响该方法的挖掘效果。

2.4 通过GO注释和KEGG富集分析挖掘基因

基于GO和KEGG数据库发展起来的David(https://david.ncifcrf.gov/)、GOEAST(http://omicslab.genetics.ac.cn/GOEAST/)、GOSim(http://www.bioconductor.org/packages/release/bioc/html/GOSim.html)、KEGGArray(http://www.kegg.jp/kegg/download/kegtools.html)、PathwaryMiner(http://apps.cytoscape.org/apps/keypathwayminer)等软件均可实现对差异表达基因的注释、富集分析和功能预测。通常认为,享有同一生物功能和调控机制的一组基因,具有相似的GO注释术语,只要算出已知基因和候选基因的GO术语相似度,就能推断出候选基因产物的近似功能。与已知功能基因GO语义的相似性比较,为未知基因功能预测提供了重要参考。尽管基因语义相似性已得到广泛研究,但具体应用时,基于GO注释两个基因产物相似性的界定仍然不清楚。为解决这个问题,2012年,王刚[71]基于OBO计划(Open biological and biomedical ontology foundry)(http://www.obofoundry.org), 采用一种公开的、正交的、实例化的语义格式,使本体之间的数据整合更加精确、界定更加清晰、实用性更强。DAVID(https://david.ncifcrf.gov/)、WebGestalt(http://www.webgestalt.org/option.php)、GO均是GO富集分析常用的数据库,可以对差异基因进行GO分类,并分析基于离散分布分类结果的显著性、错判率和富集度,得到与试验目的有显著联系的、低错判率的和高富集度的基因功能分类,再从该分类中找出同属于一个GO功能分类的变化基因,检验其统计学意义后得出变化基因主要参与的生物功能。

KEGG分析是通过计算差异表达基因通路的超几何分布概率,对每个差异表达基因通路上的所有差异表达基因进行富集,再经过统计学检验确定与某些差异表达基因显著相关的通路,进而寻找差异表达基因可能参与的细胞生化过程。通路富集分析的生物学假设是:通路上游基因表达改变,会导致下游相关基因表达改变,从而改变通路中大量基因的表达,且表达量的变化达到了富集分析的统计学显著水平。在众多差异表达基因中,很多基因在相应通路中并不是直接的相互调控关系,而是共同参与某个过程的不同环节,这些基因粗略地构成了该通路的整体轮廓。分析时,把差异表达基因集导入通路分析软件,通过超几何分布概率计算和基因差异表达显著性检验,预测差异表达基因显著富集的通路。比如,通过比较物种间转录组数据结合富集分析,成功找出了影响B. dorsalisB. correcta翅发育的基因EGFR[72],并验证了不同物种中α和β胰腺细胞亚型的转换方式[73]。当然,预测结果还需深入观察、理解某个核心通路中基因的相互作用,才能判断其中的差异表达基因是否有生物学意义,其中的某些不显著通路也值得从功能注释角度进行解读,只要结果可以解释,有生物学意义,也可以作为候选基因进行后续生物学验证分析。KEGG等数据库收录的是已有的研究结果,其中的很多通路信息远没有完善,致使许多通路只能显示大概的调控途径,而其间还有什么转录因子参与、是否还有其他代谢物生成都不知道。数据库中这些通路的完整性,也影响到pathway富集分析结果,随着研究的深入,每个通路都会有越来越多的节点被揭示出来,生命奥秘的面纱将渐渐被揭开,后续的差异基因功能分析也将越来越快捷、精准、高效。

2.5 利用BSR-Seq法定位并挖掘基因

集群分离分析(bulked segregant analysis, BSA)是将目标性状的两个极端差异表型亲本杂交,F1代自交后得到性状分离的F2代,选取F2代中表型极端差异个体组织样品分别混合构建两个DNA或mRNA池,两池间的基因序列差异片段即为所关注的基因或数量性状基因座(quantitative trait locus, QTL)可能存在的候选区域[74-75]。该方法主要适用于畜禽质量性状单基因或数量性状主效基因的定位,是快速获得目标性状主效基因或与其紧密连锁分子标记的有效方法,但对微效基因定位分析意义不大。BSR-Seq方法将BSA与RNA-Seq结合起来,实现了试验设计、测序分析、差异表达分析、基因功能分析与鉴定的全链条设计,该方法首先选择分离群中性状极端表型个体,采集组织样构建两个差异表型样本池,分别提取总RNA,进行转录组测序,根据测序得到的碱基总量(bp)与物种基因组大小的比值确定测序乘数。将转录组测序得到的clean reads与参考基因组进行比对,将比对到参考序列唯一位置的reads用于SNP发掘。然后用经典贝叶斯算法分析SNP位点,最终找出与突变表型紧密相关的基因组区段[76]。大量的RNA-Seq数据确保了特定表达模式下,对特定畜禽性状表达差异相关基因的SNP标记开发,而通过进一步精细定位和相关基因表达差异分析,则可逐渐确定候选基因及其功能。但该方法的定位结果由畜种亲本多态性、测序深度、混池个数等多个参数决定,如想得到最佳参数,则需进行多次模拟试验及参考基因组的支持。

3 展望

转录组测序技术为基因挖掘及其功能鉴定搭建了转录组学研究平台。然而,随着测序数据的爆炸式增长,如何建立高效的生物学分析方法,深度挖掘海量测序数据的生物学信息,实现功能基因的高效鉴定与生命奥秘的深度解析,越来越成为本领域急需解决的技术问题。为此,转录组数据挖掘不但从统计分析方法上进行了深入探索,在研究策略上也从最初仅利用本物种转录组数据进行基因挖掘和功能基因鉴定,向充分利用不同物种同源基因表达的相似性和差异性方向发展,为畜禽重要经济性状遗传解析及重要功能基因挖掘鉴定提供了新思路。其他物种转录组信息的引入为同一相似生物过程重要功能基因挖掘和鉴定及畜禽进化与分化研究,开阔了研究思路,也提供了更广泛的证据。随着研究的不断拓展和深入,GenBank等数据库不断得到丰富,功能基因挖掘和鉴定的基础越来越扎实。而且,测序技术也在不断发展,悄然兴起的第三代单分子测序技术,拥有比二代测序技术更快的速度、更低的价格和更高的精度,并能够测得更长的序列片段[77]。同时,光学图谱技术可基于单个DNA分子生成高分辨率、有序、完整基因组的限制性内切酶酶切图谱[78],可辅助保证序列拼接组装结果的准确性和真实性,加快畜禽重要经济性状的网络解析及重要功能基因挖掘鉴定。伴随着生命奥秘探索的逐渐深入,重复序列识别难题也将逐渐得到解决,海量测序数据的生物学意义必将会得到更高效、深度的揭示。

参考文献
[1] VELCULESCU V E, ZHANG L, ZHOU W, et al. Characterization of the yeast transcriptome[J]. Cell, 1997, 88(2): 243–251.
[2] DAL MOLIN A, BARUZZO G, DI CAMILLO B. Single-cell RNA-sequencing:assessment of differential expression analysis methods[J]. Front Genet, 2017, 8: 62. DOI: 10.3389/fgene.2017.00062
[3] 吴良涛, 郑敏, 华敏, 等. 基于RNA-Seq筛选鸭肠炎病毒感染鸭脾差异表达免疫相关基因[J]. 畜牧兽医学报, 2017, 48(2): 297–306.
WU L T, ZHENG M, HUA M, et al. Screening of differentially expressed immune-related genes from duck spleen with duck enteritis virus infection based on RNA-Seq technology[J]. Acta Veterinaria et Zootechnica Sinica, 2017, 48(2): 297–306. (in Chinese)
[4] 付言锋, 李兰, 赵为民, 等. 转录组测序分析猪胚胎附植的中期和后期卵巢差异表达基因[J]. 畜牧兽医学报, 2018, 49(9): 1879–1888.
FU Y F, LI L, ZHAO W M, et al. Differentially expressed genes analysis between mid-implantation and post-implantation in porcine ovary by RNA sequencing[J]. Acta Veterinaria et Zootechnica Sinica, 2018, 49(9): 1879–1888. (in Chinese)
[5] 浦思颖, 郑杰, 杨远潇, 等. 牦牛新鲜囊胚与玻璃化冻融囊胚转录组的比较分析[J]. 畜牧兽医学报, 2018, 49(4): 709–717.
PU S Y, ZHENG J, YANG Y X, et al. Comparative transcriptome analysis between fresh and vitrified-thawed blastocysts of the yak (Bos grunniens)[J]. Acta Veterinaria et Zootechnica Sinica, 2018, 49(4): 709–717. (in Chinese)
[6] HAAS B J, ZODY M C. Advancing RNA-Seq analysis[J]. Nat Biotechnol, 2010, 28(5): 421–423. DOI: 10.1038/nbt0510-421
[7] GAUDET P, ŠKUNCA N, HU J C, et al. Primer on the gene ontology[J]. Methods Mol Biol, 2017, 1446: 25–37. DOI: 10.1007/978-1-4939-3743-1
[8] KANEHISA M, FURUMICHI M, TANABE M, et al. KEGG:new perspectives on genomes, pathways, diseases and drugs[J]. Nucleic Acids Res, 2017, 45(D1): D353–D361. DOI: 10.1093/nar/gkw1092
[9] LI J J, HUANG H Y, BICKEL P J, et al. Comparison of D. melanogaster and C. elegans developmental stages, tissues, and cells by modENCODE RNA-Seq data[J]. Genome Res, 2014, 24(7): 1086–1101. DOI: 10.1101/gr.170100.113
[10] 陈巍.人与小鼠脑组织转录组基因表达与功能比较研究及多巴胺D5受体与胃泌素受体在促尿钠排泄中的协同互作研究[D].北京: 北京协和医学院, 2016.
CHEN W.Comparative analysis of gene expression and function between human and mouse brain tissues and the synergistic roles of cholecystokinin B and dopamine D5 receptors on the regulation of renal sodium excretion[D].Beijing: Peking Union Medical College, 2016.(in Chinese)
[11] BOLGER A M, LOHSE M, USADEL B. Trimmomatic:a flexible trimmer for Illumina sequence data[J]. Bioinformatics, 2014, 30(15): 2114–2120. DOI: 10.1093/bioinformatics/btu170
[12] BROWN J, PIRRUNG M, MCCUE L A. FQC Dashboard:integrates FastQC results into a web-based, interactive, and exrensible FASTQ quality control tool[J]. Bioinformatics, 2017, 33(19): 3137–3139. DOI: 10.1093/bioinformatics/btx373
[13] LANGMEAD B.Aligning short sequencing reads with Bowtie[J].Curr Protoc Bioinform, 2010, chapter 11(Unit11.17): Unit11.7.
[14] LANGMEAD B, SALZBERG S L. Fast gapped-read alignment with Bowtie 2[J]. Nat Methods, 2012, 9(4): 357–359. DOI: 10.1038/nmeth.1923
[15] KIM D, PERTEA G, TRAPNELL C, et al. TopHat2:accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions[J]. Genome Biol, 2013, 14(4): R36.
[16] DOBIN A, DAVIS C A, SCHLESINGER F, et al. STAR:ultrafast universal RNA-seq aligner[J]. Bioinformatics, 2013, 29(1): 15–21. DOI: 10.1093/bioinformatics/bts635
[17] WU Q, SONG J Y, SUN Y Q, et al. Transcript profiles of Panax quinquefolius from flower, leaf and root bring new insights into genes related to ginsenosides biosynthesis and transcriptional regulation[J]. Physiol Plant, 2010, 138(2): 134–149.
[18] LI R Q, YU C, LI Y R, et al. SOAP2:an improved ultrafast tool for short read alignment[J]. Bioinformatics, 2009, 25(15): 1966–1967. DOI: 10.1093/bioinformatics/btp336
[19] ZYTNICKI M, QUESNEVILLE H. S-MART, a software toolbox to aid RNA-seq data analysis[J]. PLoS One, 2011, 6(10): e25988. DOI: 10.1371/journal.pone.0025988
[20] GRABHERR M G, HAAS B J, YASSOUR M, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome[J]. Nat Biotechnol, 2011, 29(7): 644–652. DOI: 10.1038/nbt.1883
[21] HAAS B J, PAPANICOLAOU A, YASSOUR M, et al. De novo transcript sequence reconstruction from RNA-Seq using the Trinity platform for reference generation and analysis[J]. Nat Protoc, 2013, 8(8): 1494–1512. DOI: 10.1038/nprot.2013.084
[22] XU Z X, JIE H, CHEN B L, et al. Illumina-based de novo transcriptome sequencing and analysis of Chinese forest musk deer[J]. J Genet, 2017, 96(6): 1033–1040. DOI: 10.1007/s12041-017-0872-x
[23] WEINER I, ATAR S, SCHWEITZER S, et al. Enhancing heterologous expression in Chlamydomonas reinhardtii by transcript sequence optimization[J]. Plant J, 2018, 94(1): 22–31. DOI: 10.1111/tpj.2018.94.issue-1
[24] HU J F, BORITZ E, WYLIE W, et al. Stochastic principles governing alternative splicing of RNA[J]. PLoS Comput Biol, 2017, 13(9): e1005761. DOI: 10.1371/journal.pcbi.1005761
[25] PEVZNER P A, TANG H, WATERMAN M S. An eulerian path approach to DNA fragment assembly[J]. Proc Natl Acad Sci U S A, 2001, 98(17): 9748–9753. DOI: 10.1073/pnas.171285098
[26] MARDIS E R. Advances in genome biology and technology[J]. Pharmacogenomics, 2004, 5(4): 355–356. DOI: 10.1517/14622416.5.4.355
[27] 柳军涛.基于高通量RNA-Seq数据转录组拼接的算法研究[D].济南: 山东大学, 2017.
LIU J T.Algorithm studies of transcriptome assembly based on high throughput RNA-Seq data[D].Ji'nan: Shandong University, 2017. (in Chinese)
[28] 王春宇.生物高通量测序片段拼接与分子标记识别算法研究[D].哈尔滨: 哈尔滨工业大学, 2015.
WANG C Y.Research on biological high-throughput sequencing fragment assembly and molecular biomarker detection algorthms[D]. Harbin: Harbin Institute of Technology, 2015.(in Chinese)
[29] TRAPNELL C, ROBERTS A, GOFF L, et al. Differential gene and transcript expression analysis of RNA-Seq experiments with TopHat and Cufflinks[J]. Nat Protoc, 2012, 7(3): 562–578. DOI: 10.1038/nprot.2012.016
[30] LI H, HANDSAKER B, WYSOKER A, et al. The sequence alignment/map format and SAMtools[J]. Bioinformatics, 2009, 25(16): 2078–2079. DOI: 10.1093/bioinformatics/btp352
[31] LI B, DEWEY C N. RSEM:accurate transcript quantification from RNA-Seq data with or without a reference genome[J]. BMC Bioinformatics, 2011, 12: 323. DOI: 10.1186/1471-2105-12-323
[32] LOVE M I, HUBER W, ANDERS S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2[J]. Genome Biol, 2014, 15(12): 550. DOI: 10.1186/s13059-014-0550-8
[33] ROBINSON M D, MCCARTHY D J, SMYTH G K. edgeR:a Bioconductor package for differential expression analysis of digital gene expression data[J]. Bioinformatics, 2010, 26(1): 139–140. DOI: 10.1093/bioinformatics/btp616
[34] DJEBALI S, WUCHER V, FOISSAC S, et al.Bioinformatics pipeline for transcriptome sequencing analysis[M]//ØROM U. Enhancer RNAs. Methods in Molecular Biology, vol 1468.New York: Humana Press, 2017: 201-219.
[35] WANG L K, FENG Z X, WANG X, et al. DEGseq:an R package for identifying differentially expressed genes from RNA-seq data[J]. Bioinformatics, 2010, 26(1): 136–138. DOI: 10.1093/bioinformatics/btp612
[36] SHIN J Y, CHOI S H, CHOI D W, et al.Differential gene expression of RNA-Seq experiments of primo vessel in rabbit lymph[J].J Acupunct Meridian Stud, 2018, doi: 10.1016/j.jams.2018.10.008.
[37] SAMMETH M, FOISSAC S, GUIGÓ R. A general definition and nomenclature for alternative splicing events[J]. PLoS Comput Biol, 2008, 4(8): e1000147. DOI: 10.1371/journal.pcbi.1000147
[38] MOONEY M A, NIGG J T, MCWEENEY S K, et al. Functional and genomic context in pathway analysis of GWAS data[J]. Trends Genet, 2014, 30(9): 390–400. DOI: 10.1016/j.tig.2014.07.004
[39] LUO Y J, QIU Y, NA R, et al. A Golden gate and gateway double-compatible vector system for high throughput functional analysis of genes[J]. Plant Sci, 2018, 271: 117–126. DOI: 10.1016/j.plantsci.2018.03.023
[40] FABREGAT A, KORNINGER F, VITERI G, et al. Reactome graph database:Efficient access to complex pathway data[J]. PLoS Comput Biol, 2018, 14(1): e1005968. DOI: 10.1371/journal.pcbi.1005968
[41] NEAVES S R, TSOKA S, MILLARD L A C. Reactome Pengine:A web-logic API to the Homo sapiens reactome[J]. Bioinformatics, 2018, 34(16): 2856–2858. DOI: 10.1093/bioinformatics/bty181
[42] AI C, KONG L. CGPS:A machine learning-based approach integrating multiple gene set analysis tools for better prioritization of biologically relevant pathways[J]. J Genet Genomics, 2018, 45(9): 489–504. DOI: 10.1016/j.jgg.2018.08.002
[43] SUBRAMANIAN A, TAMAYO P, MOOTHA V K, et al. Gene set enrichment analysis:a knowledge-based approach for interpreting genome-wide expression profiles[J]. Proc Natl Acad Sci U S A, 2005, 102(43): 15545–15550. DOI: 10.1073/pnas.0506580102
[44] PITA-JUÁREZ Y, ALTSCHULER G, KARIOTIS S, et al. The pathway coexpression network:revealing pathway relationships[J]. PLoS Comput Biol, 2018, 14(3): e1006042. DOI: 10.1371/journal.pcbi.1006042
[45] MUÑOZ A G, KUTMON M, EIJSSEN L, et al. Pathway analysis of transcriptomic data shows immunometabolic effects of vitamin D[J]. J Mol Endocrinol, 2018, 60(2): 95–108. DOI: 10.1530/JME-17-0186
[46] YANG H L, HUO P F, HU G Z, et al. Identification of gene markers associated with metastasis in clear cell renal cell carcinoma[J]. Oncol Lett, 2017, 13(6): 4755–4761. DOI: 10.3892/ol.2017.6084
[47] XIE C, MAO X Z, HUANG J J, et al. KOBAS 2.0:a web server for annotation and identification of enriched pathways and diseases[J]. Nucleic Acids Res, 2011, 39(S2): W316–W322.
[48] YANG X, ZHU S, LI L, et al. Identification of differentially expressed genes and signaling pathways in ovarian cancer by integrated bioinformatics analysis[J]. Onco Targets Ther, 2018, 11: 1457–1474. DOI: 10.2147/OTT
[49] 郜金荣. 分子生物学[M]. 北京: 化学工业出版社, 2011.
GAO J R. Molecular biology[M]. Beijing: Chemical Industry Press, 2011. (in Chinese)
[50] WANG H J, KIRA Y, HAMURO A, et al. Differential gene expression of extracellular-matrix-related proteins in the vaginal apical compartment of women with pelvic organ prolapse[J]. Int Urogynecol J, 2018. DOI: 10.1007/s00192-018-3637-z
[51] SHANNON P, MARKIEL A, OZIER O, et al. Cytoscape:a software environment for integrated models of biomolecular interaction networks[J]. Genome Res, 2003, 13(11): 2498–2504. DOI: 10.1101/gr.1239303
[52] HAO M, JI X R, CHEN H, et al. Cell cycle and complement inhibitors may be specific for treatment of spinal cord injury in aged and young mice:transcriptomic analyses[J]. Neural Regen Res, 2018, 13(3): 518–527. DOI: 10.4103/1673-5374.226405
[53] LI Z J, LIU X L, ZHANG P P, et al. Comparative transcriptome analysis of hypothalamus-regulated feed intake induced by exogenous visfatin in chicks[J]. BMC Genomics, 2018, 19(1): 249. DOI: 10.1186/s12864-018-4644-7
[54] CARTER S L, BRECHBUHLER C M, GRIFFIN M, et al. Gene co-expression network topology provides a framework for molecular characterization of cellular state[J]. Bioinformatics, 2004, 20(14): 2242–2250. DOI: 10.1093/bioinformatics/bth234
[55] 娄鹏博.加权基因共表达网络分析(WGCNA)筛选猪卵母细胞脂肪沉积相关基因[D].成都: 四川农业大学, 2016.
LOU P B.Weighted gene coexpression network analysis (WGCNA) identify the fat deposition gene in pig oocytes[D]. Chengdu: Sichuan Agricultural University, 2016.(in Chinese)
[56] JING R X, GU L T, LI J Q, et al. A transcriptomic comparison of theca and granulosa cells in chicken and cattle follicles reveals ESR2 as a potential regulator of CYP19A1 expression in the theca cells of chicken follicles[J]. Comp Biochem Physiol Part D Genomics Proteomics, 2018, 27: 40–53. DOI: 10.1016/j.cbd.2018.04.002
[57] 苏蕊, 乔贤, 范一星, 等. 基于RNA-Seq绒山羊皮肤两种毛囊差异表达基因的分析[J]. 家畜生态学报, 2017, 38(11): 15–20.
SU R, QIAO X, FAN Y X, et al. Analysis on differentially expressed genes by RNA-Seq in two kinds of hair follicle in cashmere goat's skin[J]. Acta Ecologae Animalis Domastici, 2017, 38(11): 15–20. DOI: 10.3969/j.issn.1673-1182.2017.11.003 (in Chinese)
[58] 王旭平, 杨胜林, 柳诚刚, 等. 番鸭下丘脑组织的RNA-Seq特征分析[J]. 中国畜牧兽医, 2018, 45(3): 604–611.
WANG X P, YANG S L, LIU C G, et al. Characteristics analysis of hypothalamus of Muscovy ducks by RNA-Seq[J]. China Animal Husbandry & Veterinary Medicine, 2018, 45(3): 604–611. (in Chinese)
[59] WANG Q Q, LU W K, YANG J K, et al. Comparative transcriptomics in three Passerida species provides insights into the evolution of avian mitochondrial complex Ⅰ[J]. Comp Biochem Physiol Part D Genomics Proteomics, 2018, 28: 27–36. DOI: 10.1016/j.cbd.2018.06.002
[60] KUMAR S, STECHER G, TAMURA K. MEGA7:Molecular evolutionary genetics analysis version 7.0 for bigger datasets[J]. Mol Biol Evol, 2016, 33(7): 1870–1874. DOI: 10.1093/molbev/msw054
[61] BODENHOFER U, BONATESTA E, HOREJŠ-KAINRATH C, et al. msa:an R package for multiple sequence alignment[J]. Bioinformatics, 2015, 31(24): 3997–3999.
[62] NAKAMURA T, YAMADA K D, TOMⅡ K, et al. Parallelization of MAFFT for large-scale multiple sequence alignments[J]. Bioinformatics, 2018, 34(14): 2490–2492. DOI: 10.1093/bioinformatics/bty121
[63] 谷俊峰, 王希诚, 赵金城. 多序列渐进比对算法及其改进算法的研究与比较[J]. 大连大学学报, 2005, 26(2): 48–53.
GU J F, WANG X C, ZHAO J C. Study and comparison of multiple sequence progressive alignment algorithm and its reformative algorithms[J]. Journal of Dalian University, 2005, 26(2): 48–53. DOI: 10.3969/j.issn.1008-2395.2005.02.015 (in Chinese)
[64] NIELSEN R, HELLMANN I, HUBISZ M, et al. Recent and ongoing selection in the human genome[J]. Nat Rev Genet, 2007, 8(11): 857–868.
[65] AL-TOBASEI R, ALI A, LEEDS T D, et al. Identification of SNPs associated with muscle yield and quality traits using allelic-imbalance analyses of pooled RNA-Seq samples in rainbow trout[J]. BMC Genomics, 2017, 18(1): 582. DOI: 10.1186/s12864-017-3992-z
[66] VAN BELLEGHEM S M, ROELOFS D, VAN HOUDT J, et al. De novo transcriptome assembly and SNP discovery in the wing polymorphic salt marsh beetle Pogonus chalceus (Coleoptera, Carabidae)[J]. PLoS One, 2012, 7(8): e42605. DOI: 10.1371/journal.pone.0042605
[67] 何骋.转录组测序数据分析在玉米籽粒功能基因挖掘中的应用[D].北京: 中国农业大学, 2017.
HE C.Analysis of RNA-Seq data to identify functional genes in maize kernel[D].Beijing: China Agricultural University, 2017.(in Chinese)
[68] 陈敏惠.牛和猪基因组选择信号研究[D].北京: 中国农业大学, 2016.
CHEN M H.Selection signatures in the genomes of cattle and pigs[D].Beijing: China Agricultural University, 2016.(in Chinese)
[69] JANSSEN S F, GORGELS T G M F, BOSSERS K, et al. Gene expression and functional annotation of the human ciliary body epithelia[J]. PLoS One, 2012, 7(9): e44973. DOI: 10.1371/journal.pone.0044973
[70] BOOIJ J C, VAN SOEST S, SWAGEMAKERS S M A, et al. Functional annotation of the human retinal pigment epithelium transcriptome[J]. BMC Genomics, 2009, 10: 164. DOI: 10.1186/1471-2164-10-164
[71] 王刚.基于疾病表型的基因语义相似性分析与应用[D].西安: 西安电子科技大学, 2012.
WANG G.Analysis and application of the gene semantic similarity based on disease phenotype[D].Xi'an: Xidian University, 2012.(in Chinese)
[72] GUO S K, ZHAO Z H, LIU L J, et al. Comparative transcriptome analyses uncover key candidate genes mediating flight capacity in Bactrocera dorsalis (Hendel) and Bactrocera correcta (Bezzi) (Diptera:Tephritidae)[J]. Int J Mol Sci, 2018, 19(2): 396. DOI: 10.3390/ijms19020396
[73] TARIFEÑO-SALDIVIA E, LAVERGNE A, BERNARD A, et al. Transcriptome analysis of pancreatic cells across distant species highlights novel important regulator genes[J]. BMC Biol, 2017, 15(1): 21. DOI: 10.1186/s12915-017-0362-x
[74] 陈浣, 夏菲, 吴新儒, 等. 集群分离分析法在植物基因定位上的应用[J]. 基因组学与应用生物学, 2016, 35(6): 1546–1551.
CHEN H, XIA F, WU X R, et al. The application of Bulked Segregant analysis in plant gene mapping[J]. Genom Appl Biol, 2016, 35(6): 1546–1551. (in Chinese)
[75] HOU X G, GUO Q, WEI W Q, et al. Screening of genes related to early and late flowering in tree peony based on Bulked Segregant RNA sequencing and verification by quantitative real-time PCR[J]. Molecules, 2018, 23(3): 689. DOI: 10.3390/molecules23030689
[76] LIU S Z, YEH C T, TANG H M, et al. Gene mapping via Bulked Segregant RNA-Seq (BSR-Seq)[J]. PLoS One, 2012, 7(5): e36406. DOI: 10.1371/journal.pone.0036406
[77] EDGER P P, VANBUREN R, COLLE M, et al. Single-molecule sequencing and optical mapping yields an improved genome of woodland strawberry (Fragaria vesca) with chromosome-scale contiguity[J]. Gigascience, 2018, 7(2): 1–7.
[78] UDALL J A, DAWE R K. Is it ordered correctly? Validating genome assemblies by optical mapping[J]. Plant Cell, 2018, 30(1): 7–14. DOI: 10.1105/tpc.17.00514