Simulated taxa missing from metaDMG outputs

[libby-natola]

Hi metaDMG devs!

I've used metaDMG several times in the past year or so in conjunction with the KapCopenhagan holi pipeline to estimate damage and perform LCAs on empirical ancient metagenomic data with no issue.

Now, I'm trying to track some simulated reads through several different pipelines to check true/false positive and false negative rates. I've created some simulations of a mixture of eukaryotic/prokaryotic communities with and without damage using gargammel and they've run through the holi (and kraken, and malt) pipelines without issue, so I think the simulated files are okay. I've run into two problems running metaDMG that I'm hoping you might have insight on:

1. I'm running into this strange behavior with my files, wherein none of my animalia reads seem to be reported as hits in the metaDMG output (bacteria and viridiplantae look as expected). I can see from the read ids and reference sequence ids that reads from these taxa are definitely present in the input sam.gz files, and mapping to taxa I would expect from my bowtie2 steps, but none of these make it through to the metaDMG output. I've checked all the names, nodes, and accession2taxid files and the taxa that are missing are present in these files. I also reran an empirical sample with these taxa present through otherwise the same config and the animal taxa that go missing in the files I simulated are present in the empirical metaDMG output, which leads me to think there is nothing wrong with my system configuration or reference files.

2. I'd like to track the lca call for each read (eg what taxon was read X assigned?), but I can't seem to find the full information on that for every read in any of the output files. I'm sure this information must exist at some point for metaDMG to report the total number of reads hit per taxon, is there an argument that allows you to save temporary files that would give me this information? I expected to find this information in the data/lca/sample.lca.txt.gz files, but there are far fewer reads listed than exist in the input, so perhaps this issue is actually related to point number 1 and when metaDMG is not reporting the Animalia hits that's why I'm getting fallout of reads in the output?

I'm very grateful for any suggestions you have on how I can overcome these problems. We think metaDMG is integral to our best pipeline option and it would be great to confirm that with our simulations.

Cheers,
Libby

Below is the code I've used to run this, first unsuccessfully (expected taxa missing from output) with simulated data and then successfully with empirical data (expected taxa present in output):

`conda activate metaDMG
export LD_LIBRARY_PATH="$LD_LIBRARY_PATH:/usr/local/lib"
nam=/data/taxonomyDBs/NCBI_taxonomy/2025-01-30/names.dmp
nod=/data/taxonomyDBs/NCBI_taxonomy/2025-01-30/nodes.dmp
acc=/data/taxonomyDBs/NCBI_taxonomy/2025-01-30/nucl_gb.accession2taxid.gz

#unsuccessful
mkdir metaDMG/ && cd metaDMG
metaDMG config ../*/*merged.sam.gz --names $nam --nodes $nod --acc2tax $acc --metaDMG-cpp /data2/programs/metaDMG-cpp/metaDMG-cpp --config-file config.metaDMG.$run.yaml --bayesian --parallel-samples 1 --cores-per-sample $threads --custom-database
metaDMG compute config.metaDMG.$run.yaml
metaDMG convert --output $run.csv config.metaDMG.$run.yaml
metaDMG plot config.metaDMG.$run.yaml --output $run.plots.pdf

#successful
metaDMG config ../SC-10c-TCMP_S207_L001.merged.sam.gz --names $nam --nodes $nod --acc2tax $acc --metaDMG-cpp /data2/programs/metaDMG-cpp/metaDMG-cpp --config-file config.metaDMG.$run.yaml --bayesian --parallel-samples 1 --cores-per-sample $threads --custom-database
metaDMG compute config.metaDMG.$run.yaml
metaDMG convert --output $run.csv config.metaDMG.$run.yaml
metaDMG plot config.metaDMG.$run.yaml --output $run.plots.pdf
`

[miwipe]

Hi Libby,

The metaDMG-core has been discontinued and we are no longer supporting it. The base code with newest updates can however be found here https://github.com/metaDMG-dev/metaDMG-cpp. In essence this is almost the same, some slight modifications and improvements have been made.

While it basically contains the same outputs, it no longer has the GUI.

For tracking back your lost reads, I would go back to your bam file and check whether they mapped, to what references they mapped, and how good a fit they are to the reference. I do not know what similarity thresholds you used but these can be the reason why you don’t classify the reads. There is not an option to print intermediate files, for what has been used to calculate the LCA. but simply you sam file is this.

Then next you should double check that your taxonomy is in order, that the accessions to taxid are ok, and that you taxid to the taxonomy (names and nodes dump files are ok). This can also be the reason.

Again, we no longer use this version of metaDMG. For the use of metaDMG-cpp you can also look at our yearly PhD course day2 walk-through https://github.com/GeoGenetics/aeGenomicsCourse2024/wiki/Day-2

I hope this helps.
Mikkel