To investigate the possible aetiological agents associated with this disease, we collected bronchoalveolar lavage fluid (BALF) and performed deep meta-transcriptomic sequencing.
The clinical specimen was handled in a biosafety level 3 laboratory at Shanghai Public Health Clinical Center. Total RNA was extracted from 200 μl of BALF and a meta-transcriptomic library was constructed for pair-end (150-bp reads) sequencing using an Illumina MiniSeq as previously described.
In total, we generated 56,565,928 sequence reads that were de novo-assembled and screened for potential aetiological agents.
Of the 384,096 contigs assembled by Megahit, the longest (30,474 nucleotides (nt)) had a high abundance and was closely related to a bat SARS-like coronavirus (CoV) isolate—bat SL-CoVZC45 (GenBank accession number MG772933)—that had previously been sampled in China, with a nucleotide identity of 89.1% (Supplementary Tables 1, 2).
The genome sequence of this virus, as well as its termini, were determined and confirmed by reverse-transcription PCR (RT–PCR) and 5′/3′ rapid amplification of cDNA ends (RACE), respectively.
This virus strain was designated as WH-Human 1 coronavirus (WHCV) (and has also been referred to as ‘2019-nCoV’) and its whole genome sequence (29,903 nt) has been assigned GenBank accession number MN908947.
Remapping the RNA-sequencing data to the complete genome of WHCV resulted in an assembly of 123,613 reads, providing 99.99% genome coverage at a mean depth of 6.04× (range, 0.01–78.84×) (Extended Data Fig. 3).
The viral load in the BALF sample was estimated by qPCR to be 3.95 × 108 copies per ml (Extended Data Fig. 4).

The viral genome organization of WHCV was determined by sequence alignment to two representative members of the genus Betacoronavirus: a coronavirus associated with humans (SARS-CoV Tor2, GenBank accession number AY274119) and a coronavirus associated with bats (bat SL-CoVZC45, GenBank accession number MG772933).
The un-translational regions and open-reading frame (ORF) of WHCV were mapped on the basis of this sequence alignment and ORF prediction.
The WHCV viral genome was similar to these two coronaviruses (Fig. 1 and Supplementary Table 3).
The order of genes (5′ to 3′) was as follows: replicase ORF1ab, spike (S), envelope (E), membrane (M) and nucleocapsid (N).
WHCV has 5′ and 3′ terminal sequences that are typical of betacoronaviruses, with 265 nt at the 5′ terminal end and 229 nt at the 3′ terminal end.
The predicted replicase ORF1ab gene of WHCV is 21,291 nt in length and contained 16 predicted non-structural proteins (Supplementary Table 4), followed by (at least) 13 downstream ORFs.
Additionally, WHCV shares a highly conserved domain (LLRKNGNKG: amino acids 122–130) in nsp1 with SARS-CoV.
The predicted S, ORF3a, E, M and N genes of WHCV are 3,822, 828, 228, 669 and 1,260 nt in length, respectively.
In addition to these ORF regions, which are shared by all members of the subgenus Sarbecovirus, WHCV is similar to SARS-CoV in that it carries a predicted ORF8 gene (with a length of 366 nt) that is located between the M and N ORF genes.
The functions of WHCV ORFs were predicted on the basis of those of known coronaviruses and are described in Supplementary Table 5.
In a manner similar to SARS-CoV Tor2, a leader transcription regulatory sequence (TRS) and nine putative body TRSs could be readily identified upstream of the 5′ end of the ORF in WHCV, and the putative conserved TRS core sequence appeared in two forms—ACGAAC or CUAAAC (Supplementary Table 6).

Total RNA was extracted from the BALF sample using the RNeasy Plus Universal Mini kit (Qiagen) following the manufacturer’s instructions.
The quantity and quality of the RNA solution was assessed using a Qbit machine and an Agilent 2100 Bioanalyzer (Agilent Technologies) before library construction and sequencing.

An RNA library was then constructed using the SMARTer Stranded Total RNA-Seq kit v.2 (TaKaRa).
Ribosomal RNA depletion was performed during library construction following the manufacturer’s instructions.
Paired-end (150-bp reads) sequencing of the RNA library was performed on the MiniSeq platform (Illumina).
Library preparation and sequencing were carried out at the Shanghai Public Health Clinical Center, Fudan University, Shanghai, China.

Sequencing reads were first adaptor and quality trimmed using the Trimmomatic program.
The remaining 56,565,928 reads were assembled de novo using both Megahit (v.1.1.3) and Trinity (v.2.5.1) with default parameter settings.
Megahit generated a total of 384,096 assembled contigs (size range of 200–30,474 nt), whereas Trinity generated 1,329,960 contigs with a size range of 201–11,760 nt.

All of these assembled contigs were compared (using BLASTn and Diamond BLASTx) against the entire non-redundant (nr) nucleotide and protein databases, with e values set to 1 × 10^(−10) and 1 × 10^(−5), respectively.
To identify possible aetiological agents present in the sequencing data, the abundance of the assembled contigs was first evaluated as the expected counts using the RSEM program implemented in Trinity.
Non-human reads (23,712,657 reads), generated by filtering host reads using the human genome (human release 32, GRCh38.p13, downloaded from Gencode) by Bowtie2, were used for the RSEM abundance assessment.

As the longest contigs generated by Megahit (30,474 nt) and Trinity (11,760 nt) both showed high similarity to the bat SARS-like coronavirus isolate bat SL-CoVZC45 and were found at a high abundance (Supplementary Tables 1, 2), the longer sequence (30,474 nt)—which covered almost the whole virus genome—was used for primer design for PCR confirmation and determination of the genome termini.
Primers used for PCR, qPCR and RACE experiments are listed in Supplementary Table 8.
The PCR assay was conducted as previously described and the complete genome termini was determined using the Takara SMARTer RACE 5′/3′ kit (TaKaRa) following the manufacturer’s instructions.
Subsequently, the genome coverage and sequencing depth were determined by remapping all of the adaptor- and quality-trimmed reads to the whole genome of WHCV using Bowtie2 and Samtools.

Assemblers tested

In this study, de novo assembly of paired-end reads was performed using the current versions of eight different short-read assemblers.
We used ABySS assembler which is optimized for short reads.
The parallel version of ABySS is capable of assembling large genomes [21].
MEGAHIT is an ultra-fast and memory-efficient short-read assembler, optimized for metagenomes, also works well on generic single genome assembly of small or mammalian size [22].
Ray Meta is used for metagenome assembly and profiling [23].
SPAdes can assemble sequences from single-cell and multi-cell data types [24].
The Velvet assembler was designed for short-read sequencing data [25];

metaSPAdes is a metagenomic assembler and MetaVelvet is an extension of Velvet for metagenome assembly from short reads [26, 27].
Trinity performs de novo transcriptome assembly [28].
For every assembler mentioned above, we have used default parameters unless otherwise mentioned. K-mer lengths 21, 63, 99 were used for ABySS, Velvet, MetaVelvet and Ray Meta.
For the rest of the assemblers, default k-mer length was applied (Figure 1A, Supplementary Table S1, see Supplementary Data available online at http://bib.oxfordjournals.org/).

SARS-CoV-2 genome assembly contiguity breaks at the repeat region

We investigated the simple tandem repeat of 585 bp in the SARS-CoV-2 genome located at the 3′-end (MN908947.3:29 870–29 903).
Most of the assemblers failed to assemble the repeat region and more than 100 bp gap was created at the 3′-end in different assemblies.

DISCUSSION

In this study, we compared 16 assembler variations using eight de novo assemblers for the benchmarking of the genome assembly quality of the SARS-CoV-2 virus.
We observed two metagenomic assemblers, e.g. MEGAHIT and metaSPAdes outperformed other assemblers in regards to the genome fraction recovery, largest contig length, N50 length, NA50 length, L50 and LA50 contig number.
The fraction of genome recovery could be 10-folds different between assemblers, e.g. MEGAHIT (99%) versus MetaVelvet-K21 (10%).
Although all eight assemblers used the graph-based method for de novo assembly, the differences we observed are due to the variations in their implementation, error correction, quality thresholds and choice of other parameters.
Despite better performances by the two metagenomic assemblers, the entire viral genome was not assembled in most cases, especially at the termini of the genome.
Therefore, there is a need to develop newer assembly methods specially designed to assemble complete viral genomes.
The SARS-CoV-2 virus genome could also be assembled by aligning the reads to the reference genome or using a reference guided assembly.