Trimming reads using ivar - assembly of sars-cov-2 genome

Question

[fhsantanna]

I am following the workflow of Erin Young (https://github.com/CDCgov/SARS-CoV-2_Sequencing/tree/master/protocols/BFX-UT_ARTIC_Illumina) for assembly of the genomes of SARS-Cov-2. For sequencing the genomes, we utilized the QIAseq Library kit, which is based on ARTIC workflow, and the Illumina platform. One important step of the QIAseq kit is the fragmentation before the adapter ligation.

In the bioinformatic workflow, after mapping the reads to the reference genome, we have to trim the primers from the aligned reads to account for the true variation of the virus, and for this process we utilized ivar.

I read that it is recommended to include reads that do not present primers in their sequence, because of the fragmentation step (https://covid19.galaxyproject.org/artic/#a-galaxy-workflow-for-the-analysis-of-illumina-paired-end-sequenced-artic-amplicon-data). The problem is that when i utilize the "-e option" of ivar, which includes reads without primer, the consensus sequence presents more Ns than that generated without the "-e option".

Why is it happening? Since more reads are being included, I wondered that I would have a better consensus.

Here are the commands:

bwa mem -t 6 ../genome/mn908947.fasta ../../fastq/0044_L001_ds.bf5a13d8963047738a70006321c2778e/0044_S8_L001_R1_001.fastq.gz ../../fastq/0044_L001_ds.bf5a13d8963047738a70006321c2778e/0044_S8_L001_R2_001.fastq.gz | samtools sort | samtools view -F 4 -o 0044.sorted.eopt.bam

ivar trim -e -i 0044.sorted.eopt.bam -b ../artic_primers/nCoV-2019.primer.bed -p 0044.primertrim

samtools sort 0044.primertrim.bam  -o 0044.primertrim.sorted.eopt.bam 

samtools mpileup -A -d 6000000 -B -Q 0 --reference ../genome/mn908947.fasta 0044.primertrim.sorted.eopt.bam  | ivar consensus -p 0044_consensus.fas -n N

Here is the comparison of the coverages:

#reads without primers excluded
MN908947.3 (29.9Kbp)
>  90.00% │ ██████████████████████████████████████████████████████████████████████████████████  ██████████████████████████████████ ████████ │ Number of reads: 426385
>  80.00% │ ██████████████████████████████████████████████████████████████████████████████████▁ ██████████████████████████████████ ████████▆│ 
>  70.00% │▅███████████████████████████████████████████████████████████████████████████████████▅██████████████████████████████████▂█████████│ Covered bases:   29.6Kbp
>  60.00% │█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████│ Percent covered: 98.98%
>  50.00% │█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████│ Mean coverage:   1.73e+03x
>  40.00% │█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████│ Mean baseQ:      37.2
>  30.00% │█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████│ Mean mapQ:       60
>  20.00% │█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████│ 
>  10.00% │█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████│ Histo bin width: 231bp
>   0.00% │█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████│ Histo max bin:   100%
          1        2.3K      4.6K      6.9K      9.2K     11.6K     13.9K     16.2K     18.5K     20.8K     23.1K     25.4K              29.9K  


#e option
MN908947.3 (29.9Kbp)
>  90.00% │ ██████████████████████████████████████████████████████████████████████████████████ █████████████████████████████████████████████│ Number of reads: 673402
>  80.00% │▅██████████████████████████████████████████████████████████████████████████████████▅█████████████████████████████████████████████│ 
>  70.00% │█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████│ Covered bases:   29.8Kbp
>  60.00% │█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████│ Percent covered: 99.68%
>  50.00% │█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████│ Mean coverage:   2.76e+03x
>  40.00% │█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████│ Mean baseQ:      37.2
>  30.00% │█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████│ Mean mapQ:       60
>  20.00% │█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████│ 
>  10.00% │█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████│ Histo bin width: 231bp
>   0.00% │█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████│ Histo max bin:   100%
          1        2.3K      4.6K      6.9K      9.2K     11.6K     13.9K     16.2K     18.5K     20.8K     23.1K     25.4K              29.9K

Here are the consensus files: https://we.tl/t-Cdxf1Zimuf

[GenoMax]

iVar uses primer positions supplied in a BED file to soft clip primer sequences from an aligned and sorted BAM file.

So in theory it should be doing nothing to reads that don't have primers? So whatever is happening is because of some other program?

[fhsantanna]

Yes, if I utilize the -e option reads that don't have primers are kept. Since amplicons are fragmented during QIASeq procedure, some DNAs of the library pool probably will not have primer sites. Reads without primers are removed by default by ivar. The only program that could influence the procedure is samtools.

[Arup Ghosh]

Inspect the ivar processed bam files in IGV or other genome browsers and also check ivar the log file.

[fhsantanna]

This was actually a good idea. I identified a region that was different between the assembly approaches. This region presents a depth mode of 3X, it was included in the original protocol, but removed when I utilized the -e option of ivar. Interestingly, the first two bases of this region presents a depth of 69 X and 15X in the orginal protocol, but 9x and 5X using the -e option.

Answer 1

[fhsantanna]

I found the problem. I was actually comparing results of different versions of ivar, 1.1 and 1.3. The last version was introducing more N's in the consensus, because it demands a depth of 10x for calling a base.