-
Notifications
You must be signed in to change notification settings - Fork 34
Metagenomic Co Assembly Proof of Concept
By Tomer Altman, Connor Morgan-Lang, and Ryan McLaughlin**
(** Big thanks to Ryan for doing a ton of work to put all the data together!)
We set out to demonstrate that when you co-assemble read sets from many samples from a common source, you can recover more complete genomes, and genomes of rare members of the community. To this end, we fetched all 50 runs of PRJNA602689, which contains runs SRR10951656 and SRR10951660; the "low coverage Infectious Bronchitis Virus in Sus scrofa" study. So while we did not perform assembly for all of the benchmark datasets, we zeroed in on a particular use case to see if we could recover additional viral contigs. Below we present our approach, our findings, and instructions on fetching the data resources.
Samples were downloaded from the SRA using the prefetch (v2.8.0) and fastq files extracted using fastq-dump (v2.8.0). The Sus scrofa reference genome (GCF_000003025.6) was downloaded from NCBI. The Sus scrofa genome was indexed for bmfilter (v.3.102.4) and srprism (2.4.24-alpha) and a blastdb was created for blastn (2.5.0+), all for use by bmtagger (1.1.0). Paired-end reads for all samples were run through bmtagger to remove contaminating reads from the host (Sus s.). Blacklisted reads were filtered from all fastqs using seqmagick (v0.8.0). Reads were trimmed for quality using bbduk (v38.79). Finally all QC'ed samples were concatenated then co-assembled with MegaHIT (v1.2.9). Binning was preformed on contigs from the co-assembly with MetaBat2 (v2.12.1). Bin quality, presence and quality of recovered viral contigs, was assessed using CheckV (v0.6.0). We ran NCBI Blast in blastn mode against a local copy of the NT database, restricting the search to only viral entries, to annotate the contigs.
The entire workflow except for the Blast annotation was performed at a UBC computing grid, on a machine with 8 cores and approximately 64 GB RAM.
Read filtering and trimming took approximately 45 minutes and <16 GB RAM. The combined QC'ed read pairs numbered 30,152,318, and took up 15.3 GB on disk. It took 127 minutes to run the co-assembly, with an approximate RAM utilization of 53 GB and resulted in XXXX contigs of length greater than XXXX. The binning took 20 minutes and generated XXXX bins. Running CheckV on all of the bins took <1 minute and <1 GB RAM per bin.
Quality metric output from MegaHIT for the co-assembly:
2395 contigs, total 1818119 bp, min 500 bp, max 15999 bp, avg 759 bp, N50 702 bp
Time elapsed: 2346.165702 seconds
Quality metric output from MegaHIT for run SRR10951656:
2020-06-04 06:37:36 - 62 contigs, total 48807 bp, min 501 bp, max 2010 bp, avg 787 bp, N50 760 bp
2020-06-04 06:37:36 - ALL DONE. Time elapsed: 25.971600 seconds
Quality metric output from MegaHIT for run SRR10951660:
2020-06-04 06:39:34 - 49 contigs, total 38468 bp, min 502 bp, max 3121 bp, avg 785 bp, N50 783 bp
2020-06-04 06:39:34 - ALL DONE. Time elapsed: 28.095327 seconds
The following table lists the name of the organism, and the cumulative length of contigs with annotations to that organism, for organisms with a cumulative length greater than 1,000:
Infectious bronchitis virus 27753
Dickeya phage phiDP10.3 15708
Avian avulavirus 1 15205
Porcine astrovirus 4 10996
Moraxella phage Mcat16 8915
uncultured human fecal virus 8261
Mamastrovirus 3 6978
Pseudomonas phage PN05 4099
Epstein-Barr virus 2652
Pseudomonas phage vB_Pae_CF53a 1336
Porcine bocavirus 1221
Porcine rotavirus 1172
We discuss a select number of these sequences below.
We recovered a bin of IBV with 27,753 bp across three contigs of length 15,999, 8,971, 2,783, respectively. While the right target size for a coronavirus, CheckV declared the largest contig of medium-quality, and the other contigs of low quality, despite the Blast results showing high identity:
# Fields: query acc.ver, subject acc.ver, subject tax id, subject com names, query length, subject length, alignment length, % identity, % query coverage per subject, bit score, evalue
k119_5330 MH878976.1 11120 Infectious bronchitis virus 2783 27467 2783 99.928 100 5011 0.0
k119_5774 FJ904716.1 11120 Infectious bronchitis virus 15999 27629 16002 93.670 100 24288 0.0
k119_6018 MH878976.1 11120 Infectious bronchitis virus 8971 27467 8885 99.989 99 16019 0.0
CheckV flagged the following bins as being of high quality:
checkv_bin_12/quality_summary.tsv:k119_5808 15205 1.0 6 6 0 High-quality High-quality 99.74 AAI-based 0.0 No
checkv_bin_17/quality_summary.tsv:k119_6451 6458 1.0 3 3 0 High-quality High-quality 99.35 AAI-based 0.0 No
checkv_bin_47/quality_summary.tsv:k119_5360 6670 1.0 3 3 0 High-quality High-quality 100.0 AAI-based 0.0 No
Looking up their annotation from Blast:
# Fields: query acc.ver, subject acc.ver, subject tax id, subject com names, query length, subject length, alignment length, % identity, % query coverage per subject, bit score, evalue
k119_5808 MH996950.1 11176 Avian avulavirus 1 15205 15147 15138 99.908 99 27237 0.0
k119_6451 KP747574.1 1239567 Mamastrovirus 3 6458 6402 6400 95.828 99 10356 0.0
k119_5360 MK613068.1 1105379 Porcine astrovirus 4 6670 6733 6770 77.947 99 5383 0.0
So we recovered two high-quality genomes, one for Avian avulavirus 1 and one for Mamastrovirus 3, which have high identity to the reference database. And we recovered a third high-quality genome, Porcine astrovirus 4, which has 77.95 % identity to the closest reference sequence in the database.