Metagenome assembly and binning¶
In this tutorial you'll learn how to inspect assemble metagenomic data and retrieve draft genomes from assembled metagenomes
We'll use a mock community of 20 bacteria sequenced using the Illumina HiSeq. In reality the data were simulated using InSilicoSeq.
The 20 bacteria in the dataset were selected from the Tara Ocean study that recovered 957 distinct Metagenome-assembled-genomes (or MAGs) that were previsouly unknown! (full list on figshare )
Getting the Data¶
mkdir -p ~/data
cd ~/data
curl -O -J -L https://osf.io/th9z6/download
curl -O -J -L https://osf.io/k6vme/download
chmod -w tara_reads_R*
Quality Control¶
we'll use FastQC to check the quality of our data, as well as sickle for trimming the bad quality part of the reads. If you need a refresher on how and why to check the quality of sequence data, please check the Quality Control and Trimming tutorial
mkdir -p ~/results
cd ~/results
ln -s ~/data/tara_reads_* .
fastqc tara_reads_*.fastq.gz
Question
What is the average read length? The average quality?
Question
Compared to single genome sequencing, what graphs differ?
Now we'll trim the reads using sickle
sickle pe -f tara_reads_R1.fastq.gz -r tara_reads_R2.fastq.gz -t sanger \ -o tara_trimmed_R1.fastq -p tara_trimmed_R2.fastq -s /dev/null
Question
How many reads were trimmed?
Assembly¶
Megahit will be used for the assembly.
megahit -1 tara_trimmed_R1.fastq -2 tara_trimmed_R2.fastq -o tara_assembly
the resulting assenmbly can be found under tara_assembly/final.contigs.fa
.
Question
How many contigs does this assembly contain?
Binning¶
First we need to map the reads back against the assembly to get coverage information
ln -s tara_assembly/final.contigs.fa . bowtie2-build final.contigs.fa final.contigs bowtie2 -x final.contigs -1 tara_reads_R1.fastq.gz -2 tara_reads_R2.fastq.gz | \ samtools view -bS -o tara_to_sort.bam samtools sort tara_to_sort.bam -o tara.bam samtools index tara.bam
then we run metabat
runMetaBat.sh -m 1500 final.contigs.fa tara.bam
mv final.contigs.fa.metabat-bins1500 metabat
Question
How many bins did we obtain?
Checking the quality of the bins¶
The first time you run checkm
you have to create the database
sudo checkm data setRoot ~/.local/data/checkm
checkm lineage_wf -x fa metabat checkm/ checkm bin_qa_plot -x fa checkm metabat plots
Question
Which bins should we keep for downstream analysis?
Note
checkm can plot a lot of metrics. If you have time, check the manual and try to produce different plots
Warning
if checkm fails at the phylogeny step, it is likely that your vm doesn't have enough RAM. pplacer requires about 35G of RAM to place the bins in the tree of life.
In that case, execute the following
cd ~/results
curl -O -J -L https://osf.io/xuzhn/download
tar xzf checkm.tar.gz
checkm qa checkm/lineage.ms checkm
then plot the completeness
checkm bin_qa_plot -x fa checkm metabat plots
and take a look at plots/bin_qa_plot.png