Which aligner for MAGs

[SilasK]

I wanted to change the aligner used in atlas recently.

For a long time, there was an error in bbmap which let it ignore the minimum identity threshold.

I thought usign minimap would be a good idea. https://lh3.github.io/2018/04/02/minimap2-and-the-future-of-bwa

I implemented minimap mapping for contigs #577 and genomes #578

When tested I see that now only 50% of reads are mapped to the genomes using the test data, where as almost all reads got mapped using bbsplit. The test data is a simulated metagenome from the 3 small genomes. The mutation rate used to generate the reads could be too high. There is als no parameter to modify the alignment with minimap2 e.g. min id.

Does someone want to compare the two mappers on real data?

@Waschina or @jmtsuji Do you have any idea?

[SilasK]

I try also bwa mem2 #582
It seems the main difference bwa is local and minimap global

[SilasK]

The problem is more with pileup!

Is this an error of counting paired-end reads only once or is there a threshold for why the reads are not counted with pileup?

Here is the log from pileup and the output from stamtool stats.

pileup output:

java -ea -Xmx51G -cp /gpfs/nihs/microbiome/Software/anaconda/Anaconda3-ATLAS/ATLAS_DATABASES/conda_envs/fe227e748e1ba3565890f0e6ee62c2d6_/opt/bbmap-39.01-0/current/ jgi.CoveragePileup in=genomes/alignments/sample1.bam threads=8 -Xmx51G covstats=genomes/alignments/coverage/sample1.tsv.gz fastaorf=genomes/all_orfs.fasta outorf=genomes/alignments/orf_coverage/sample1.tsv.gz concise=t physical=t bincov=genomes/alignments/coverage_binned/sample1.tsv.gz
Executing jgi.CoveragePileup [in=genomes/alignments/sample1.bam, threads=8, -Xmx51G, covstats=genomes/alignments/coverage/sample1.tsv.gz, fastaorf=genomes/all_orfs.fasta, outorf=genomes/alignments/orf_coverage/sample1.tsv.gz, concise=t, physical=t, bincov=genomes/alignments/coverage_binned/sample1.tsv.gz]

Set threads to 8
COUNT_GC disabled for physical coverage mode.
USE_SECONDARY disabled for physical coverage mode.
Found samtools 1.16.1
Could not find sambamba.
Reads: 1030435
Mapped reads: 518947
Mapped bases: 56689415
Ref scaffolds: 190
Ref bases: 3428138

Percent mapped: 50.362
Percent proper pairs: 49.608
Average coverage: 16.537
Average coverage with deletions: 32.886
Standard deviation: 55.444
Percent scaffolds with any coverage: 33.68
Percent of reference bases covered: 76.15

Time: 6.598 seconds.

samtools stats output:

reads mapped and paired: 1026640 # paired-end technology bit set + both mates mapped
reads unmapped: 28131
reads properly paired: 1022258 # proper-pair bit set
reads paired: 1058008 # paired-end technology bit set
reads duplicated: 0 # PCR or optical duplicate bit set
reads MQ0: 5221 # mapped and MQ=0
reads QC failed: 0
non-primary alignments: 0
supplementary alignments: 558
total length: 116150988 # ignores clipping
total first fragment length: 58071722 # ignores clipping
total last fragment length: 58079266 # ignores clipping
bases mapped: 113064706 # ignores clipping
bases mapped (cigar): 112666387 # more accurate
bases trimmed: 0
bases duplicated: 0
mismatches: 126624 # from NM fields
error rate: 1.123884e-03 # mismatches / bases mapped (cigar)
average length: 109
average first fragment length: 110
average last fragment length: 110
maximum length: 120
maximum first fragment length: 120
maximum last fragment length: 120
average quality: 33.5
insert size average: 218.1
insert size standard deviation: 41.9
inward oriented pairs: 510394
outward oriented pairs: 815
pairs with other orientation: 98
pairs on different chromosomes: 2013
percentage of properly paired reads (%): 96.6

[SilasK]

The problem is probably not at all as serious as I thought. I use pileup with physical coverage. E.g. the bases between the two read pairs are counted as covered. However the result is that only one fragment per properly paired input reads are counted.

[Waschina]

Hi @SilasK

I just wanted to test the minimap2 solution, but the new rules refer to scripts from the github snakemake-wrapper repository. The computation nodes in our cluster do not have internet access so the rules index_genomes and align_reads_to_genomes do fail :(

e.g. from the cluster log:

Failed to open source file https://github.com/snakemake/snakemake-wrappers/raw/v1.19.0/bio/minimap2/aligner/environment.yaml
ConnectionError: HTTPSConnectionPool(host='github.com', port=443): Max retries exceeded with url: /snakemake/snakemake-wrappers/raw/v1.19.0/bio/minimap2/aligner/environment.yaml (Caused by NewConnectionError('<urllib3.connection.HTTPSConnection object at 0x2b589c337190>: Failed to establish a new connection: [Errno 101] Network is unreachable')), attempt 1/3 failed - retrying in 3 seconds...

[Waschina]

Quick follow-up question: At which stage are wrappers installed? I can't find the snakemake-wrappers listed in any of environment yamls.

[SilasK]

Actually, I think using a dev atlas where it is installed with pip --editable is not a bad idea. The only thing that could be outdated is the tools in the atlasenv.yml especially snakemake (there are so many new versions all the time).

I think it's odd that you got the connection error in the cluster log. Could you try it with the latest snakemake version? and if the error persists we should raise an issue at the snakemake page.

In theory, the wrappers are installed like any other conda package by the main snakemake. Your error refers to the fact that snakemake tries to access the conda.yml on GitHub to install it. However, each job executed on the cluster is a new snakemake run, which would install all dependencies.

[Waschina]

Hi @SilasK

Sorry for getting back to this only now... We have tried reinstalling the atlas dev version from github and tried with the latest snakemake version. We still have this issue, and it occurs also at other rules that refer to snakemake-wrappers (e.g. align_reads_to_contigs in the assembly workflow).

This seems to be a known issue in snakemake for a while now: snakemake/snakemake#1157

I'll try with the "--wrapper-prefix" flag that you mentioned earlier.

[Waschina]

Using the --wrapper-prefix 'git+file://path/to/your/local/clone' worked, thank you!

It would be better if snakemake handled the wrappers and the installation of the environments in the same way, regardless of whether the wrapper repository is stored locally or on Github. I'll try to come up with a minimal example and post it to the snakemake github repo.

[SilasK]

Thank you for following up. I find it also anoying to use the best way by using the wrappers and then having more problems. Feel free do tag me if you have a snakemake issue.