Describe the bug
Some read alignments are shifted after ivar trim
(v1.2) compared to untrimmed read alignments from Minimap2 (v2.17-r941) and NGMLR (v0.2.7).
To Reproduce
I aligned the Nanopore reads to the reference genome (MN908947.3) with Minimap2 (and NGMLR), trimmed primer sequences and low quality bases
$ minimap2 -ax map-ont -t16 ref.fa reads.fq | samtools sort -@16 | samtools view -F4 -b -o reads.sort.bam
$ ivar trim -b artic-ncov2019/primer_schemes/nCoV-2019/V2/nCoV-2019.bed -p reads.trim -i reads.sort.bam -q 1 -m 30 -s 4
Note: for some reason many of the reads had very low average quality scores for this run and we were mostly interested in trimming the primer sequences.
I viewed and compared alignments in IGV (v2.8.2 installed from BioConda) and noticed some big differences between the untrimmed and trimmed alignments of some reads in certain regions.
Top alignment is untrimmed Minimap2 and bottom is iVar trimmed alignment:
![Screenshot_20200415_121302-IGV-ivar-minimap2-comparison-multiple-reads](https://user-images.githubusercontent.com/4203424/79369712-04385f00-7f17-11ea-984e-cf3bc017db92.png)
For example, here's the same read at position MN908947.3:29,383 where in the original alignment the base is G (QV=27) and in the trimmed alignment the base is A (QV=11):
Before ivar trim
read alignment information in IGV
Read name = 04853829-6653-42bb-b303-351fbd160119
Read length = 515bp
----------------------
Mapping = Primary @ MAPQ 60
Reference span = MN908947.3:29,363-29,646 (+) = 284bp
Cigar = 89S15M1D11M1I3M1D20M1I14M5I3M2I2M2I18M1I33M1D7M1I7M2I12M1I28M2I8M5I3M1D12M2I3M4I12M1I6M2D5M1D23M1D31M120S
Clipping = Left 89 soft; Right 120 soft
----------------------
s1 = 79
s2 = 0
NM = 51
AS = 314
de = 0.1145
rl = 0
cm = 9
nn = 0
tp = P
ms = 314Location = MN908947.3:29,383
Base = G @ QV 27
After ivar trim
Read name = 04853829-6653-42bb-b303-351fbd160119
Read length = 515bp
----------------------
Mapping = Primary @ MAPQ 60
Reference span = MN908947.3:29,380-29,647 (+) = 268bp
Cigar = 104S11M1I3M1D20M1I14M5I3M2I2M2I18M1I33M1D7M1I7M2I12M1I28M2I8M5I3M1D12M2I3M4I12M1I6M2D5M1D23M1D31M120S
Clipping = Left 104 soft; Right 120 soft
----------------------
s1 = 79
s2 = 0
NM = 51
AS = 314
de = 0.1145
rl = 0
cm = 9
nn = 0
tp = P
ms = 314
Hidden tags: XALocation = MN908947.3:29,383
Base = A @ QV 11
It looks like the iVar trimmed alignment has been shifted when I zoom in on the read in both untrimmed and trimmed alignments:
![Screenshot_20200415_112320-IGV-ivar-minimap2-comparison](https://user-images.githubusercontent.com/4203424/79369646-e8cd5400-7f16-11ea-923e-a7dd38346748.png)
Note that the position is incremented by 1 in the trimmed read alignment relative to the untrimmed alignment.
Expected behavior
I expected that the alignments should remain the same aside from primer and quality trimming.
Workstation:
- OS: Arch Linux x86_64 (Kernel Release: 5.5.9-arch1-2)
- CPU: Intel(R) Core(TM) i9-9980HK CPU @ 2.40GHz
- iVar Version 1.2 from BioConda
Additional context
Conda env YAML ($ conda env export
) for running iVar, Minimap2, NGMLR and samtools:
channels:
- conda-forge
- bioconda
- defaults
dependencies:
- _libgcc_mutex=0.1=conda_forge
- _openmp_mutex=4.5=0_gnu
- bzip2=1.0.8=h516909a_2
- ca-certificates=2020.4.5.1=hecc5488_0
- curl=7.69.1=h33f0ec9_0
- htslib=1.9=h4da6232_3
- krb5=1.17.1=h2fd8d38_0
- libcurl=7.69.1=hf7181ac_0
- libdeflate=1.5=h516909a_0
- libedit=3.1.20170329=hf8c457e_1001
- libgcc-ng=9.2.0=h24d8f2e_2
- libgomp=9.2.0=h24d8f2e_2
- libssh2=1.8.2=h22169c7_2
- libstdcxx-ng=9.2.0=hdf63c60_2
- minimap2=2.17=h8b12597_1
- ncurses=6.1=hf484d3e_1002
- ngmlr=0.2.7=he860b03_1
- openssl=1.1.1f=h516909a_0
- samtools=1.9=h10a08f8_12
- sniffles=1.0.11=hdbcaa40_1
- tclap=1.2.1=h470a237_1
- tk=8.6.10=hed695b0_0
- xz=5.2.5=h516909a_0
- zlib=1.2.11=h516909a_1006