黄芩高通量转录组测序数据组装和分析

王德富，成媛媛，杨 锋，牛二波，牛颜冰

（山西农业大学生命科学学院，山西太谷030801）

黄芩高通量转录组测序数据组装和分析

王德富，成媛媛，杨 锋，牛二波，牛颜冰

（山西农业大学生命科学学院，山西太谷030801）

采用高通量测序技术对黄芩的转录组进行测序，获得了29 099 899条reads数据，拼接后得到53 353条Unigene；将所获得的Unigene与COG，GO，KEGG，Swiss-Prot，NR这5个公共数据库进行比对，结果发现，分别有10 756，21 950，8 101，20 339，29 288条Unigene可比对到以上5个数据库中；已注释的Unigene与COG数据库比对后按功能可分为25类；根据GO功能可分为三大类57个亚类；经过与KEGG数据库比对后按照代谢通路可分为116类；利用GetORF软件进行ORF预测，获得长度大于300 nt的ORF共20 552个；通过SSR分析，共获得5 658个SSR标记。获得的转录组信息可为今后进行黄芩分子标记的开发和关键基因的克隆及功能分析等研究提供基础数据。

黄芩；转录组测序；生物信息学分析

黄芩（Scutellaria baicalensis Georgi）为唇形科黄芩属植物的干燥根[1-3]，具有清热燥湿、凉血安胎、解毒等功效，常用于治疗胸闷呕恶、湿热痞满、泻痢、胎动不安等症[4-6]，临床疗效极为广泛。随着黄芩在药物开发中的广泛应用，其需求量逐年上升，为了满足日益膨胀的市场需求和减轻对黄芩野生资源的掠夺式采挖压力，采用生物技术大规模合成有效成分，已成为药物生产和新药开发的重要途径。但由于黄芩目前还未进行全基因组测序，相关基因资料不全，导致黄芩中有效成分生物合成途径方面的研究比较滞后。而高通量转录组测序技术的出现，为研究黄芩活性成分生物代谢途径提供了重要的基因资源，并为开展黄芩功能基因组学研究提供了全新的思路和方法[7-8]。

转录组分析可在不知道物种基因信息的情况下，直接从RNA水平分析细胞代谢规律，揭示基因表达与一些生命现象之间的关系，进而对细胞进行修饰和改造[9-11]。该法特别适合基因信息十分匮乏的物种[12-14]，并已在人类[15]、酵母[16]、玉米[17]、拟南芥[18]、大豆[19]、蝴蝶[20]等多个物种中得到广泛应用。在中药材方面，郝大程等[21]、李铁柱等[22]、李滢等[23]、郭溆[24]分别对虎杖根、杜仲果实和叶片、丹参根、铁皮石斛和人参进行了转录组测序，并发现了一些参与药效成分生物合成相关的转录本和关键酶序列。这些研究结果说明，RNA-Seq技术在药用植物研究中也得到了极其广泛的应用。

本研究利用Illumina HiSeqTM2500对黄芩进行转录组测序和分析，从黄芩转录组中识别出与黄芩苷生物合成相关的关键基因，并发现了大量的简单重复序列（SSR）。这些重要基因和SSR序列的发现，为进一步克隆其全长、研究其功能和进行分子标记的开发提供了基础数据，同时也为药效成分的生物合成研究奠定了基础。

1 材料和方法

1.1 材料

2014年8月在山西绛县2年生黄芩种植基地，分别采集最新长出的黄芩根、茎、叶和花，每个组织选取10株取样，用蒸馏水冲洗干净后，置于液氮中冷冻保存，备用。

1.2 方法

1.2.1 总RNA提取及检测 使用TRIzol（Invitrogen）法提取黄芩4种不同组织的总RNA，并测定其质量。若RIN（RNA完整值）大于8.5，28S/18S大于2，OD260/OD280在1.8～2.2，OD260/OD230大于1.8，则样品合格，然后将质量较好的RNA样本等比例混合后进行转录组测序。

1.2.2 黄芩cDNA文库的建立和测序 cDNA文库构建和转录组测序工作委托北京百迈客生物科技有限公司完成。

1.2.3 黄芩转录组数据处理和组装 将测序得到的原始图像数据转化为原始读序（rawreads），去除其中的接头序列和低质量的reads，即可获得高质量的干净读序（clean reads）。利用Trinity[25]软件对clean reads进行de novo组装。

1.2.4 功能注释和分类 首先通过BLASTX将得到的Unigene分别与核酸和蛋白数据库进行比对（E-value＜10-5），获得最佳注释信息；接着用Blast2GO和WEGO软件对所有GO注释信息进行GO功能分类统计[26]；用COG数据库进行Unigene功能注释和分类，了解基因功能分布特性；用KEGG数据库对Unigene进行功能注释，分析某个Unigene在细胞代谢通路中的功能定位[27]。

1.2.5 ORF预测与SSR位点搜索 利用GetORF软件对Unigene进行ORF（Open ReadingFrame）预测，选取预测结果中最长的序列作为该Unigene最终的ORF。使用MISA（http://pgrc.ipk-gatersleben. de/misa/misa.htmlt/）软件对筛选得到的1 kb以上的Unigenes做SSR序列分析，检测Unigene序列中SSR的分布特征。

2 结果与分析

2.1 RNA质量检测

利用Agilent 2100生物分析仪（Agilent Technology，USA）和NanoDrop分光光度计（Thermo，USA）对提取的黄芩RNA进行检测，结果表明，其RIN值大于8.5，OD260/OD280在1.8～2.2，OD260/OD230大于1.8，28S/18S大于2，说明所得到的RNA质量较好，能够满足后续的转录组测序试验（表1）。

表1 RNA质量检测结果

2.2 黄芩测序数据质量评价及序列组装分析

经HiSeqTM2500高通量测序获得29 099 899条读序（clean data），总核苷酸数达7 326 522 924个，数据量为7.33 G，GC含量为45.87%，Q30（测序错误率＜1%）值为90.743%（图1-a）；同时碱基分布分析表明，除了由于在read 5′端前十几个碱基存在明显的偏向性而导致前端波动较大外，每个测序循环中AT和GC的含量在整个测序过程中基本呈水平线，含量稳定不变（图1-b），说明所得到的测序数据质量较好，可进行后续分析。

利用Trinity软件对上述获得的clean reads进行组装，得到4538775条contigs，平均长149.51 bp，N50为148 bp；然后再通过contigs之间的相似性以及双末端信息将其组装得到107 533条transcripts，平均长2 038 bp，N50为1 257.69 bp；进一步的去冗组装获得53 353条Unigenes，平均长1 467 bp，N50为797.64 bp（表2）。所得到transcript和Unigene的N50分别达到1 257.69，797.64 bp，表明组装片段的完整性比较高。

表2 黄芩转录组数据组装统计结果

2.3 Unigene的功能注释

表3 Unigenes在5个数据库中的分布情况

将组装得到的53 353条Unigene与已知的COG，GO，KEGG，Swiss-Prot和NR这5个公共数据库进行序列比对。从表3可以看出，Unigene在5个不同数据库中的同源比对数目不同，分别为10 756，21 950，8 101，20 339，29 288条。总共所获得的同源比对信息数量为29 382条，占Unigene总数量的55.07%，但仍有23 971条（44.93%）序列定位不清楚，这可能与公共数据库中没有黄芩全基因组序列有关。

从图2-A可以看出，比对到NR数据库中的29 288条Unigenes，有40%的E值分布在10-50～10-5，20%位于 10-100～10-50，12%位于 10-150～10-100，小于10-150的占5%。从匹配的物种来源分析（图2-B），发现注释到番茄中Unigene占26%，注释到葡萄中Unigene占25%，注释到其他物种中的Unigene相对较少，如注释到大豆和麻风树中的Unigene均为3%，蓖麻、野草莓和樱桃的Unigene均为4%，毛果杨为6%，黄瓜为2%，其余23%注释到其他物种。

2.4 Unigene的功能分类

由表3可知，有21 950条Unigene（41.14%）获得GO数据库注释。这些Unigene被分在上述3大类的57个亚类中（图3），其中，生物过程分为25个亚类，获得的注释信息最多，共有65749条Unigene，占全部注释信息的44.33%，这里涉及代谢过程、细胞过程和对刺激的应激效应的Unigene比较多，分别有14 564，13 416，6 291条；其次是细胞成分，分为16个亚类，共有57 243条Unigene，占全部注释信息的38.60%，涉及到细胞部分、细胞和细胞器的Unigene比较多，分别有13 930，13 802，11 606条；分子功能的注释信息最少，共有25 315条Unigene，占17.07%，这其中的催化活性与结合排在前2位，分别有11 123，10 106条Unigene。

另外，对获取的Unigene进行了COG数据库功能注释，并归入不同的分类中，共有15 159条Unigene归入了25种COG分类中（图4）。其中，一般功能预测类包含Unigene最多（2 817条），其次是转录（1 329条）、复制、重组和修复（1 313条）、信号转导机制（1 122条）、翻译、核糖体结构和生物合成（1 111条）；涉及细胞活性和细胞核结构的Unigene较少，分别仅有13，5条；没有发现细胞外结构的功能基因。另外，有453条Unigene未被注释上，该注释与其他物种注释结果比较相似。

2.5 Unigene代谢途径分析

将所获得的53 353条Unigene映射到KEGG数据库中，有8 101条（15.18%）Unigene得到注释，并归入116种KEGG代谢途径中。从表4可看出，涉及核糖体代谢途径的Unigene数量最多（410条/ 5.06%），其次是氧化磷酸化（287条/3.54%）、内质网蛋白加工（270条/3.33%）、剪接体（252条/3.11%）和嘌呤代谢（245条/3.02%）等代谢通路。但归入生物素代谢和油菜素内酯合成等中的Unigene较少，数量不超过10条。

2.6 ORF与SSR位点分析

为进一步解析黄芩转录组编码的基因信息，利用GetORF软件预测了已得到注释的Unigene，共得到53 028个ORF，其中，长度大于300 nt的Unigene最多，占38.81%（20 552条）（图5）。另外，本研究还利用MISA软件找出5 658个SSR标记位点，其中，双碱基型重复最为丰富（2 928个/51.75%），其他类型依次为：单碱基型（1 540个/27.22%）、三碱基型（1 127个/19.92%）、四碱基型（37个/0.65%），五碱基型和六碱基型分别为11，15个，各占0.19%和0.27%（表5）。黄芩SSR位点的分析，将为黄芩的遗传标记研究提供非常重要的物质资源和依据。

表4 黄芩转录组Unigene KEGG代谢途径分类

续表4

表5 SSR分析统计结果

3 讨论

目前，RNA-seq技术已被广泛应用于各物种的转录组分析中[15，11，17，20-21]。本试验采用高通量测序技术对黄芩转录组进行测序，经过数据组装共得到53 353条Unigenes，平均序列长度为797.64 bp，其中，长度在1 kb以上的Unigenes有13 439条。经过生物信息学分析，成功获得功能注释的Unigene有29 382条（55.07%），但仍有23 971条（44.93%）Unigene未能获得注释。分析其原因可能有2个：一是本试验得到的Unigene序列长度较短，超过1/2的长度分布在500 nt以下，难以获得同源性比对，这在一定程度上增加了基因功能注释的难度；二是基因数据库中的生物信息暂时缺乏，一些表达不丰富的基因可能无法获得准确的功能注释。本试验所获得的黄芩转录组数据，可丰富黄芩的基因组信息和公共数据库信息，为进一步挖掘和鉴定出更多的黄芩功能基因奠定基础。

本研究利用生物信息学手段对所有获得的53 353条Unigenes进行了功能分类和代谢途径分析。有10 756条Unigenes得到COG功能注释，并归入25种COG分类中，但仍有453条Unigenes未得到功能注释，这可能与黄芩的基因组信息不完整有关。随着更多物种基因组和转录组测序结果的出现，使公共数据库不断得到补充和完善，这些功能未知的基因将被准确注释，并可能从中发现一些新的功能基因。另外，通过KEGG数据库分析，发现了一些与黄芩有效成分合成有关的关键酶序列，这不仅为基因的克隆和功能研究提供了一定的基础数据，还为进一步研究该有效成分生物合成途径中关键基因的调控机制奠定基础，同时也为应用生物学技术来合成黄芩中的有效成分提供可行性。

高通量测序技术与传统测序方法相比，操作简单，效率高，能够挖掘出大量的SSR资源。本研究利用MISA软件查找黄芩转录组测序数据，共发现了5 658个SSR位点，而且发现黄芩转录组SSR重复碱基以双碱基型重复最多，占所有SSR的51.75%，这与其他研究结果类似[28-29]。黄芩SSR位点的发现可为黄芩分子标记的开发、群体遗传多样性分析、遗传连锁图谱构建等后续试验研究奠定理论基础。本试验对黄芩转录组进行了初步的探究，弥补了黄芩基因组信息十分缺乏的局面，为将来进行黄芩分子生物学方面的相关研究打下基础，但要想深入了解黄芩的更多信息，今后仍需要加快黄芩转录组学和基因组学的相关研究步伐。

［1］李锡文，Hedge I C.中国植物志[M].北京：科学出版社&圣路易斯：密苏里植物园出版社，1994.

［2］黄爽.神农本草经[M].北京：中医古籍出版社，1982.

［3］中华人民共和国药典委员会.中华人民共和国药典 [M].北京：化学工业出版社，2005.

［4］中华本草编委.中华本草[M].上海：上海科学技术出版社，1998.

［5］Huang J M，Wang C J，Chou F P，et al.Protective effect of baicalin on tert-butyl hydroperoxide-induced rat hepatotoxicity[J].Arch Toxicol，2005，79（2）：102-109.

［6］Zhang D Y，Wu J，Ye F，et al.Inhibition of cancer cell proliferation and prostaglandin E2 synthesis byScutellaria baicalensis[J].Cancer Res，2003，63（14）：4037-4043.

［7］Marioni J C，Mason C E，Mane S M，et al.RNA-seq:An assessment oftechnical reproducibilityand comparison with gene expression arrays[J].Genome Res，2008，18（9）：1509-1517.

［8］Fullwood M J，Wei C L，Liu E T，et al.Next-generation DNA sequencing of paired-end tags（PET）for transcriptome and genome analyses[J].Genome Res，2009，19（4）：521-532.

［9］井赵斌，魏琳，俞靓，等.转录组测序及其在牧草基因资源发掘中的应用前景[J].草业科学，2011，28（7）：1364-1369.

［10］ Jewett M C，Oliveira A P，Patil K R，et al.The role of high-throughput transcriptome analysis in metabolic engineering [J].Biotechnol Bioproc Eng，2005，10：385-399.

［11］Donson J，Fang Y，Espiritu-Santo G，et al.Comprehensive gene expression analysis by transcript profiling[J].Plant Mol Biol，2002，48（1/2）：75-97.

［12］Sultan M，Schulz MH，Richard H，et al.A global viewof gene activity and alternative splicing by deep sequencing of the human transcriptome[J].Science，2008，321：956-960.

［13］Wang E T，Sandberg R，Luo S，et al.Alternative isoform regulation in human tissue transcriptomes[J].Nature，2008，456：470-476.

［14］Birzele F，Schaub J，Rust W，et al.Into the unknown:Expression profiling without genome sequence information in CHO by next generation sequencing[J].Nucleic Acids Res，2010，38（12）：3999-4010.

［15］Pan Q，Shai O，Lee L J，et al.Deep surveyingof alternative splicing complexity in the human transcriptome by high throughput sequencing[J].Nat Genet，2008，40（12）：1413-1415.

［16］Nagalakshmi U，Wang Z，Waern K，et al.The transcriptional landscape ofthe yeast genome defined byRNA sequencing[J].Science，2008，320：1344-1349.

［17］Kakumanu A，AmbavaramMM，Klumas C，et al.Effects of drought on gene expression in maize reproductive and leaf meristem tissue revealed byRNA-Seq[J].Plant Physiol，2012，160（2）：846-867.

［18］Filichkin S A，Priest H D，Givan S A，et al.Genome-wide mapping of alternative splicing in Arabidopsis thaliana[J].Genome Res，2010，20（1）：45-58.

［19］WangL，Cao C，Ma Q，et al.RNA-Seq analyses of multiple meristems of soybean:novel and alternative transcripts，evolutionary and functional implications[J].BMCPlant Biol，2014，14：169.

［20］Vera J C，Wheat C W，Fescemyer H W，et al.Rapid transcriptome characterization for a nonmodel organismusing454 pyrosequencing [J].Mol Ecol，2008，17（7）：1636-1647.

［21］郝大程，马培，穆军，等.中药植物虎杖根的高通量转录组测序及转录组特性分析 [J].中国科学：生命科学，2012，42（5）：398-412.

［22］李铁柱，杜红岩，刘慧敏，等.杜仲果实和叶片转录组数据组装及基因功能注释 [J].中南林业科技大学学报，2012，32（11）：122-130.

［23］李滢，孙超，罗红梅，等.基于高通量测序454 GS FLX的丹参转录组学研究[J].药学学报，2010，45（4）：524-529.

［24］郭溆.基于转录组测序的石斛生物碱和人参皂苷生物合成相关基因的发掘、克隆及鉴定[D].北京：北京协和医学院，2013.

［25］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.

［26］Götz S，García-Gómez J M，Terol J，et al.High-throughput functional annotation and data miningwith the Blast2GOsuite[J].Nucleic Acids Res，2008，36（10）：3420-3435.

［27］ Kanehisa M，Goto S.KEGG：Kyoto encyclopedia of genes and genomes[J].Nucleic Acids Res，2000，28：27-30.

［28］郑纪伟.柳树转录组高通量测序及SSR标记开发研究 [D].南京：南京林业大学，2013.

［29］陈浩东.达尔文氏棉旱胁迫转录组测序、EST-SSR开发及高密度遗传图谱构建[D].北京：中国农业科学院，2013.

Transcriptome Data Assembly and Analysis ofScutellaria baicalensisthrough High-throughput Sequencing

WANGDefu，CHENGYuanyuan，YANGFeng，NIUErbo，NIUYanbing
（College ofLife Sciences，Shanxi Agricultural University，Taigu 030801，China）

To study the transcriptome data of Scutellaria baicalensis Georgi，the root,stem,flower and leaf organs were used as the experimental material for the transcriptome sequencing,and analyzed by bioinformatics method.The transcriptome library of Scutellaria baicalensis contained 29 099 899 reads，and 53 353 Unigenes were obtained by assembling the Scaffolds in the transcriptome library; Unigene in the transcriptome of Scutellaria baicalensis could be divided into 25 classes according to the function by comparing Unigene and the COG database.The Unigene GO functions in the transcriptome library were classificated into 3 categories and 57 branches; Unigene in the transcriptome could be divided into116 classes accordingtothe metabolic pathway.Meanwhile,the paper alsogot number of 20 552 ORF using the GetORF software and got a number of 5 658 SSR markers by SSR analysis.Data presented in this study will constitute an important genomic resource for Scutellaria baicalensis Georgi and laya solid foundation for future gene clone and regulation research about the biosynthesis ofbaicalin.

Scutellaria baicalensis Georgi;transcriptome sequencing;bioinformatics analysis

S567.23＋9

A

1002-2481（2016）08-1065-08

10.3969/j.issn.1002-2481.2016.08.04

2016-07-22

国家自然科学基金项目（31540050）；山西农业大学科技创新基金项目（20142-10）；山西农业大学引进人才科研启动基金项目（2014YJ05）

王德富（1983-），男，甘肃庆阳人，讲师，博士，主要从事分子植物病毒学研究工作。牛颜冰为通信作者。