1. 引言
马铃薯(Solanum tuberosum L.)属茄科茄属作物，具有良好的适应性和均衡的营养，已经在世界约160个国家和地区广泛种植，是世界第四大粮食作物 [1]。我国马铃薯种植面积和产量均居世界第一，但单产量远不及世界发达国家水平，而病毒侵染引起的种薯退化是导致减产的重要原因之一 [2]。马铃薯病毒病会导致其块茎产量和品质的下降，长期以来给广大马铃薯产业工作者造成了巨大损失，也制约了我国产业发展。
病毒在细胞水平侵染马铃薯，因此一旦感病，很难用药物进行治疗，只能主要采取预防的办法，而预防病毒病的关键是对各种马铃薯病毒的快速检测。然而，侵染马铃薯的病毒多达40种以上，而且均为RNA病毒，用传统的酶联免疫(ELISA)和反转录聚合酶链式反应(RT-PCR)等方法需要投入大量的人力和时间，而且上述方法只能检测数十种主要病毒。
近年来，高通量测序发展迅速，已经被广泛应用于包括植物病毒生物学研究的众多领域 [3] [4] [5]。2011年，马铃薯基因组测序已基本完成，给许多马铃薯科研工作者带来了新的机遇。转录组(transcriptome)是指某种生理条件下细胞中产生的全部转录物(包括mRNA和非编码RNA)的总和 [6]。转录水平调控是生物体最重要的调控方式之一。近几年，转录组学研究是科学家研究的热门领域。通过高通量测序对样品的基因表达水平进行分析被称为RNA-seq，而在所测得的数据中不仅包含着样品的转录本，也含有那些感染寄主的RNA病毒序列信息，因此利用植物转录组测序数据，过滤掉寄主的序列，可以对病毒资源进行挖掘 [7]。本实验室在进行马铃薯低温糖化以及结薯机制的相关研究中已经进行了大量的转录组测序 [8] [9] [10]，利用这些现有的数据，可以快速而低成本的获得潜在的马铃薯病毒序列信息，从而“变废为宝”，为这些马铃薯病毒病的检测和防治工作奠定了重要的基础。
RNA-seq测序得到非常多很短的读段(reads)，想要获得病毒基因组序列，最核心的问题是将这些短的reads拼接组装成长的转录本，因此需要对RNA-seq数据进行从头组装(de novo assembly)。从头组装，是指没有参考基因组的帮助下，将reads拼接组装成真实DNA序列的过程。目前，基因组的从头组装主要基于以下两种算法：Overlap-Layout-Consensun (OLC)算法 [11] 和De Bruijn图算法 [12]。且后者更加适合reads较短、数据量较大的二代测序数据进行拼接组装 [13]。基于此算法开发的软件有：Trinity、SOAPdenovo、IDBA-UD等。因此，本研究利用一系列实验室相关研究获得的测序数据，在通过比对过滤掉马铃薯基因组数据后，用上述3种软件进行了病毒序列数据的组装，通过比较组装效果，确定从RNA-seq测序数据中发掘病毒序列的最佳分析软件。
2. 材料与方法
2.1. 数据收集
研究用到的相关测序数据来自本实验室其它研究积累的转录组测序数据，分别随机选择了6个PE125数据样本(chang-1, chang-2, CW2-1, CW2-2, duan-1, duan-2)和6个PE150数据样本(AC142, PVS, E20, E108, Ri-1, Ri-2)。
本研究使用的马铃薯参考基因组序列由国际马铃薯基因组测序协会(International Potato Genome Sequencing Consortium, PGSC)公布的对DM1-3 (S. phureja)进行全基因组测序的基因组序列及其注释文件，版本号为6.1 (http://solanaceae.plantbiology.msu.edu/pgsc_download.shtml)。相关病毒序列来自NCBI公共数据库下载的病毒标准序列数据库(ftp://ftp.ncbi.nih.gov/refseq/release/viral/) [14]。
2.2. 生物信息数据分析方法
2.2.1. 组装前数据预处理
在进行数据组装前，首先使用FastQC (v0.11.3)对数据进行质量评估，使用Trimmomatic [15] 进行质控，得到clean_data；接着，将clean_data数据比对到马铃薯参考基因组，以区分来源于马铃薯转录本的reads和非来源于马铃薯转录本的reads，即unmapped reads，作为进行组装方法比较的核心数据。
2.2.2. 使用Trinity软件进行组装
对于上一步得到的unmapped reads，即非来源于马铃薯转录本的数据，用Trinity (v2.1.0)软件进行重头组装：首先，对于未匹配到马铃薯基因组的数据，即上步得到的unmapped.bam文件，使用samtools (v0.1.18)根据bam文件的flag标签，提取双端reads，得到unmapped_reads1.fastq和unmapped_reads2.fastq。本研究所使用的命令如下：
samtools view unmapped.bam|awk '{if ($2~/69/) print >$1\n$10\n+\n$NF}' > unmapped_reads1.fastq
samtools view unmapped.bam|awk '{if ($2~/133/) print >$1\n$10\n+\n$NF}' > unmapped_reads2.fastq
得到双端reads后，使用Trinity (v2.1.0)进行重头组装。本研究所使用的命令如下：Trinity --seqType fq --max_memory 10G --left unmapped_reads1.fastq --right unmapped_reads2.fastq --CPU 10
2.2.3. 使用SOAPdenovo软件进行组装
使用SOAPdenovo (v2.04)进行组装，输入文件不能用上一步samtools提取的双端数据，因为samtools提取的双端reads是全部未匹配的reads，没有考虑reads的配对情况，会出现只有某一端，没有另一端而无法配对的情况。因此使用bam2fastq (v1.1.0)软件提取双端reads，该软件默认会去掉无法配对的reads；另外，SOAPdenovo组装需要建立一个配置文件，配置文件指定了输入文件的路径、数据特征和程序运行的参数。此处由于会设置average insert size参数，所以需要先用Picard (v1.119)计算average insert size，命令如下：java-jar/home/software/picard-tools-1.119/CollectInsertSizeMetrics.jar I = accepted_hits.bam O = insertsize.txt H = insertsize.pdf。最后，建好配置文件以后，可以运行SOAPdenovo的主程序。本研究使用命令如下(以k-mer取105为例)：SOAPdenovo-127mer all-s config.txt-o out_K105-K 105-p 10-R & > soapdenovo.log。
2.2.4. 使用IDBA-UD软件进行组装
IDBA-UD与SOAPdenovo都需要使用bam2fastq软件提取的双端数据；组装前，要先对双端数据进行合并，合并成一个文件：fq2fa--merge--filter unmapped_reads1.fastq unmapped_reads2.fastq unmapped_reads.fa；本研究使用命令如下(以125 bp reads为例)：idba_ud -r unmapped_reads.fa--mink 69--maxk 107--step 2 -o./denovo--num_threads 10。
2.2.5. 本地blast进行同源性比对
使用本地化blast+ (v2.2.31)将组装得到的contig与所构建的数据库III进行比对，研究组装结果中的病毒多样性，这些contig中可能会有已报道的对马铃薯有侵染能力的病毒，也可能会发现未报道的对马铃薯有侵染能力的病毒；再使用本地化blast+将组装得到的contig与所构建的数据库II进行比对，研究组装结果中contig具体与马铃薯病毒的哪个亚型或分离物同源性最高。
3. 结果与分析
3.1. 不同k-mer下软件的组装效果比较
本研究选择了contig N50和最长contig两个指标，作为衡量组装效果的指标，同时也作为选择最适k-mer的评判标准。N50反映的是一个组装结果整体的序列长度，N50越大，说明组装得到长序列越多，组装效果越好。另外，由于本研究进行的是病毒基因组组装，希望得到最长contig越长越好，因此本研究还选择最长contig作为最适k-mer选择的另一个指标。
本研究对SOAPdenovo和IDBA-UD两种软件在不同k-mer大小下，N50和max contig进行统计。如图1A和B所示，对于PE125的数据(Chang-1和CW2-1)，随k-mer从69增加到87-91左右，contig N50逐渐增大，k-mer值在91-95左右急剧下降，随后逐渐增大，因此Chang-1和CW2-1分别选择87和91为最佳的k-mer。而对于PE150的数据(AC142和PVS) (图1C和图1D)，随着k-mer逐渐增大，contig N50逐渐增大，在k-mer取127时达到最大值，因此最佳k-mer为127。但是，对于max contig这一指标，无论是PE125和PE150数据，其变化规律不明显(图1A-D)。
使用IDBA-UD软件进行从头组装，支持的最大k-mer值是123。无论是PE125还是PE150的数据，conitg N50、max contig与k-mer之间存在着相似的规律：对于PE125数据(Chang-1和CW2-1)，Contig N50在k-mer达到71时有一个明显上升，随后保存缓慢线性增长，k-mer达到105后达到最大值，随后有所回落。而max contig的长度也随k-mer值增大呈现一定的增长态势，以CW2-1为例(图1F)，k-mer值从93增加到95时，max contig长度增加300 bp左右，而Chang-1 (图1E)在k-mer值处于75~107之间，max contig保存不变。因此，为了获得最佳组装效果，该软件对PE125数据，k-mer选择105为宜。对于PE150数据(AC142和PVS) (图1G和图1H)，contig N50随k-mer的变化规律与PE125数据相似，只是一般在k-mer值121左右达到最大值。而max contig长度在contig N50极值范围内变不大。因此，对于PE150的reads推荐使用115~121作为最优k-mer的选择范围。
综上所述，对于SOAPdenovo软件，组装效果随着K-mer变化无明显变化规律，需要视实际组装情况选择合适K-mer。而对于IDBA-UD软件组装效果与k-mer取值的变化规律非常明显，即K-mer越大，组装效果越好。因此，IDBA-UD软件因此更容易选择最合适的k-mer值，而且组装效果更加容易让人信服。
3.2. 不同组装软件效果比较
本研究分别使用Trinity、SOAPdenovo、IDBA-UD软件对随机选择的PE125、PE150各6组数据进
Figure 1. Effect of the k-mer on the de novo assembly quality using SOAPdenovo and IDBA-UD
图1. SOAPdenovo和IDBA-UD软件k-mer值对组装质量的影响
行组装，比较组装效果，选择最适合作为病毒基因组组装的软件。为了衡量软件组装效果，选择了contig N50、contig平均长度、contig中位数长度、最长contig、大于1000的contig的数量、组装结果总contig数目以及大于1000 bp的contig数目占总contig数目的百分比作为指标，其中前3个指标都反映的是组装得到的序列的整体长度，这3个指标数值越大，表示组装结果中长序列越多，组装效果越好。对三个组装软件的组装结果如表1所示。
以前3个指标进行考查，对于PE125的数据，这3个指标都是IDBA-UD的最大，Trinity次之，SOAPdenovo最小，以CW2-1数据样本为例，IDBA-UD软件的3个指标的值分别是425 bp、422 bp、387 bp，Trinity对应的分别是321 bp、325 bp、277 bp，SOAPdenovo对应的分别是276 bp、179 bp、138 bp。对于PE150的数据来说，3个软件对应的3个指标的值相差不大，以AC142数据样本为例，IDBA-UD的3个指标分别是540 bp、543 bp、470 bp，Trinity对应的分别是444 bp、424 bp、341 bp，SOAPdenovo对应的分别是454 bp、459 bp、431 bp。对于最长contig这一指标，3个软件组装结果相差比较大，但是Trinity软件的该指标都是最优的，以CW2-1数据样本为例，使用Trinity软件进行组装最长contig达到8469 bp，而IDBA-UD和SOAPdenovo软件分别只有3661 bp和1484 bp，远远小于Trinity软件。对于大于1000的contig的数量、组装结果总contig数目以及大于1000 bp的contig数目占总contig数目的百分比这3个指标进行比较，Trinity软件优势非常明显，使用该软件组装得到的大于1000 bp的contig数量以及组装得到的总contig数量都是最多的，且远远高于另外2个软件，以CW2-1数据样本为例，Trinity软件组装结果中共有23,571条序列，其中大于1000 bp的有198条，对于IDBA-UD这2个值分别是1737和26，SOAPdenovo这2个值分别是2882和2。这反映了Trinity软件的组装结果中包含了丰富的序列信息且长序列最多，这为后续进行序列信息挖掘提供了基础。
Table 1. Comparison of the assembly quality for SOAPdenovo, IDBA-UD and Trinity
表1. SOAPdenovo、IDBA-UD和Trinity 3个软件组装结果比较
3.3. 组装得到病毒序列
使用Trinity软件对处理后的数据进行组装，比对病毒数据库后对结果进行整理分析。在收集的数据中共找到了10种植物病毒，9种RNA病毒和1种DNA病毒，总计2129条植物病毒序列，其中314条比对长度超过1000 bp，大部分序列比对结果的一致性都在95%以上，其中114条覆盖90%以上的基因组，接近数据库中的病毒全长序列。
4. 讨论
4.1. 不同软件的k-mer选择
使用SOAPdenovo对序列进行组装，k-mer的取值对基因组组装结果影响较大。本研究参考了其他生物信息领域科研人员的建议，k-mer设置范围取reads长度的55%~85%，且为奇数。因此，对于PE125的数据，k-mer选择69~107，对于PE150的数据，k-mer选择83~127。
使用IDBA-UD (v1.2.0)进行组装不需要设置不同大小的k-mer分步运行，该软件可以从小的k-mer开始到大的k-mer进行递增计算，多个k-mer一次就可运行完成。合并后的数据再进行组装，设置最大、最小的k-mer参数，与上一步选择相同的k-mer，但是由于IDBA-UD软件支持的最大k-mer是123，因此对于150 bp的reads最大k-mer只能选择到123。
对于Trinity软件，由于该软件不需要手动设置k-mer大小，该软件可以自动计算最优k-mer，因此本研究并未对Trinity软件进行上述研究。
4.2. 病毒基因组装
在本研究中利用SOAPdenovo、IDBA-UD、Trinity三种软件对测序数据进行组装，通过对组装效果的衡量对比，我们发现Trinity软件运行病毒序列组装，信息最为丰富，长序列最多，但运行时间长；而SOAPdenovo和IDBA-UD这两款软件耗时都比较短，但是组装效果不佳，序列信息少且长序列少，所以在追求组装质量的前提下，本研究建议选择Trinity进行全部RNA-seq数据的组装。但从在耗时以及消耗资源的角度看，SOAPdenovo和IDBA-UD有非常明显的优势，相同线程、内存资源条件下，组装同一个样本，Trinity所消耗的时间是SOAPdenovo和IDBA-UD的几十倍。综上所述，Trinity软件组装效果最好，Smith等人和Sparks等人分别对艺神袖蝶和茶翅蝽进行转录组测序，使用Trinity进行数据拼接组装，都发现和鉴定了新病毒 [16] [17]。因此推荐使用Trinity软件进行基于转录组数据的病毒基因组组装。本研究中，利用Trinity软件对未匹配到马铃薯参考基因组的数据进行组装，组装结果比对病毒数据库，共发现9种RNA病毒和1种DNA病毒。其中有常见的6种马铃薯病毒，所占比例最大且为90.94%。此外还发现了TuMV等几种目前尚未报道对马铃薯有侵染能力的病毒，表明了使用Trinity组装马铃薯病毒序列的可靠性。但是，如果考虑到时间成本的话推荐使用IDBA-UD软件。而SOAPdenovo组装“碎片”较多，不适合病毒基因组组装。
利用生物学分析软件构建基于RNA-seq数据的马铃薯病毒基因组挖掘平台，分析RNA-seq数据中病毒多样性，研究病毒变异、进化，发现和鉴定可能存在的新病毒或者尚未报道的对马铃薯有侵染能力的病毒，为马铃薯与病毒互作研究、马铃薯抗性材料筛选奠定了基础，对马铃薯抗病毒育种具有重要意义。
4.3. 测序数据深度挖掘
高通量测序技术由于其快速、准确、低成本的特点，已经被广泛应用于生物学研究中。近年来，越来越多的植物基因组测序完成，与马铃薯相关的转录组数据快速增长。同时测序技术也在不断创新，二代测序技术由于其通量高、精度高、信息丰富、成本低的特点已经广泛应用于生命科学研究的各个领域，产生了大量的转录组测序数据。但是，大量的测序数据在满足了特定的生物意义分析之后，就被闲置起来，十分可惜。但是，我们可以从其它角度进行深度挖掘，从这些闲置的数据中发现潜藏的有价值信息。
基于以上理论和技术基础，在本研究中收集了实验室数据库中的转录组数据，进行深度挖掘，挖掘其中的病毒序列，分析能够侵染马铃薯的病毒类型、病毒多样性。但由于测序最初的目的是针对宿主马铃薯进行研究，所以大部分数据是宿主序列，病毒序列丰度不足，因此导致组装到完整病毒序列的数量有限。本研究针对转录组数据构建了基于RNA-seq数据的病毒组装平台，该平台能够用于分析RNA-seq数据中病毒序列多样性，能发现一些新的、有研究价值的病毒序列，具有良好的应用前景。
致谢
作者感谢陈汝豪、李春燕等同学为本研究提供的相关测序数据，感谢罗鸣博士在分析方法上提出的意见和建议。感谢农业农村部马铃薯生物学与生物技术重点实验室为本研究提供了良好的科研平台。
基金项目
国家自然科学基金项目(31971989)；农业部华中作物有害生物综合治理重点实验室/农作物重大病虫草害防控湖北省重点实验室开放基金课题(2019ZTSJJ4)；湖北省农业科技创新中心项目(2019-620-003-001)；国家大学生创新创业基金(202010504041)。