contDeam.pl failing in detecting deamination

[NicRamb]

Hi! I'm experiencing some problems when running contDeam.pl on some BAM files. It prints the nan error and apparently fails to detect deamination (since the average depth of coverage is not that low, between 2X and 12X. If I try running the single subparts, it only works by omitting the -endo option in bam2prof, then the other parts (contDeam, endoCalle, and mtCont) work and the final results seem reasonable. I don't know whether this solution is correct or I'm just adding more errors. Do you have any ideas about the cause of this issue?

[grenaud]

Ciao Nicola! Thank you for your interest! Could you post the output of bam2prof? I suspect that this has to do with the lack of fragments to get a sufficiently strong statistical signal. The -endo option conditions on one end being deaminated to measure the other one. This fails if you don't have enough damage or aDNA fragments.

[NicRamb]

Hi! Thank you for your answer! Here it is the contDeam.pl command and error for a sample with an average depth of coverage of ~12X:

contDeam.pl --library double --lengthDeam 5 --out testsample testsample.bam

running cmd schmutzi/src/bam2prof -length 5 -endo -double -5p testsample.endo.5p.prof -3p testsample.endo.3p.prof testsample.bam running cmd schmutzi/src/contDeam -deamread -deam5p testsample.endo.5p.prof -deam3p testsample.endo.3p.prof -log testsample.cont.deam testsample.bam Utils.cpp: destringify() Unable to convert string="-nan" system cmd schmutzi/src/contDeam -deamread -deam5p testsample.endo.5p.prof -deam3p testsample.endo.3p.prof -log testsample.cont.deam testsample.bam failed: 256 at schmutzi/src/contDeam.pl line 22.

While here it is the bam2prof command for the same sample with the -endo option:

schmutzi/src/bam2prof -length 5 -endo -double -5p testsample_withEndo.endo.5p.prof -3p testsample_withEndo.endo.3p.prof testsample.bam

and here without the -endo option:

schmutzi/src/bam2prof -length 5 -double -5p testsample_NOendo.endo.5p.prof -3p testsample_NOendo.endo.3p.prof testsample.bam

Here both the outputs from bam2prof:

bam2prof_output.zip

[grenaud]

Ciao! I see very little damage on the 5' end, it this normal? Which protocol is this? Also, why does G->A go up after the 2nd base at the 3' end?

[NicRamb]

The protocol is as in Lindo et al., 2017 and Scheib et al., 2018 (Illinois University sections). Damages on 5' and 3' ends are as in figure S1 in Lindo et al., 2017 (443 and 302 samples). The NEBNext Ultra DNA library prep kit was used with blunt-end repair. Honestly, I'm not sure why the damage is so low on 5' end, could it be due to single-end sequencing? Thanks again for your answers!

[grenaud]

Ok I would try to determine why the damage patterns at the 5' end are so low before moving forward, perhaps it's a bioinformatics processing problem. You mean single-end sequencing or single-strand library prep?

[NicRamb]

I mean single-end sequencing. I will go through the bioinformatic process and try to figure this out, thanks again for your answers and your time!

[grenaud]

I can tell you how to get around that problem but I think it's a good first step to check if everything worked well on your end.

[NicRamb]

Hi, sorry for the late in the answer. I went through the bioinformatic process and it seems ok. I think the low damage could be due to the blunt/TA protocol with the use of USER enzyme. Thanks again for your answers and advice.

[grenaud]

ok, if I were you, I would be sure that a non-USER treated library does have some damage. If so, I would run: Software/schmutzi/src/bam2makeSchmutziHaplogrep2.pl with the option --nodeam and put some prior for contamination like 0.2

[ywakiyama]

Hello,

I have a question related to this Issue. Is it possible to use schmutzi for analysis of mtDNA reads processed by NEBNext Ultra DNA library prep kit (which Nicola said above and I am using)? In the library preparation using NEBNext Ultra DNA library prep kit, USER is used in order to cut adapter (please see this picture https://www.nebj.jp/jp/info/E7370_Illumina_UltraDNAWorkflow_0114.png) . If a fragment has Uracil at the 5' end, USER also remove it, 5' end's adapter is off and the library is not duplicated. As a result, the damage patterns at the 5' end are not detected.

[grenaud]

Hello! Is there any damage to speak of? What people in Copenhagen usually do is have a USER treated library and a non-USER to make sure it is ancient DNA. Also, what is the mitochondrial coverage? Did you align to the mt in isolation.

[ywakiyama]

Thank you for your prompt reply! In my output of mapDamage, G->A at the 3' end is seen, but C->T at the 5' end hardly seen (The plots show C->T are almost parallel to the horizontal axis) in all samples. I aligned the reads to rCRS using BWA and the average depth of mitochondrial coverage is more than 24X.

[grenaud]

Ok, so you will not be able to use contDeam.pl. It needs some damage to drive the signal.

However, you can use the rest of the software (schmutzi.pl). I would use:

schmutzi.pl --contprior 0.3 --notusepredC

There is a risk that if there is no damage, you do not have enough statistical signal to call the endogenous space. endoCaller essentially relies on this damage to call the endogenous base.

[ywakiyama]

I am sorry to know that I will not be able to use contDeam.pl. So, I will try schimtuzi.pl. Thank you for your kind instructon!

[grenaud]

Well it is not a problem but the absence of damage does not allow us to say that this is bona fide ancient DNA.

On Fri, Mar 12, 2021 at 3:52 PM ywakiyama @.***> wrote:

I am sorry to know that I will not be able to use contDeam.pl. So, I will try schimtuzi.pl. Thank you for your kind instructon!

— You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub https://github.com/grenaud/schmutzi/issues/10#issuecomment-797538314, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAQRNI3BWTEHWETWCKCMHTTTDITCBANCNFSM4ROLY42A .