BACKGROUND
increases in global trade and movement are placing significant pressure on post entry quarantine systems, with an increase in the frequency of incursions of pathogens causing the emergence of diseases and pests that are both difficult and costly to eradicate and control  <cit> .
the challenge of maximising the benefits of global trade whilst minimising the negative impacts of biosecurity threats is one faced by most nations  <cit> .
historically, the geographical isolation of australia and new zealand, coupled with stringent quarantine screening measures, has provided protection from the introduction of exotic pests and pathogens that have the potential to harm human health, agriculture, the environment and the economy.
plant biosecurity is defined as “a set of measures designed to protect crops from emergency plant pests at national, regional and individual farm level”  <cit> .
the diagnosis of viral pathogens is a crucial component of plant biosecurity surveillance, required to prevent the potential introduction of exotic plant viruses and viroids.
existing ‘specific’ serological and molecular detection methods such as enzyme-linked immunosorbent assay , polymerase chain reaction , or nucleic acid spot hybridization, while highly sensitive, are limited by their ability to detect only known plant viruses/viroids.
these methods lack the capacity to detect unknown, poorly characterised or highly variable viral pathogens  <cit> .
furthermore the host range of many viral pathogens is not defined and known exotic viruses/viroids could be missed if these infect new plant species for which standard screening assays are not applied.
if pathogens are not initially detected via these methods, more ‘investigational’ techniques may be applied, such as electron microscopy, host plant inoculation, or pcr using degenerate primers  <cit> .
the time and effort taken to screen imported plants using these existing methods has a direct economic impact, with plants that are currently imported into australia and new zealand spending up to two years in quarantine .
recent studies have demonstrated both the detection of viral pathogens and the identification of novel viruses by the deep sequencing of small rnas  of plant species .
rna silencing is a natural anti-viral defence system present in plants, insects and invertebrates that recognise dsrna viral genomes and/or viral intermediate sequences, and cleave them into small interfering rnas  of 21-24 nt in length  <cit> .
these virus-derived sirnas  accumulate in the small rna fraction of host plants making it amenable to identify viruses through a next generation sequencing  approach, even at extremely low viral titres and in symptomless infections  <cit> .
small rna ngs screening of viral pathogens is more cost- and time-effective compared with current detection methods.
the bottleneck for the uptake of ngs technology for routine surveillance and diagnosis of viral sequences is the lack of an automated bioinformatics pipeline that enables users to evaluate, scrutinize and modify all key steps of the analysis workflow including de novo assembly parameters with access to intermediate outputs  <cit> .
this is key to optimise and increase the confidence in the de novo assembled sequences with similarity to viral genomes.
two such automated pipelines have recently been made available, one designed to detect known and novel viruses through de novo assembly , and the other, limited to extract virnas through mapping onto a priori known viral reference genomes  <cit> .
virfind is available as a web-based graphical front-end interface, with users completing a sequence submission form and uploading sequence files via the virfind ftp server.
users are able to set their own parameters at several selected stages of the pipeline, such as choosing to map to a host reference genome, the expected value for blastn and blastx, or choosing to search conserved domains.
l miozzi and v pantaleo  <cit>  developed a pipeline to extract virnas through the open source galaxy web-based platform  <cit> .
this pipeline uses a reference guided approach to map reads to reference viral or viroid genomes, therefore virnas that correspond to viruses or viroids not present in public databases will be missed.
users upload a fastq file, and filter the results of the mapping to specific viruses or viroids of interest.
the resulting sam file can then be downloaded, and visualised using the java standalone tool misis  <cit> .
an automated, yet customisable bioinformatics pipeline for the detection and screening of viruses and viroids is required for the adoption of ngs technology by agencies without established bioinformatics expertise.
previously, we presented yabi, an analysis workflow environment that is able to create and reuse workflows, as well as manage large amounts of raw and processed data in a secure and flexible environment  <cit> .
yabi is accessed via a simple ‘drag and drop’ web-based environment by researchers without bioinformatics expertise or through the yabi command line for advanced users.
individual tools can be configured and easily incorporated into sophisticated workflows in real time.
importantly, comprehensive provenance for each workflow is kept, including input files and the parameters used for each tool, enabling researchers to track previous analyses performed and share optimised workflows with others.
files can be managed across different high performance computing storage infrastructures.
this study presents a novel automated internet-based bioinformatics toolkit for the detection of viruses and viroids utilising the online research environment yabi.
this toolkit offers users the flexibility to process small rna-seq samples using existing optimised workflows and/or to further customise available tools or add new tools into the web-based analytical environment.
it is envisaged that this resource will significantly reduce post entry quarantine ongoing costs and quarantine lead times.
furthermore, the web-based bioinformatics toolkit is customized to meet quarantine expectations, facilitate collaborations and inform policy makers.
methods
sample collection, rna extraction and ngs sequencing
imported plants and positive control samples were grown in quarantine glasshouse facilities until sample collection.
plants were grown under natural lighting with a daytime temperature of approximately 22 °c.
for each plant sample one or more leaves were collected prior to rna extraction.
total rna and/or small rna enriched fraction  were extracted from approximately 10 mg of tissue using the mirvana microrna isolation kit  following manufactures instructions.
collected samples were stored at -80 °c within quarantine facilities until shipped to the beijing genomics institute .
libraries were prepared using the truseq small rna sample prep kit  and sequenced with 50 bp single-end  reads deep sequencing of collected small rna samples  on an illumina hiseq <dig>  we sequenced  <dig> quarantined plant samples .
small rna-seq datasets has been submitted to the short read archive  under the bioproject prjna <dig>
selection of small rna assembler and scaffolding tools
we compared velvet  <cit> , spades  <cit> , abyss  <cit>  and soapdenovo  <cit>  assemblers using twelve selected small rna-seq samples collected from distinct plant species generated in this study .
we tested de novo assembly using individual kmer lengths of  <dig> ,  <dig> ,  <dig>  and  <dig>  as well as combined kmer sets of  <dig> , <dig>  and  <dig> , <dig>  .
assembled contigs were further scaffolded using cap <dig> using optimised parameters for short overlaps   <cit> .
additionally, merging and scaffolding of contigs produced by two or all three assemblers were also evaluated.
assembly statistics were calculated using the quality assessment tool for genome assemblies   <cit> .
overview of the automated viral diagnosis and surveillance toolkit
the viruses and viroids surveillance and diagnosis  bioinformatics toolkit was developed utilising yabi  <cit> , an open source internet-based analytical environment that allows for the customisation of tools and scripts into analysis workflows  <cit> .
yabi has five tabs, namely, ‘jobs’, ‘design’, ‘files’, ‘account’ and ‘admin’ tabs, where the later is only visible to a person or group responsible for the maintenance and further customisation of the yabi platform  <cit> .
the ‘jobs’ tab allows visualising and downloading results from prior jobs.
the ‘design’ tab enables re-use of existing optimised workflows, design of modified versions of existing workflows, and the construction of new analysis workflows.
the ‘files’ tab present files and directories of all available backend resources  to the user  <cit> .
the ‘account’ tab enables a user to easily modify their password information to their yabi account.
the ‘admin’ tab facilitates the management and addition of new computational tools into the yabi environment.
new features of the yabi platform include: i) save and share workflows; ii) fetch data from public repositories; iii) submission of processed data to specialised databases such as national or international patient registries; and iv) enables ‘bioinformatics on demand’ analyses through the deployment of cloud instance at the beginning of a computational workflow and its obliteration at the final step of the data processing and analysis workflow.
the vsd toolkit has three versions of the ‘virus and viroid detection’ workflows , with users able to choose from three subsets of small rna read lengths .
existing automated workflows can be reused or modified and saved .
additional workflows such as the ‘detecting novel viroids’ and ‘mapping reads onto a reference genome’ are also available, and can be run as a separate job, or added to the ‘virus and viroid detection’ workflows .fig.
<dig> workflow for the bioinformatics vsd toolkit for the discovery of viruses and viroids.
three versions of the workflow are available, with users able to choose from three options of read lengths .
★ indicates tools where users may change options or parameters as desired.
✚ indicates tools and workflows users can add or remove to the pipeline as by default these are not part of the automated vsd pipeline.
★ ✚ both ★ and ✚ options are available to users.
*results from this optional step are not used in downstream analyses
virus and viroid detection workflow
files of small rna reads in fastq format  are first uploaded through the ‘files’ tab in yabi.
files may be uploaded directly from a personal computer or transferred to a yabi directory from another high performance computing storage location.
once the file is uploaded, users then navigate to the ‘design’ tab, where they choose the saved workflow of interest .
users are also able to build their own workflows, by simply dragging and dropping tools into the workflow area.
the first stage of the workflow is the ‘select file’ tool.
the fastq file of interest is then selected.
if an adapter trimming step needs to be performed, users can add in the ‘fastx_clipper’ tool  to the workflow, and perform quality control checks using the ‘fastqc’ tool  <cit> .
reads then undergo quality trimming through the content dependent read trimming tool condetri , which trims and removes reads with low quality scores  <cit> .
minimum read length is set to 18 nt.
reads of the desired read length are then extracted through the ‘extract_reads_21–25 nt’ , ‘extract_reads_21–22 nt’ , or ‘extract_reads_ <dig> nt’ tools with de novo assemblies of contigs performed with spades  with kmer sizes set to  <dig>  and  <dig>  <cit> .
overlapping spades contigs are then merged with cap <dig>   <cit> .
contigs greater than or equal to 40 nt are then extracted using an in-house python script , and blast searched against databases generated from all plant, virus and viroid sequences populated by the entrez search query’s for viruses “txid <dig> ”, not cellular organisms “txid <dig> ”, viroids “txid <dig> ” and plants “txid <dig> ”.
for blastn, the task is set to blastn short, and for both blastn and blastx, the maximum number of aligned sequences  is set to  <dig> and the expected value  set to 1e− <dig>  blastn and blastx results are written out in a customised tabular format and extracted through an in-house script  into ‘plant’ or ‘virus and viroid’ blast output files.
contigs are also extracted into fasta output files through an in-house script  into subsets with a blast hit to plant or virus and viroid sequences, and contigs with no blast hits.
the blast output files are then parsed through an in-house script, which produces several output files in csv format including all blast results in tabular format , blast results with alignment lengths >39 nt, and a summary file which reports the genbank id of the virus or plant hit, the name of the plant or virus, the number of contig hits, the average percent sequence identity of the hit to the virus, the alignment length, the length of the virus or plant sequence, and the percentage coverage by contigs of a virus or plant sequence.
these statistics are calculated using the bioconductor’s genomicranges package .
detecting novel viroids workflow
the output file ‘no_hits_contigs.fasta’ from the ‘extract contigs’ tool represent sequences with no blastn and blastx sequence similarity to viral and plant sequences.
these sequences are further filtered to extract contigs with lengths between 200– <dig> b that are typical for viroids.
an in-house script is then utilised to evaluate the sequence similarity and overlap of both 5′-end and 3′-end of selected sequences.
sequences with overlapping ends with 100% sequence similarity are reported as putative circular viroid candidates.
users are recommended to further inspect identified candidate circular sequences.
for example evaluate sequence similarity to non-coding rna databases such as rfam  <cit>  and mirbase  <cit> , which are not part of the vsd toolkit.
viroids from the pospiviroidae  and some from the avsunviroidae  form hairpin-like rna secondary structures  <cit> .
thus, such confirmation can be evaluated in filtered candidate novel viroid circular sequences using rnafold  <cit> .
furthermore, the expression of candidate sequences passing all filtering steps can be evaluated in multiple tissues of the infected plant and/or its progeny to validate i) the de novo assembled circular sequence, and ii) provide independent evidence of its expression in multiple tissues and/or individuals.
mapping small rnas onto reference genome workflow
mapping of virnas onto identified viral genomes from similarity searches typically provide a broader coverage of the viral pathogen sequence as compared to de novo assembled contigs.
quality trimmed reads or the subset of reads  can be mapped against a reference genome of choice  through bowtie <dig>  <cit> .
optionally, the samtools suite   <cit>  can be used to sort and index aligned reads.
the resulting alignment file  can then be downloaded and viewed using the java standalone tool misis  <cit> .
unique features of the vsd toolkit
the major unique feature of our vsd toolkit as compared to virfind  <cit> , is the ability to exclusively use 21–22 nt small rna reads for the de novo assembly of viral sequences.
assembly of viral sequences with this set of reads directly reflect the plant endogenous antiviral response mediated by dicer <dig> and dicer <dig>  <cit> .
additionally, we provide an assembly pipeline that uses 24 nt small rnas overlapping the expected size for endogenous heterochromatin and transposon derived sirnas  <cit> .
this 24 nt pipeline identifies viral sequences potentially integrated in the host genome, particularly if they are not detected using 21–22 nt pipeline.
finally, we provide a 21–25 nt pipeline for users to compare their outputs against other published work that typically use a broad range of small rnas and/or compare with the results from the targeted 21–22 nt and 24 nt pipelines.
another unique feature of the vsd toolkit is the ability to modify the parameter options for most of the individual steps in the workflow .
additionally, the vsd toolkit uses optimised spades de novo assembly settings that yield improved results as compared to other tested assemblers .
the similarity screening of viruses in vsd toolkit is run in parallel using both blastn and blastx  <cit>  for all de novo assembled contigs against viruses and viroid sequences in the nt and nr databases , respectively.
the top five database hits for each de novo assembled contig are reported improving coverage of specific isolate/strain viral sequences and/or preventing false negative results when a top viral hit is annotated as “synthetic sequence”.
virfind runs blast screening in a staggered manner, reporting first nucleotide top hits against viral sequences, and then for contigs with negative blastn results, a blastx screening is conducted reporting the best hit  <cit> .
finally, the vsd toolkit also provides a list of potential viroid-like circular sequences with no sequence similarity to any nucleotide sequence in public databases.
deployment of the toolkit
during optimisation and testing of the bioinformatics toolkit, the workflows were run on a dynamic sge cluster located on amazon web services  elastic compute cloud , which allows compute nodes  to be easily added or removed as required.
yabi and the bioinformatics toolkit may be deployed on a variety of high performance computing resources.
the source code is available from https://github.com/muccg/yabi.
RESULTS
selection of de novo assembler for viral sequences
preliminary de novo assembly of  <dig> quarantine samples sequenced in this study  were evaluated with velvet  <cit> , spades  <cit> , abyss  <cit>  and soapdenovo  <cit>  assemblers using individual kmer lengths of  <dig> ,  <dig> ,  <dig>  and  <dig>  as well as kmer sets of  <dig>   <dig>   <dig>  and  <dig>   <dig>   <dig>   <dig> .
we found that spades k15-17- <dig> and spades k15-17-19- <dig> coupled with cap <dig>  <cit>  scaffolding produced the longest assembled sequences .
furthermore, we found that spades k15-17- <dig> yielded a larger number of assembled bases than spades 15-17-19- <dig> .
spades uses a unique approach to progressively build the assembled contigs using first the shortest kmer size in the first round, and then it builds upon the results of the previous round to continue to assemble sequences using additional user-defined longer kmer sizes .
this feature is not available to velvet, abyss and soapdenovo yielding a large fraction of redundant contigs assembled by individual kmers sizes that account for the increased total number of assembled bases prior and after cap <dig> scaffolding , respectively.
overall we found that spades k15-17- <dig> coupled with cap <dig> scaffolding produced the best results as compared to other tools and kmer settings tested .fig.
<dig> comparison of spades, velvet, abyss and soapdenovo assemblies using various kmer settings prior and after cap <dig> scaffolding.
twelve samples were selected for the comparisons .
a average longest assembled contigs for each kmer setting using 21–22 nt reads are shown.
see text for kmer settings used.
assemblies after cap <dig> scaffolding are denoted with ‘_cap3’.
b average total number of assembled bases for each kmer setting using 21–22 nt reads is shown.
c average longest assembled contigs for each kmer setting using 21–25 nt reads are shown.
d average total number of assembled bases for each kmer setting using 21–25 nt reads is shown
de novo assembly of viral genomes using small rna subsets
typically de novo assembly of viral sequences has been conducted using  <dig> to 25 nt long reads.
in this study, we compared the de novo assembly outcomes of viral sequences using three subsets of small rna reads: i) 21–25 nt, ii) 21–22 nt, and iii) 24 nt long reads.
we identified contigs with sequence similarity to known viral genomes in  <dig> samples .
in agreement with the known plant antiviral response activities of dicer <dig> and dicer <dig> enzymes  <cit> , most of the identified viral sequences were assembled using 21–22 nt reads .
only a citrus endogenous pararetrovirus was more effectively assembled using 24 nt long reads than with 21–22 nt reads .
endogenous pararetrovirus sequences  belonging to the plant virus family caulimoviridae have been discovered in the genomes of a wide range of angiosperms preferentially integrated into at dinucleotide repeats  <cit> .
transcription of eprv sequences along with flanking repeats may trigger a dicer3-mediated silencing pathway that cleaves repeats and heterochromatin sequences into 24 nt long small interference rnas  <cit> .
users are encouraged to run in parallel 21–22 nt and 24 nt assembly workflows to screen for distinct types of viruses.fig.
<dig> comparison of viral genome coverage using subsets of small rna reads.
a percentage coverage of viral genomes by contigs assembled using 21–25 nt, 21–22 nt or 24 nt small rnas.
s1 = vitis vinifera; s2 = prunus persica; s3 = rubus idaeus; s9 = miscanthus sinensis; s10 = citrus medica; s12 = citrus sp.
; s17 = pennisetum advena; s21 = actinidia deliciosa; s22 = nicotiana tabacum; grspav = grapevine rupestris stem pitting-associated virus; gvb = grapevine virus b; gfkv = grapevine fleck virus; glrav-3 = grapevine leafroll-associated virus; hsvd = hop stunt viroid; pnrsv = prunus necrotic ringspot virus; rbdv = raspberry bushy dwarf virus; rynv = rubus yellow net virus; spsmv-1 = sweetpotato symptomless mastrevirus 1; smov = strawberry mottle virus; srmv = sorghum moasaic virus; ctv = citrus tristeza virus; cvd-vi = citrus viroid vi; mcdv = maize chlorotic dwarf virus; acvb = actinidia virus b; ctv = tvcv = tobacco vein-clearing virus.
b mapping of high quality reads with no mismatches onto a reference grapevine fleck virus
gvba
pnrsv rna2b
pnrsv rna3b
rbdv rna <dig> c
miscanthus sinensis
cvd-vid
nicotiana tabacum
virus name: grspav grapevine rupestris stem pitting-associated virus, gvb grapevine virus b, gfkv grapevine fleck virus, glrav- <dig> grapevine leafroll-associated virus, hsvd hop stunt viroid, pnrsv prunus necrotic ringspot virus, rbdv raspberry bushy dwarf virus, rynv rubus yellow net virus, pcv  persimmon cryptic virus , acv  arisotelia chilensis virus  <dig> , spsmv- <dig> sweetpotato symptomless mastrevirus  <dig>  smov strawberry mottle virus, srmv sorghum mosaic virus, ctv citrus tristeza virus, cvd-vi citrus viroid vi, mcdv maize chlorotic dwarf virus, acvb actinidia virus b, tvcv tobacco vein-clearing virus, pvy potato virus y
apositive control for gvb identified using biological indexing at post entry quarantine
bpeq positive control for pnrsv detected using elisa
cpeq positive control for rbdv identified using elisa
dpeq positive control for ctv and cvd-vi detected using pcr-based assays
to assign sequence similarity of de novo assembled contigs to plant, viruses and viroid sequences using either blastn or blastx against nt and nr databases, respectively, a threshold of e-value 1e− <dig> is applied.
the mean average alignment length for blastn and blastx hits assembled using the 21–22 nt pipeline was  <dig>  bp and  <dig>  bp, respectively .
the minimal length of de novo assembled contigs that are compared against public databases is 40 bp.
we found that for blastn alignment lengths of less than 50 bp the average nucleotide sequence identity was  <dig> % ranging from  <dig> % to up to 100% .
in contrast, we identified two contigs with the lowest nucleotide sequence identities of  <dig> % and  <dig> % that had an alignment length of 944 bp and 178 bp, respectively .
similar results were observed for de novo assembled contigs produced using either the 21–25 nt  or 24 nt  pipelines.
our approach identifies viral sequences with high nucleotide sequence similarity to known viruses deposited in the nt database.
to identify more divergent viral sequences the blastx similarity to the nr database is used.
we identified that the minimal protein alignment length reported for an e-value of 1e− <dig> was  <dig> amino acids corresponding to 84 bp.
inspection of protein alignment shorter than  <dig> amino acids for de novo assembled contigs produced using the 21–22 nt pipeline showed that the average amino acid sequence identity was  <dig> % ranging from  <dig> % to 100% .
similar blastx results were observed for contigs assembled using the 21–25 nt pipeline .
although most samples processed by the 24 nt pipeline also showed similar blastx results, we identified a pennisetum advena  and a citrus latifolia  that did not yield a hit against viral proteins .
most viral sequences found in this study have sequence similarity to ssrna  viruses .
to evaluate the length distribution of small rnas making up the viral sequences we mapped quality trimmed reads onto a reference grapevine fleck virus  genome, and showed that out of  <dig>  mapped reads with zero mismatches  <dig> %,  <dig> % and  <dig> % had lengths of 21 nt, 22 nt and 24 nt, respectively .
we also aligned reads onto a citrus tristeza virus  genome, yielding  <dig>  million mapped reads with no mismatches, of which  <dig> %,  <dig> % and  <dig> % were 21 nt, 22 nt and 24 nt in length, respectively .
these findings correlate with the notion that dicer <dig> is the main antiviral response enzyme in leaf tissues  <cit> .
we detected three viroids including hop stunt viroid  and citrus viroid vi  in citrus, and hsvd in grapevine .
interestingly, all three viroid sequences were similarly assembled using either 21–22 nt or 24 nt small rna subsets.
mapping of high quality trimmed reads with no mismatches onto the reference hsvd genome  yielded  <dig>  mapped reads, of these  <dig> %,  <dig> % and  <dig> % had a length of 21 nt, 22 nt and 24 nt, respectively .
the identified viroids in this study are predicted to form a hairpin-like rna secondary structure that may be recognised not only by the antiviral response pathway but also by other plant small rna pathways .
the possible intervention of two or more plant small rna pathways may result in the cleavage of viroids into 21 nt, 22 nt and 24 nt long sequences.
furthermore, to identify new to science viroids the complete assembly of their genome using 21–22 nt or 24 nt small rna subsets can be applied as filtering criteria.
to evaluate the fraction of small rna reads making up de novo assembled sequences we mapped high quality adaptor clipped reads onto assembled contigs using bowtie  <cit>  with up to  <dig> mismatches.
we found on average that  <dig> % of the reads were not mapped onto assembled contigs  .
inspection of the mapped reads, showed that on average  <dig> %,  <dig> %,  <dig> % and  <dig> % of these reads were anchored onto assembled contigs with sequence similarity to plants, viruses, sequences with no similarity to plants or viral sequences  and viroids, respectively .
the fraction of reads anchored onto viral sequences ranged from  <dig> % to up to  <dig> %.
we then aimed to define the minimum amount of small rna data that should be collected for the surveillance and diagnosis of viral sequences.
as a case study we used a prunus persica  sample and tested the diagnosis of the tripartite prunus necrotic ringspot virus  using subsets of 1 m, 2 m, 5 m, and 10 m high quality small rna reads.
we found that when using 10 m reads for de novo assembly we detected  <dig> %,  <dig> % and  <dig> % of the pnrsv rna <dig>  rna <dig> and rna <dig>  respectively .
the use of smaller subsets of reads resulted in limited assembly of pnrsv rna <dig> and rna <dig> sequences.
to increase the confidence in the detected viral sequences we propose to use at least 10 m small rna reads per quarantined sample.
it is critical to detect exotic viral sequences in imported quarantined plant samples to protect domestic plant industries, native plant biodiversity and prevent economic and social impacts to the broader community.
assembly issues using combined 21–25 nt small rnas
collected small rnas from plant host samples contain the products of several biological pathways including antiviral response, heterochromatin and transposon silencing, and microrna biogenesis.
we found that conducting de novo assembly using 21–25 nt small rna sequences resulted in a partial assembly of a potato virus y  genome .
two contigs were assembled, one encoding most of the polyprotein region and another encoding part of the 3′end region of the genome.
the shorter contig included 184 bp at the 5′-end with no sequence similarity to the terminal region of the pvy polyprotein, instead reverse complementation of these 184 bp produced a match to the beginning of the pvy genome.
in contrast to these findings, de novo assembly using only a 21–22 nt subset of small rnas resulted in the complete assembly of the pvy sequence including  <dig> and  <dig> extra bases in their 5′-end and 3′-end terminal regions, respectively .
these findings suggest that improved and accurate assembly results are obtained using specific 21–22 nt small rnas produced by the plant antiviral response pathway  <cit> .fig.
<dig> assembly of a potato virus y genome using 21–25 nt and 21–22 nt reads
comparison of toolkit outcomes with virfind
we subjected  <dig> raw rna-seq datasets to evaluate the performance of vsd toolkit as compared to virfind  <cit> .
summary of results with the virfind pipeline are outlined in additional file  <dig>  virfind identified viral sequences in all  <dig> samples, of these in  <dig> samples the vsd toolkit found the same viral sequences , while  <dig> samples showed viral hits unique to the virfind pipeline .
inspection of the later samples showed hits to viral sequences from plants, algae, fungi, invertebrate and vertebrates hosts .
we aligned the adaptor clipped reads onto the assembled virfind contigs for these samples and found reads from 19 nt to up to 29 nt aligned onto contigs for these samples highlighting the broad spectrum of reads that are used in the virfind pipeline .
to test if these unique hits in the virfind pipeline may potentially be chimeric sequences with host plant sequences, we conducted a blastx screening against a combined plant and virus protein database and reported the top hit.
additional file  <dig> shows that  <dig> % of the virfind contigs annotated as viral sequences have a top hit to a plant or an insect virus.
these findings highlight further the susceptibility to assemble chimeric sequences and/or plant sequences with similarity to viruses  when using a broad spectrum of small rna read lengths collected from plant hosts.
comparison of the viral sequences commonly assembled by virfind and our approach showed that the vsd toolkit assembled an average of  <dig> % and  <dig> % more of the  <dig> viral sequences found in  <dig> plant samples using the 21–25 nt and 21–22 nt workflows, respectively .
the largest increases in the assembly of the known viral sequences were for the citrus endogenous pararetrovirus , grapevine freck virus , potato virus y  and prunus necrotic ringspot virus rna <dig> .
only in the case of the grapevine leafroll-associated virus  <dig>  did virfind assemble a greater  proportion of the sequences than the vsd toolkit  .
both the vsd toolkit and virfind detected similarly all three viroids sequences .
in the case of the citrus medica endogenous pararetrovirus  the vsd toolkit assembled  <dig> %,  <dig> % and  <dig> % of the genome using 21–25 nt, 21–22 nt and 24 nt assembly workflows, respectively.
in contrast, virfind only assembled  <dig> % of the cmeeprv genome sequence .
overall we show that vsd toolkit produced improved and more accurate results than virfind.
it is critical to utilise the specific by-products of the plant immune defence pathway against viruses and viroids to assemble highly accurate viral sequences that reflect the active host antiviral response.
detection of exotic viral pathogens at the border is critical to safeguard plant industries and their access to national and international trade markets.
CONCLUSIONS
we have implemented an automated viral surveillance and diagnosis toolkit using the yabi web-based analytical environment that enables improved detection of viruses and viroids pathogens.
we found that all single stranded rna  viruses found in this study were assembled using exclusively 21–22 nt small rnas, while viroids were equally assembled using 21–22 nt or 24 nt subsets.
the use of specific small rna subsets increases specificity of the identified viral sequences as it reflects the antiviral response activity of the plant hosts.
additional files
additional file 1: adaptor clipped and quality trimmed read statistics of small rna-seq datasets.
clean reads = adaptor clipped small rna reads provided by bgi service provider; qc = quality control consisting of adaptor clipping and poor base call trimming; n.a.
= not available.
additional file 2: yabi viral surveillance and diagnosis  toolkit user-guide.
additional file 3: boxplot comparison of spades assemblies using individual kmers, sets of kmers and cap <dig> scaffolding.
additional file 4: paired sample t-test of contigs produced by spades, velvet, abyss and soapdenovo assemblers.
our optimised ‘spades k15-17-19_cap3’   assembly was compared against all other kmer settings and/or tools.
a) spades - longest assembled contigs ; b) spades - total number of assembled bases for contigs > =100 nt ; c) spades - longest assembled contigs ; d) spades - total number of assembled bases for contigs > =100 nt ; e) comparison of longest assembled contig using velvet against sp22_k15-17-19_cap3; f) comparison of longest assembled contig using abyss against sp22_k15-17-19_cap3; g) comparison of longest assembled contig using soapdenovo against sp22_k15-17-19_cap <dig>
additional file 5: blastn and blastx sequence similarity statistics for de novo assembled contigs using the 21–22 nt , 21–25 nt  and 24 nt  pipelines, respectively.
twelve small rna-seq samples generated in this study with nucleotide similarity to viral pathogens were compared.
see additional file  <dig> for details of the selected samples.
additional file 6: statistics of the viral sequences assembled using 21–25 nt, 21–22 nt and 24 nt small rna subsets.
*virus name: grspav = grapevine rupestris stem pitting-associated virus; gvb = grapevine virus b; gfkv = grapevine fleck virus; glrav-3 = grapevine leafroll-associated virus; hsvd = hop stunt viroid; pnrsv = prunus necrotic ringspot virus; rbdv = raspberry bushy dwarf virus; rynv = rubus yellow net virus; spsmv-1 = sweetpotato symptomless mastrevirus 1; smov = strawberry mottle virus; msimv = miscanthus sinensis mosaic virus; ctv = citrus tristeza virus; cvd-vi = citrus viroid vi; mcdv = maize chlorotic dwarf virus; acvb = actinidia virus b; tvcv = tobacco vein-clearing virus.
**virus type: rt = retro-transcribing virus; ssrna + = ssrna positive-strand virus; ssdna = single-stranded dna virus; dsrna = double-strand rna virus; v = viroid.
additional file 7: mapping of high quality adaptor-clipped and quality trimmed small rnas with no mismatches onto the reference genomes: a) citrus tristesa virus .
domains of the ctv genomes are denoted.
b) hop stunt viroid .
additional file 8: predicted rna secondary structure of viroids found in this study.
minimal free energy rna secondary structure encoding base-pair probabilities are shown for s1_hsvd , s10_cvd-vi , and s10_hsvd .
additional file 9: distribution of mapped and unmapped small rna reads for quarantined samples generated in this study.
a) percentage of unmapped and mapped reads onto de novo assembled contigs with sequence similarity to plants; viruses, viroids and unknown are shown.
b) distribution of the percentage of mapped reads for each sample along with the total number of distinct viral sequences  are shown.
details for each sample can be found in additional file  <dig>
additional file 10: de novo assembly of pnrsv rna <dig>  rna <dig> and rna <dig> viral sequences using five randomly generated subsets of 1 m, 2 m, 5 m, 10 m and all adaptor clipped small rna reads.
additional file 11: summary of virfind screening of small rna-seq datasets collected from quarantined plants.
a) virfind blastn results.
identified viruses/viroids were also identified using our vsd toolkit.
b) virfind blastx results.
identified viruses/viroids were only found using the virfind pipeline.
hit to viruses from vertebrates, invertebrates, fungi or algae are shown in red background.
additional file 12: virfind blastx results for contigs with sequence simalrity to viruses or viroids.
additional file 13: mapping of adaptor clipped reads onto virfind assembled contigs for samples s <dig>  s <dig>  s <dig>  s <dig>  s <dig> and s <dig>  details of samples can be found in additional file  <dig>
additional file 14: virfind contigs annotated with hits to viruses were run through the yabi vsd toolkit pipeline.
top blastx hits  to a combined plant and virus protein database are reported.
virfind assembled contigs with a top hit to a plant virus are highlighted in green background.
additional file 15: comparison of genome coverage of known viral sequences by contigs assembled by the vsd toolkit and virfind analysis pipelines.
abbreviations
abyssa parallel assembler for short read sequence data
acv
arisotelia chilensis virus  <dig>
acvbactinidia virus b
awsamazon web services
blastbasic local alignment search tool
condetria content dependent read trimmer for illumina data
ctv
citrus tristeza virus
cvd-vicitrus viroid vi
ec2elastic compute cloud
elisaenzyme-linked immunosorbent assay
eprvendogenous pararetrovirus sequences
gfkvgrapevine fleck virus
glrav-3grapevine leafroll-associated virus
grspavgrapevine rupestris stem pitting-associated virus
gvbgrapevine virus b
hsvdhop stunt viroid
k15de novo assembly using kmer size of 15
k15-17-19de novo assembly using a combined kmer set of  <dig>   <dig> and 19
k15-17-19-21de novo assembly using a combined kmer set of  <dig>   <dig>   <dig> and 21
k17de novo assembly using kmer size of 17
k19de novo assembly using kmer size of 19
k21de novo assembly using kmer size of 21
mcdvmaize chlorotic dwarf virus
ngsnext generation sequencing
pcrpolymerase chain reaction
pcv persimmon cryptic virus
peqpost entry quarantine
pnrsvprunus necrotic ringspot virus
pvypotato virus y
rbdvraspberry bushy dwarf virus
rna-seqrna sequencing
rynvrubus yellow net virus
samsequence alignment/map format
sirnasmall interference small rna
smovstrawberry mottle virus
soapdenovoshort oligonucleotide analysis package for short read de novo assembly
spadesa genome assembly algorithm
spsmv-1sweetpotato symptomless mastrevirus 1
srmvsorghum mosaic virus
ssrnasingle-strand rna
tvcvtobacco vein-clearing virus
velvetalgorithm package for de novo genome assembly and short read alignment
virfindvirus find pipeline
virnavirus-derived sirna
vsdviral surveillance and diagnosis toolkit
yabiyet another bioinformatics interface

 increases in global trade and movement are placing significant pressure on post entry quarantine systems, with an increase in the frequency of incursions of pathogens causing the emergence of diseases and pests that are both difficult and costly to eradicate and control. the challenge of maximising the benefits of global trade whilst minimising the negative impacts of biosecurity threats is one faced by most nations. historically, the geographical isolation of australia and new zealand, coupled with stringent quarantine screening measures, has provided protection from the introduction of exotic pests and pathogens that have the potential to harm human health, agriculture, the environment and the economy. plant biosecurity is defined as “a set of measures designed to protect crops from emergency plant pests at national, regional and individual farm level”. the diagnosis of viral pathogens is a crucial component of plant biosecurity surveillance, required to prevent the potential introduction of exotic plant viruses and viroids. existing ‘specific’ serological and molecular detection methods such as enzyme-linked immunosorbent assay, polymerase chain reaction, or nucleic acid spot hybridization, while highly sensitive, are limited by their ability to detect only known plant viruses/viroids. these methods lack the capacity to detect unknown, poorly characterised or highly variable viral pathogens. furthermore the host range of many viral pathogens is not defined and known exotic viruses/viroids could be missed if these infect new plant species for which standard screening assays are not applied. if pathogens are not initially detected via these methods, more ‘investigational’ techniques may be applied, such as electron microscopy, host plant inoculation, or pcr using degenerate primers. the time and effort taken to screen imported plants using these existing methods has a direct economic impact, with plants that are currently imported into australia and new zealand spending up to two years in quarantine. recent studies have demonstrated both the detection of viral pathogens and the identification of novel viruses by the deep sequencing of small rnas of plant species. rna silencing is a natural anti-viral defence system present in plants, insects and invertebrates that recognise dsrna viral genomes and/or viral intermediate sequences, and cleave them into small interfering rnas of 21-24 nt in length. these virus-derived sirnas accumulate in the small rna fraction of host plants making it amenable to identify viruses through a next generation sequencing approach, even at extremely low viral titres and in symptomless infections. small rna ngs screening of viral pathogens is more cost- and time-effective compared with current detection methods. the bottleneck for the uptake of ngs technology for routine surveillance and diagnosis of viral sequences is the lack of an automated bioinformatics pipeline that enables users to evaluate, scrutinize and modify all key steps of the analysis workflow including de novo assembly parameters with access to intermediate outputs. this is key to optimise and increase the confidence in the de novo assembled sequences with similarity to viral genomes. two such automated pipelines have recently been made available, one designed to detect known and novel viruses through de novo assembly, and the other, limited to extract virnas through mapping onto a priori known viral reference genomes. virfind is available as a web-based graphical front-end interface, with users completing a sequence submission form and uploading sequence files via the virfind ftp server. users are able to set their own parameters at several selected stages of the pipeline, such as choosing to map to a host reference genome, the expected value for blastn and blastx, or choosing to search conserved domains. l miozzi and v pantaleo developed a pipeline to extract virnas through the open source galaxy web-based platform. this pipeline uses a reference guided approach to map reads to reference viral or viroid genomes, therefore virnas that correspond to viruses or viroids not present in public databases will be missed. users upload a fastq file, and filter the results of the mapping to specific viruses or viroids of interest. the resulting sam file can then be downloaded, and visualised using the java standalone tool misis. an automated, yet customisable bioinformatics pipeline for the detection and screening of viruses and viroids is required for the adoption of ngs technology by agencies without established bioinformatics expertise. previously, we presented yabi, an analysis workflow environment that is able to create and reuse workflows, as well as manage large amounts of raw and processed data in a secure and flexible environment. yabi is accessed via a simple ‘drag and drop’ web-based environment by researchers without bioinformatics expertise or through the yabi command line for advanced users. individual tools can be configured and easily incorporated into sophisticated workflows in real time. importantly, comprehensive provenance for each workflow is kept, including input files and the parameters used for each tool, enabling researchers to track previous analyses performed and share optimised workflows with others. files can be managed across different high performance computing storage infrastructures. this study presents a novel automated internet-based bioinformatics toolkit for the detection of viruses and viroids utilising the online research environment yabi. this toolkit offers users the flexibility to process small rna-seq samples using existing optimised workflows and/or to further customise available tools or add new tools into the web-based analytical environment. it is envisaged that this resource will significantly reduce post entry quarantine ongoing costs and quarantine lead times. furthermore, the web-based bioinformatics toolkit is customized to meet quarantine expectations, facilitate collaborations and inform policy makers. methods sample collection, rna extraction and ngs sequencing imported plants and positive control samples were grown in quarantine glasshouse facilities until sample collection. plants were grown under natural lighting with a daytime temperature of approximately 22 °c. for each plant sample one or more leaves were collected prior to rna extraction. total rna and/or small rna enriched fraction were extracted from approximately 10 mg of tissue using the mirvana microrna isolation kit following manufactures instructions. collected samples were stored at -80 °c within quarantine facilities until shipped to the beijing genomics institute. libraries were prepared using the truseq small rna sample prep kit and sequenced with 50 bp single-end reads deep sequencing of collected small rna samples on an illumina hiseq we sequenced quarantined plant samples. small rna-seq datasets has been submitted to the short read archive under the bioproject prjna selection of small rna assembler and scaffolding tools we compared velvet, spades, abyss and soapdenovo assemblers using twelve selected small rna-seq samples collected from distinct plant species generated in this study. we tested de novo assembly using individual kmer lengths of,, and as well as combined kmer sets of, and,. assembled contigs were further scaffolded using cap using optimised parameters for short overlaps. additionally, merging and scaffolding of contigs produced by two or all three assemblers were also evaluated. assembly statistics were calculated using the quality assessment tool for genome assemblies. overview of the automated viral diagnosis and surveillance toolkit the viruses and viroids surveillance and diagnosis bioinformatics toolkit was developed utilising yabi, an open source internet-based analytical environment that allows for the customisation of tools and scripts into analysis workflows. yabi has five tabs, namely, ‘jobs’, ‘design’, ‘files’, ‘account’ and ‘admin’ tabs, where the later is only visible to a person or group responsible for the maintenance and further customisation of the yabi platform. the ‘jobs’ tab allows visualising and downloading results from prior jobs. the ‘design’ tab enables re-use of existing optimised workflows, design of modified versions of existing workflows, and the construction of new analysis workflows. the ‘files’ tab present files and directories of all available backend resources to the user. the ‘account’ tab enables a user to easily modify their password information to their yabi account. the ‘admin’ tab facilitates the management and addition of new computational tools into the yabi environment. new features of the yabi platform include: i) save and share workflows; ii) fetch data from public repositories; iii) submission of processed data to specialised databases such as national or international patient registries; and iv) enables ‘bioinformatics on demand’ analyses through the deployment of cloud instance at the beginning of a computational workflow and its obliteration at the final step of the data processing and analysis workflow. the vsd toolkit has three versions of the ‘virus and viroid detection’ workflows, with users able to choose from three subsets of small rna read lengths. existing automated workflows can be reused or modified and saved. additional workflows such as the ‘detecting novel viroids’ and ‘mapping reads onto a reference genome’ are also available, and can be run as a separate job, or added to the ‘virus and viroid detection’ workflows.fig. workflow for the bioinformatics vsd toolkit for the discovery of viruses and viroids. three versions of the workflow are available, with users able to choose from three options of read lengths. ★ indicates tools where users may change options or parameters as desired. ✚ indicates tools and workflows users can add or remove to the pipeline as by default these are not part of the automated vsd pipeline. ★ ✚ both ★ and ✚ options are available to users. *results from this optional step are not used in downstream analyses virus and viroid detection workflow files of small rna reads in fastq format are first uploaded through the ‘files’ tab in yabi. files may be uploaded directly from a personal computer or transferred to a yabi directory from another high performance computing storage location. once the file is uploaded, users then navigate to the ‘design’ tab, where they choose the saved workflow of interest. users are also able to build their own workflows, by simply dragging and dropping tools into the workflow area. the first stage of the workflow is the ‘select file’ tool. the fastq file of interest is then selected. if an adapter trimming step needs to be performed, users can add in the ‘fastx_clipper’ tool to the workflow, and perform quality control checks using the ‘fastqc’ tool. reads then undergo quality trimming through the content dependent read trimming tool condetri, which trims and removes reads with low quality scores. minimum read length is set to 18 nt. reads of the desired read length are then extracted through the ‘extract_reads_21–25 nt’, ‘extract_reads_21–22 nt’, or ‘extract_reads_ nt’ tools with de novo assemblies of contigs performed with spades with kmer sizes set to and. overlapping spades contigs are then merged with cap. contigs greater than or equal to 40 nt are then extracted using an in-house python script, and blast searched against databases generated from all plant, virus and viroid sequences populated by the entrez search query’s for viruses “txid ”, not cellular organisms “txid ”, viroids “txid ” and plants “txid ”. for blastn, the task is set to blastn short, and for both blastn and blastx, the maximum number of aligned sequences is set to and the expected value set to 1e− blastn and blastx results are written out in a customised tabular format and extracted through an in-house script into ‘plant’ or ‘virus and viroid’ blast output files. contigs are also extracted into fasta output files through an in-house script into subsets with a blast hit to plant or virus and viroid sequences, and contigs with no blast hits. the blast output files are then parsed through an in-house script, which produces several output files in csv format including all blast results in tabular format, blast results with alignment lengths >39 nt, and a summary file which reports the genbank id of the virus or plant hit, the name of the plant or virus, the number of contig hits, the average percent sequence identity of the hit to the virus, the alignment length, the length of the virus or plant sequence, and the percentage coverage by contigs of a virus or plant sequence. these statistics are calculated using the bioconductor’s genomicranges package. detecting novel viroids workflow the output file ‘no_hits_contigs.fasta’ from the ‘extract contigs’ tool represent sequences with no blastn and blastx sequence similarity to viral and plant sequences. these sequences are further filtered to extract contigs with lengths between 200– b that are typical for viroids. an in-house script is then utilised to evaluate the sequence similarity and overlap of both 5′-end and 3′-end of selected sequences. sequences with overlapping ends with 100% sequence similarity are reported as putative circular viroid candidates. users are recommended to further inspect identified candidate circular sequences. for example evaluate sequence similarity to non-coding rna databases such as rfam and mirbase, which are not part of the vsd toolkit. viroids from the pospiviroidae and some from the avsunviroidae form hairpin-like rna secondary structures. thus, such confirmation can be evaluated in filtered candidate novel viroid circular sequences using rnafold. furthermore, the expression of candidate sequences passing all filtering steps can be evaluated in multiple tissues of the infected plant and/or its progeny to validate i) the de novo assembled circular sequence, and ii) provide independent evidence of its expression in multiple tissues and/or individuals. mapping small rnas onto reference genome workflow mapping of virnas onto identified viral genomes from similarity searches typically provide a broader coverage of the viral pathogen sequence as compared to de novo assembled contigs. quality trimmed reads or the subset of reads can be mapped against a reference genome of choice through bowtie. optionally, the samtools suite can be used to sort and index aligned reads. the resulting alignment file can then be downloaded and viewed using the java standalone tool misis. unique features of the vsd toolkit the major unique feature of our vsd toolkit as compared to virfind, is the ability to exclusively use 21–22 nt small rna reads for the de novo assembly of viral sequences. assembly of viral sequences with this set of reads directly reflect the plant endogenous antiviral response mediated by dicer and dicer. additionally, we provide an assembly pipeline that uses 24 nt small rnas overlapping the expected size for endogenous heterochromatin and transposon derived sirnas. this 24 nt pipeline identifies viral sequences potentially integrated in the host genome, particularly if they are not detected using 21–22 nt pipeline. finally, we provide a 21–25 nt pipeline for users to compare their outputs against other published work that typically use a broad range of small rnas and/or compare with the results from the targeted 21–22 nt and 24 nt pipelines. another unique feature of the vsd toolkit is the ability to modify the parameter options for most of the individual steps in the workflow. additionally, the vsd toolkit uses optimised spades de novo assembly settings that yield improved results as compared to other tested assemblers. the similarity screening of viruses in vsd toolkit is run in parallel using both blastn and blastx for all de novo assembled contigs against viruses and viroid sequences in the nt and nr databases, respectively. the top five database hits for each de novo assembled contig are reported improving coverage of specific isolate/strain viral sequences and/or preventing false negative results when a top viral hit is annotated as “synthetic sequence”. virfind runs blast screening in a staggered manner, reporting first nucleotide top hits against viral sequences, and then for contigs with negative blastn results, a blastx screening is conducted reporting the best hit. finally, the vsd toolkit also provides a list of potential viroid-like circular sequences with no sequence similarity to any nucleotide sequence in public databases. deployment of the toolkit during optimisation and testing of the bioinformatics toolkit, the workflows were run on a dynamic sge cluster located on amazon web services elastic compute cloud, which allows compute nodes to be easily added or removed as required. yabi and the bioinformatics toolkit may be deployed on a variety of high performance computing resources. the source code is available from https://github.com/muccg/yabi.  selection of de novo assembler for viral sequences preliminary de novo assembly of quarantine samples sequenced in this study were evaluated with velvet, spades, abyss and soapdenovo assemblers using individual kmer lengths of,, and as well as kmer sets of and. we found that spades k15-17- and spades k15-17-19- coupled with cap scaffolding produced the longest assembled sequences. furthermore, we found that spades k15-17- yielded a larger number of assembled bases than spades 15-17-19-. spades uses a unique approach to progressively build the assembled contigs using first the shortest kmer size in the first round, and then it builds upon the results of the previous round to continue to assemble sequences using additional user-defined longer kmer sizes. this feature is not available to velvet, abyss and soapdenovo yielding a large fraction of redundant contigs assembled by individual kmers sizes that account for the increased total number of assembled bases prior and after cap scaffolding, respectively. overall we found that spades k15-17- coupled with cap scaffolding produced the best results as compared to other tools and kmer settings tested.fig. comparison of spades, velvet, abyss and soapdenovo assemblies using various kmer settings prior and after cap scaffolding. twelve samples were selected for the comparisons. a average longest assembled contigs for each kmer setting using 21–22 nt reads are shown. see text for kmer settings used. assemblies after cap scaffolding are denoted with ‘_cap3’. b average total number of assembled bases for each kmer setting using 21–22 nt reads is shown. c average longest assembled contigs for each kmer setting using 21–25 nt reads are shown. d average total number of assembled bases for each kmer setting using 21–25 nt reads is shown de novo assembly of viral genomes using small rna subsets typically de novo assembly of viral sequences has been conducted using to 25 nt long reads. in this study, we compared the de novo assembly outcomes of viral sequences using three subsets of small rna reads: i) 21–25 nt, ii) 21–22 nt, and iii) 24 nt long reads. we identified contigs with sequence similarity to known viral genomes in samples. in agreement with the known plant antiviral response activities of dicer and dicer enzymes, most of the identified viral sequences were assembled using 21–22 nt reads. only a citrus endogenous pararetrovirus was more effectively assembled using 24 nt long reads than with 21–22 nt reads. endogenous pararetrovirus sequences belonging to the plant virus family caulimoviridae have been discovered in the genomes of a wide range of angiosperms preferentially integrated into at dinucleotide repeats. transcription of eprv sequences along with flanking repeats may trigger a dicer3-mediated silencing pathway that cleaves repeats and heterochromatin sequences into 24 nt long small interference rnas. users are encouraged to run in parallel 21–22 nt and 24 nt assembly workflows to screen for distinct types of viruses.fig. comparison of viral genome coverage using subsets of small rna reads. a percentage coverage of viral genomes by contigs assembled using 21–25 nt, 21–22 nt or 24 nt small rnas. s1 = vitis vinifera; s2 = prunus persica; s3 = rubus idaeus; s9 = miscanthus sinensis; s10 = citrus medica; s12 = citrus sp.; s17 = pennisetum advena; s21 = actinidia deliciosa; s22 = nicotiana tabacum; grspav = grapevine rupestris stem pitting-associated virus; gvb = grapevine virus b; gfkv = grapevine fleck virus; glrav-3 = grapevine leafroll-associated virus; hsvd = hop stunt viroid; pnrsv = prunus necrotic ringspot virus; rbdv = raspberry bushy dwarf virus; rynv = rubus yellow net virus; spsmv-1 = sweetpotato symptomless mastrevirus 1; smov = strawberry mottle virus; srmv = sorghum moasaic virus; ctv = citrus tristeza virus; cvd-vi = citrus viroid vi; mcdv = maize chlorotic dwarf virus; acvb = actinidia virus b; ctv = tvcv = tobacco vein-clearing virus. b mapping of high quality reads with no mismatches onto a reference grapevine fleck virus gvba pnrsv rna2b pnrsv rna3b rbdv rna c miscanthus sinensis cvd-vid nicotiana tabacum virus name: grspav grapevine rupestris stem pitting-associated virus, gvb grapevine virus b, gfkv grapevine fleck virus, glrav- grapevine leafroll-associated virus, hsvd hop stunt viroid, pnrsv prunus necrotic ringspot virus, rbdv raspberry bushy dwarf virus, rynv rubus yellow net virus, pcv persimmon cryptic virus, acv arisotelia chilensis virus, spsmv- sweetpotato symptomless mastrevirus smov strawberry mottle virus, srmv sorghum mosaic virus, ctv citrus tristeza virus, cvd-vi citrus viroid vi, mcdv maize chlorotic dwarf virus, acvb actinidia virus b, tvcv tobacco vein-clearing virus, pvy potato virus y apositive control for gvb identified using biological indexing at post entry quarantine bpeq positive control for pnrsv detected using elisa cpeq positive control for rbdv identified using elisa dpeq positive control for ctv and cvd-vi detected using pcr-based assays to assign sequence similarity of de novo assembled contigs to plant, viruses and viroid sequences using either blastn or blastx against nt and nr databases, respectively, a threshold of e-value 1e− is applied. the mean average alignment length for blastn and blastx hits assembled using the 21–22 nt pipeline was  bp and  bp, respectively. the minimal length of de novo assembled contigs that are compared against public databases is 40 bp. we found that for blastn alignment lengths of less than 50 bp the average nucleotide sequence identity was % ranging from % to up to 100%. in contrast, we identified two contigs with the lowest nucleotide sequence identities of % and % that had an alignment length of 944 bp and 178 bp, respectively. similar results were observed for de novo assembled contigs produced using either the 21–25 nt or 24 nt pipelines. our approach identifies viral sequences with high nucleotide sequence similarity to known viruses deposited in the nt database. to identify more divergent viral sequences the blastx similarity to the nr database is used. we identified that the minimal protein alignment length reported for an e-value of 1e− was amino acids corresponding to 84 bp. inspection of protein alignment shorter than amino acids for de novo assembled contigs produced using the 21–22 nt pipeline showed that the average amino acid sequence identity was % ranging from % to 100%. similar blastx results were observed for contigs assembled using the 21–25 nt pipeline. although most samples processed by the 24 nt pipeline also showed similar blastx results, we identified a pennisetum advena and a citrus latifolia that did not yield a hit against viral proteins. most viral sequences found in this study have sequence similarity to ssrna viruses. to evaluate the length distribution of small rnas making up the viral sequences we mapped quality trimmed reads onto a reference grapevine fleck virus genome, and showed that out of mapped reads with zero mismatches %, % and % had lengths of 21 nt, 22 nt and 24 nt, respectively. we also aligned reads onto a citrus tristeza virus genome, yielding million mapped reads with no mismatches, of which %, % and % were 21 nt, 22 nt and 24 nt in length, respectively. these findings correlate with the notion that dicer is the main antiviral response enzyme in leaf tissues. we detected three viroids including hop stunt viroid and citrus viroid vi in citrus, and hsvd in grapevine. interestingly, all three viroid sequences were similarly assembled using either 21–22 nt or 24 nt small rna subsets. mapping of high quality trimmed reads with no mismatches onto the reference hsvd genome yielded mapped reads, of these %, % and % had a length of 21 nt, 22 nt and 24 nt, respectively. the identified viroids in this study are predicted to form a hairpin-like rna secondary structure that may be recognised not only by the antiviral response pathway but also by other plant small rna pathways. the possible intervention of two or more plant small rna pathways may result in the cleavage of viroids into 21 nt, 22 nt and 24 nt long sequences. furthermore, to identify new to science viroids the complete assembly of their genome using 21–22 nt or 24 nt small rna subsets can be applied as filtering criteria. to evaluate the fraction of small rna reads making up de novo assembled sequences we mapped high quality adaptor clipped reads onto assembled contigs using bowtie with up to mismatches. we found on average that % of the reads were not mapped onto assembled contigs. inspection of the mapped reads, showed that on average %, %, % and % of these reads were anchored onto assembled contigs with sequence similarity to plants, viruses, sequences with no similarity to plants or viral sequences and viroids, respectively. the fraction of reads anchored onto viral sequences ranged from % to up to %. we then aimed to define the minimum amount of small rna data that should be collected for the surveillance and diagnosis of viral sequences. as a case study we used a prunus persica sample and tested the diagnosis of the tripartite prunus necrotic ringspot virus using subsets of 1 m, 2 m, 5 m, and 10 m high quality small rna reads. we found that when using 10 m reads for de novo assembly we detected %, % and % of the pnrsv rna rna and rna respectively. the use of smaller subsets of reads resulted in limited assembly of pnrsv rna and rna sequences. to increase the confidence in the detected viral sequences we propose to use at least 10 m small rna reads per quarantined sample. it is critical to detect exotic viral sequences in imported quarantined plant samples to protect domestic plant industries, native plant biodiversity and prevent economic and social impacts to the broader community. assembly issues using combined 21–25 nt small rnas collected small rnas from plant host samples contain the products of several biological pathways including antiviral response, heterochromatin and transposon silencing, and microrna biogenesis. we found that conducting de novo assembly using 21–25 nt small rna sequences resulted in a partial assembly of a potato virus y genome. two contigs were assembled, one encoding most of the polyprotein region and another encoding part of the 3′end region of the genome. the shorter contig included 184 bp at the 5′-end with no sequence similarity to the terminal region of the pvy polyprotein, instead reverse complementation of these 184 bp produced a match to the beginning of the pvy genome. in contrast to these findings, de novo assembly using only a 21–22 nt subset of small rnas resulted in the complete assembly of the pvy sequence including and extra bases in their 5′-end and 3′-end terminal regions, respectively. these findings suggest that improved and accurate assembly results are obtained using specific 21–22 nt small rnas produced by the plant antiviral response pathway.fig. assembly of a potato virus y genome using 21–25 nt and 21–22 nt reads comparison of toolkit outcomes with virfind we subjected raw rna-seq datasets to evaluate the performance of vsd toolkit as compared to virfind. summary of results with the virfind pipeline are outlined in additional file virfind identified viral sequences in all samples, of these in samples the vsd toolkit found the same viral sequences, while samples showed viral hits unique to the virfind pipeline. inspection of the later samples showed hits to viral sequences from plants, algae, fungi, invertebrate and vertebrates hosts. we aligned the adaptor clipped reads onto the assembled virfind contigs for these samples and found reads from 19 nt to up to 29 nt aligned onto contigs for these samples highlighting the broad spectrum of reads that are used in the virfind pipeline. to test if these unique hits in the virfind pipeline may potentially be chimeric sequences with host plant sequences, we conducted a blastx screening against a combined plant and virus protein database and reported the top hit. additional file shows that % of the virfind contigs annotated as viral sequences have a top hit to a plant or an insect virus. these findings highlight further the susceptibility to assemble chimeric sequences and/or plant sequences with similarity to viruses when using a broad spectrum of small rna read lengths collected from plant hosts. comparison of the viral sequences commonly assembled by virfind and our approach showed that the vsd toolkit assembled an average of % and % more of the viral sequences found in plant samples using the 21–25 nt and 21–22 nt workflows, respectively. the largest increases in the assembly of the known viral sequences were for the citrus endogenous pararetrovirus, grapevine freck virus, potato virus y and prunus necrotic ringspot virus rna. only in the case of the grapevine leafroll-associated virus did virfind assemble a greater proportion of the sequences than the vsd toolkit. both the vsd toolkit and virfind detected similarly all three viroids sequences. in the case of the citrus medica endogenous pararetrovirus the vsd toolkit assembled %, % and % of the genome using 21–25 nt, 21–22 nt and 24 nt assembly workflows, respectively. in contrast, virfind only assembled % of the cmeeprv genome sequence. overall we show that vsd toolkit produced improved and more accurate results than virfind. it is critical to utilise the specific by-products of the plant immune defence pathway against viruses and viroids to assemble highly accurate viral sequences that reflect the active host antiviral response. detection of exotic viral pathogens at the border is critical to safeguard plant industries and their access to national and international trade markets.  we have implemented an automated viral surveillance and diagnosis toolkit using the yabi web-based analytical environment that enables improved detection of viruses and viroids pathogens. we found that all single stranded rna viruses found in this study were assembled using exclusively 21–22 nt small rnas, while viroids were equally assembled using 21–22 nt or 24 nt subsets. the use of specific small rna subsets increases specificity of the identified viral sequences as it reflects the antiviral response activity of the plant hosts. additional files additional file 1: adaptor clipped and quality trimmed read statistics of small rna-seq datasets. clean reads = adaptor clipped small rna reads provided by bgi service provider; qc = quality control consisting of adaptor clipping and poor base call trimming; n.a. = not available. additional file 2: yabi viral surveillance and diagnosis toolkit user-guide. additional file 3: boxplot comparison of spades assemblies using individual kmers, sets of kmers and cap scaffolding. additional file 4: paired sample t-test of contigs produced by spades, velvet, abyss and soapdenovo assemblers. our optimised ‘spades k15-17-19_cap3’ assembly was compared against all other kmer settings and/or tools. a) spades - longest assembled contigs; b) spades - total number of assembled bases for contigs > =100 nt; c) spades - longest assembled contigs; d) spades - total number of assembled bases for contigs > =100 nt; e) comparison of longest assembled contig using velvet against sp22_k15-17-19_cap3; f) comparison of longest assembled contig using abyss against sp22_k15-17-19_cap3; g) comparison of longest assembled contig using soapdenovo against sp22_k15-17-19_cap additional file 5: blastn and blastx sequence similarity statistics for de novo assembled contigs using the 21–22 nt, 21–25 nt and 24 nt pipelines, respectively. twelve small rna-seq samples generated in this study with nucleotide similarity to viral pathogens were compared. see additional file for details of the selected samples. additional file 6: statistics of the viral sequences assembled using 21–25 nt, 21–22 nt and 24 nt small rna subsets. *virus name: grspav = grapevine rupestris stem pitting-associated virus; gvb = grapevine virus b; gfkv = grapevine fleck virus; glrav-3 = grapevine leafroll-associated virus; hsvd = hop stunt viroid; pnrsv = prunus necrotic ringspot virus; rbdv = raspberry bushy dwarf virus; rynv = rubus yellow net virus; spsmv-1 = sweetpotato symptomless mastrevirus 1; smov = strawberry mottle virus; msimv = miscanthus sinensis mosaic virus; ctv = citrus tristeza virus; cvd-vi = citrus viroid vi; mcdv = maize chlorotic dwarf virus; acvb = actinidia virus b; tvcv = tobacco vein-clearing virus. **virus type: rt = retro-transcribing virus; ssrna + = ssrna positive-strand virus; ssdna = single-stranded dna virus; dsrna = double-strand rna virus; v = viroid. additional file 7: mapping of high quality adaptor-clipped and quality trimmed small rnas with no mismatches onto the reference genomes: a) citrus tristesa virus. domains of the ctv genomes are denoted. b) hop stunt viroid. additional file 8: predicted rna secondary structure of viroids found in this study. minimal free energy rna secondary structure encoding base-pair probabilities are shown for s1_hsvd, s10_cvd-vi, and s10_hsvd. additional file 9: distribution of mapped and unmapped small rna reads for quarantined samples generated in this study. a) percentage of unmapped and mapped reads onto de novo assembled contigs with sequence similarity to plants; viruses, viroids and unknown are shown. b) distribution of the percentage of mapped reads for each sample along with the total number of distinct viral sequences are shown. details for each sample can be found in additional file additional file 10: de novo assembly of pnrsv rna rna and rna viral sequences using five randomly generated subsets of 1 m, 2 m, 5 m, 10 m and all adaptor clipped small rna reads. additional file 11: summary of virfind screening of small rna-seq datasets collected from quarantined plants. a) virfind blastn results. identified viruses/viroids were also identified using our vsd toolkit. b) virfind blastx results. identified viruses/viroids were only found using the virfind pipeline. hit to viruses from vertebrates, invertebrates, fungi or algae are shown in red background. additional file 12: virfind blastx results for contigs with sequence simalrity to viruses or viroids. additional file 13: mapping of adaptor clipped reads onto virfind assembled contigs for samples s s s s s and s details of samples can be found in additional file additional file 14: virfind contigs annotated with hits to viruses were run through the yabi vsd toolkit pipeline. top blastx hits to a combined plant and virus protein database are reported. virfind assembled contigs with a top hit to a plant virus are highlighted in green background. additional file 15: comparison of genome coverage of known viral sequences by contigs assembled by the vsd toolkit and virfind analysis pipelines. abbreviations abyssa parallel assembler for short read sequence data acv arisotelia chilensis virus acvbactinidia virus b awsamazon web services blastbasic local alignment search tool condetria content dependent read trimmer for illumina data ctv citrus tristeza virus cvd-vicitrus viroid vi ec2elastic compute cloud elisaenzyme-linked immunosorbent assay eprvendogenous pararetrovirus sequences gfkvgrapevine fleck virus glrav-3grapevine leafroll-associated virus grspavgrapevine rupestris stem pitting-associated virus gvbgrapevine virus b hsvdhop stunt viroid k15de novo assembly using kmer size of 15 k15-17-19de novo assembly using a combined kmer set of and 19 k15-17-19-21de novo assembly using a combined kmer set of and 21 k17de novo assembly using kmer size of 17 k19de novo assembly using kmer size of 19 k21de novo assembly using kmer size of 21 mcdvmaize chlorotic dwarf virus ngsnext generation sequencing pcrpolymerase chain reaction pcv persimmon cryptic virus peqpost entry quarantine pnrsvprunus necrotic ringspot virus pvypotato virus y rbdvraspberry bushy dwarf virus rna-seqrna sequencing rynvrubus yellow net virus samsequence alignment/map format sirnasmall interference small rna smovstrawberry mottle virus soapdenovoshort oligonucleotide analysis package for short read de novo assembly spadesa genome assembly algorithm spsmv-1sweetpotato symptomless mastrevirus 1 srmvsorghum mosaic virus ssrnasingle-strand rna tvcvtobacco vein-clearing virus velvetalgorithm package for de novo genome assembly and short read alignment virfindvirus find pipeline virnavirus-derived sirna vsdviral surveillance and diagnosis toolkit yabiyet another bioinformatics interface 
many studies reported the use of a broad range of small rnas encompassing the product sizes of several dicer enzymes involved in distinct biological pathways. we present a fully automated viral surveillance and diagnosis web-based bioinformatics toolkit that provides a flexible, user-friendly, robust and scalable interface for the discovery and diagnosis of viral pathogens. we have implemented an automated viral surveillance and diagnosis bioinformatics toolkit that produces improved viruses and viroid sequence assemblies. the vsd toolkit provides several optimised and reusable workflows applicable to distinct viral pathogens. optimised analysis workflows were customised within the yabi web-based analytical environment. existing post entry quarantine screening protocols rely on time-consuming biological indicators and/or molecular assays that require knowledge of infecting viral pathogens. we envisage that this resource will facilitate the surveillance and diagnosis viral pathogens in plants, insects and invertebrates. most of the viruses were assembled using only 21–22 nt long virus-derived sirnas, except for one citrus endogenous pararetrovirus that was more efficiently assembled using 24 nt long virnas. among the selected samples, we identified contigs with sequence similarity to viruses and viroids in samples. we also found that de novo assembly using 21–25 nt small rnas can result in chimeric assemblies of viral sequences and plant host sequences. plants have developed the ability to recognise and respond to viral infections through dicer-like enzymes that cleave viral sequences into specific small rna products. here we optimise the assembly of viral sequences by using specific small rna subsets. such non-specific assemblies can be resolved by using 21–22 nt or 24 nt small rnas subsets. we sequenced the small rna fractions of plants held at quarantine glasshouse facilities in australia and new zealand.
 detection and preventing entry of exotic viruses and viroids at the border is critical for protecting plant industries trade worldwide. existing post entry quarantine screening protocols rely on time-consuming biological indicators and/or molecular assays that require knowledge of infecting viral pathogens. plants have developed the ability to recognise and respond to viral infections through dicer-like enzymes that cleave viral sequences into specific small rna products. many studies reported the use of a broad range of small rnas encompassing the product sizes of several dicer enzymes involved in distinct biological pathways. here we optimise the assembly of viral sequences by using specific small rna subsets.  we sequenced the small rna fractions of plants held at quarantine glasshouse facilities in australia and new zealand. benchmarking of several de novo assembler tools yielded spades using a kmer of to produce the best assembly outcomes. we also found that de novo assembly using 21–25 nt small rnas can result in chimeric assemblies of viral sequences and plant host sequences. such non-specific assemblies can be resolved by using 21–22 nt or 24 nt small rnas subsets. among the selected samples, we identified contigs with sequence similarity to viruses and viroids in samples. most of the viruses were assembled using only 21–22 nt long virus-derived sirnas, except for one citrus endogenous pararetrovirus that was more efficiently assembled using 24 nt long virnas. all three viroids found in this study were fully assembled using either 21–22 nt or 24 nt virnas. optimised analysis workflows were customised within the yabi web-based analytical environment. we present a fully automated viral surveillance and diagnosis web-based bioinformatics toolkit that provides a flexible, user-friendly, robust and scalable interface for the discovery and diagnosis of viral pathogens.  we have implemented an automated viral surveillance and diagnosis bioinformatics toolkit that produces improved viruses and viroid sequence assemblies. the vsd toolkit provides several optimised and reusable workflows applicable to distinct viral pathogens. we envisage that this resource will facilitate the surveillance and diagnosis viral pathogens in plants, insects and invertebrates. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords bioinformaticsplant biosecuritynext generation sequencingplant viruses and viroidsquarantinevirnasvirus diagnosisyabismall rna-seqworkflowsplant biosecurity cooperative research centrecrc 2064barrero roberto a. issue-copyright-statement© the author 2017 