Skip to content

Day 4

Jump to bottom
Mateu Menéndez-Serra edited this page Sep 24, 2024 · 19 revisions

In this tutorial, we will replicate some current analyses using the new tools we have developed. By the end of this tutorial, you will be able to:

  • Map ancient metagenomic reads to a database of reference genomes to understand the sample's composition.
  • Identify specific genomes in the ancient metagenome.
  • Authenticate genomes based on damage patterns.

Session 1: Setting the Stage

First, let's create the necessary folder structure in our working directory (wdir):

cd ~/course/wdir
mkdir -p day4
cd day4

ln -s ~/course/data/day4/ data

Now, activate the environment for today's session:

conda activate env4

Before moving forward, we need to install a missing package:

pip install tabview

If everything went as expected, we will have our system ready for today.

Session 2 | Extension, dereplication and mapping

Our approach for analysing ancient microbial data is a bit different from what is the norm in the field. After QC'ing our sequences we follow the following steps:

extension -> dereplication -> mapping -> reassign -> filtering -> damage estimation

One of the first steps we do is to extend the reads by doing very gentle assemblies, but before proceeding let's get some statistics for our fastQ files:

seqkit stats -j 5 -T data/fastq/PRI-TJPGK-CATN-160-162.fq.gz | csvtk -t pretty

ℹ️ seqkit, csvtk and taxonkit are very useful tools. Check them at https://bioinf.shenwei.me/

And extend the reads using bbmap. We will need to do two different steps, first calculate how much memory we need so the extension is deterministic and then do the extension using the estimates:

MEM=$(loglog.sh seed=1234 k=17 in=data/fastq/PRI-TJPGK-CATN-160-162.fq.gz ignorebadquality 2> >(grep Cardinality) \
    | awk -vP=0.5 -vB=16 -vH=3 '{{print int( (((B*H)/8)*$2)/P )}}')

tadpole.sh -Xmx10G \
    k=17 \
    in=data/fastq/PRI-TJPGK-CATN-160-162.fq.gz \
    out=PRI-TJPGK-CATN-160-162.extended.fq.gz \
    mode=extend \
    ibb=f \
    prefilter=0 \
    el=100 er=100 \
    threads=5 \
    overwrite=true \
    trimends=9 \
    ecc=f ecco=f \
    filtermem="${MEM}" \
    conservative=t \
    ignorebadquality

Get the stats of PRI-TJPGK-CATN-160-162.extended.fq.gz

Once the reads have been extended, we will dereplicate them:

seqkit rmdup -j 5 -s -o PRI-TJPGK-CATN-160-162.extended.derep.fq.gz PRI-TJPGK-CATN-160-162.extended.fq.gz

Get the stats of PRI-TJPGK-CATN-160-162.extended.derep.fq.gz

Before mapping, we will get the original reads using the extended-dereplicated as a guide. Although we could use this set of reads for mapping, at the moment we prefer not to (any ideas why?):

filterbyname.sh in=data/fastq/PRI-TJPGK-CATN-160-162.fq.gz out=PRI-TJPGK-CATN-160-162.mapping.fastq.gz names=PRI-TJPGK-CATN-160-162.extended.derep.fq.gz threads=5 overwrite=t include=t

So finally we have the reads we can use for mapping. Although we use Bowtie2 we follow a slightly different strategy for mapping. First let's activate the conda environment that we need:

conda activate mapping

And let's do the mapping:

bowtie2 -p 5 -k 100 \
    -N 1 -D 5 -R 1 -L 22 -i S,0,2.50 \
    --np 1 --mp "1,1" --rdg "0,1" --rfg "0,1" --score-min "L,0,-0.1" \
    -x data/db/aegenomics.db \
    -q PRI-TJPGK-CATN-160-162.mapping.fastq.gz --no-unal \
    | samtools view -F 4 -b \
    | samtools sort -@ 32 -m 8G -o PRI-TJPGK-CATN-160-162.sorted.bam

We can have a look at the specific parameters at the bowtie2 manual: https://bowtie-bio.sourceforge.net/bowtie2/manual.shtml


Expand for an explanation of the parameters 
There's a lot going on in this last command. 
First, with -k we ask for a maximum of 100 alignments per read. This is because many of those reads are going to be mapping equally well to multiple references. 

Seed: We will do a fast search by using -N 1 -D 5 -R 1 -L 22 -i S,0,2.50. -i S,0,2.5 sets the interval function f to f(x) = 1 + 2.5 * sqrt(x), where x is the read length, which set how many seed substrings will be generated.
Those are the preset parameters for a fast search in Bowtie2, but with the difference that we are using -N 1 to allow a mismatch in the seed to increase sensitivity.


The parameters --np 1 --mp "1,1" --rdg "0,1" --rfg "0,1" determine the reward/penalization values.
We can adjust the % identity in --score-min "L,0,-0.1", where the -0.1 == 90, 0.05 == 95%


This is just an example for the tutorial, in real life you might want to use:


recommended, good compromise between sensitivity, specificity and speed: -D 15 -R 2 -N 1 -L 22 -i S,1,1.15 --np 1 --mp "1,1" --rdg "0,1" --rfg "0,1" --score-min "L,0,-0.1"


faster: -D 10 -R 2 -N 1 -L 22 -i S,0,2.50 --np 1 --mp "1,1" --rdg "0,1" --rfg "0,1" --score-min "L,0,-0.1"


more sensitive, but ~10X slower: -D 15 -R 2 -N 1 -L 20 -i S,1,1.15  --np 1 --mp "1,1" --rdg "0,1" --rfg "0,1" --score-min "L,0,-0.1"


super sensitive, but ~20X slower (Reduce -L if you want to be more sensitive): -D 15 -R 3 -N 1 -L 20 -i S,1,0.5  --np 1 --mp "1,1" --rdg "0,1" --rfg "0,1" --score-min "L,0,-0.1"

Now we have the reads mapped and ready for the next step!

Session 3 | Reassign and Filtering

Now is time to process the BAM file, filter out references with uneven coverages and estimate multiple metrics, including taxonomic abundances.

First we will remove duplicates from the BAM files:

And we will use sambamba to remove duplicates. Sambamba has similar heuristics than Picard, but is way faster:

sambamba markdup -r -t 5 -p PRI-TJPGK-CATN-160-162.sorted.bam PRI-TJPGK-CATN-160-162.sorted.rmdup.bam

At this point, duplicates are defined by the alignment position of the reads.

Then we reassign reads to the reference they belong using an E-M algorithm that takes into account the alignment score.

filterBAM reassign \
    -r data/misc/aegenomics.db.len.map \
    --bam PRI-TJPGK-CATN-160-162.sorted.rmdup.bam \
    -t 5 \
    -i 0 \
    -m "8G" \
    -o PRI-TJPGK-CATN-160-162.reassigned.bam \
    -n 3

Expand for an explanation of the printed output 
Found 2,689 reference sequences
BAM file looks good.
Resolving multi-mapping reads...
Loading BAM file
IO Threads: 1 | Processing Threads: 5
Found 2,689 reference sequences
Removing references with less than 3 reads...
Kept 56,487,020 alignments

Keeping 2,623 references


(...)


Iter: 1 - R: (removed alignments on the iteration) 44,750,470 | U: (remaining unique mapping reads) 933,240 | NU: 1,441,255 | L: (remaining alingments) 10,274,358


(...)


Iter: 17 - R: (removed alignments on the iteration) 0 | U: (remaining unique mapping reads) 3,052,259 | NU: 3,054,882 | L: (remaining alingments) 3,052,403
::: No more alignments to remove. Stopping.
References: 1,321 | Reads: 3,052,331 | Alignments: 3,052,403
Unique mapping reads: 3,052,259 | Multimapping reads: 72


(...)


Removing references with less than 3...
Total refs/reads combination: 3,052,050
Total references: 1,052

Then we run the filtering step, which estimates several metrics for each reference in our BAM file and filters out those that do not meet the defined criteria:

filterBAM filter \
    --bam PRI-TJPGK-CATN-160-162.reassigned.bam \
    -N \
    -r data/misc/aegenomics.db.len.map \
    -A 92 \
    -a 94 \
    -n 100 \
    -b 0.75 \
    -B 0.01 \
    -t 5 \
    --sort-memory 8G \
    --include-low-detection \
    --stats PRI-TJPGK-CATN-160-162.stats.tsv.gz \
    --stats-filtered PRI-TJPGK-CATN-160-162.stats-filtered.tsv.gz \
    --bam-filtered PRI-TJPGK-CATN-160-162.filtered.bam

Expand for an explanation of the parameters 
There's a lot going on in this last command:
We use -N to sort the filtered BAM files by name, so it can be used by metaDMG. 
The `-r` specifies where we can find the non-concatenated length of the references for the abundance estimation. 
With `-A` we will only keep those reads with at least 92% ANI and `-a`  will keep those references with at least 94% of average ANI. 
The `-n` will keep only those references with at least 100 reads mapping.

And... what are the -b and -B accounting for?

Let's deactivate the mapping environment:

conda deactivate

Let's have a look at the results of the filtering:

zcat PRI-TJPGK-CATN-160-162.stats-filtered.tsv.gz \
  | csvtk cut -t -T -f "reference,n_reads,read_ani_mean,read_ani_std,coverage_mean,breadth,exp_breadth,breadth_exp_ratio,norm_entropy,norm_gini,cov_evenness,tax_abund_tad" \
  | csvtk grep -r -t -v -f reference -p _plas -p _mito \
  | csvtk sort -t -T -k "n_reads:Nr" \
  | tabview -

We will explore four different examples to understand the effect of each filter. We will get the references GCA_002781685.1, IMGVR_UViG_3300027782_000260 and GCA_014380485.1.

printf "GCA_002781685.1\nIMGVR_UViG_3300027782_000260\nGCA_014380485.1" > ref-list.txt
getRPercId --bam PRI-TJPGK-CATN-160-162.sorted.rmdup.bam --reference-list ref-list.txt --threads 5 --sort-memory 8G

Activate the environment with the necessary tools:

conda activate env2

And let's explore the coverage patterns:

bamcov -w 0 -m GCA_002781685.1.bam 
bamcov -w 0 -m GCA_014380485.1.bam
bamcov -w 0 -m IMGVR_UViG_3300027782_000260.bam

Session 4 | Damage

No that we have our set of refined references, is time to authenticate them, this means to look for damage patterns. We recently developed a new method that can estimate damage over thousands of taxa (LCA mode), references (local mode) or provide a global estimate (global mode).

Let's first activate the needed environment:

conda activate metaDMG

We will use the local mode approach which allows retrieving damage estimates for the references present in the bam file. For this, we first run the getdamage command. This will create the mismatch matrices which we be use as an input of the dfit function.

metaDMG-cpp getdamage PRI-TJPGK-CATN-160-162.filtered.bam \
    --out_prefix PRI-TJPGK-CATN-160-162 \
    --threads 5 \
    --min_length 30 \
    --print_length 25 \
    --run_mode 1

Next, we run the dfit function which performs numerical optimization of the deamination frequencies based on the mismatch matrix (PRI-TJPGK-CATN-160-162.bdamage.gz) to estimate four parameters: A,q,c,phi.

metaDMG-cpp dfit PRI-TJPGK-CATN-160-162.bdamage.gz \
    --threads 5 \
    --bam PRI-TJPGK-CATN-160-162.filtered.bam \
    --showfits 2 \
    --out_prefix PRI-TJPGK-CATN-160-162
conda deactivate metaDMG

Session 6 | Visualisation

Let's load the R environment and lunch a r terminal attached to the VSC workspace:

conda activate r

radian

We will open the R scripts in the folder ~/course/data/day4/src in VS Code and load the r extension.

Clone this wiki locally