deamSim terminated successfuly but 0 sequences were written

[hyl317]

Hi,

I am trying to use the deamSim sub-program to simulate some DNA post-mortem damage by giving a BAM file as the input. The command I used is

./deamSim -mapdamage ../examplesMapDamage/results_Ust_Ishim/misincorporation.txt single -b [path to my BAM file]

However, the program terminates immediately with the following message:
Program ./deamSim terminated succesfully, wrote 0 sequences

This is apparently not a successful termination. I am wondering what have gone wrong. My BAM file is a modern BAM file downloaded SGDP. Hope you can help me figure this out! Thanks in advance.

[grenaud]

you need to put a bam file, you put
[path to my BAM file]

[hyl317]

Hi Gabriel,

Thanks for the prompt reply! Do you mean I need to supply a input stream instead of a bam file name?

I tried three commands,
samtools view -u [path to my BAM] | ./deamSim -mapdamage ../examplesMapDamage/results_Ust_Ishim/misincorporation.txt single -b

and

./deamSim -mapdamage ../examplesMapDamage/results_Ust_Ishim/misincorporation.txt single -b <(samtools view /mnt/archgen/users/yilei/Data/hapCon_sim/SGDPcontaminant/B_French-3.chrX.bam)

./deamSim -mapdamage ../examplesMapDamage/results_Ust_Ishim/misincorporation.txt single -b <(samtools view -u [path to my BAM])

None of the above works. The first gives "Unable to open the file -b", the second "Could not open input BAM file /dev/fd/63" and the third "Program ./deamSim terminated succesfully, wrote 0 sequences".

I am wondering what is the correct way to supply a BAM file?

[grenaud]

./deamSim -mapdamage ../examplesMapDamage/results_Ust_Ishim/misincorporation.txt single -b <(samtools view /mnt/archgen/users/yilei/Data/hapCon_sim/SGDPcontaminant/B_French-3.chrX.bam)

Why not just write;
./deamSim -mapdamage ../examplesMapDamage/results_Ust_Ishim/misincorporation.txt single -b /mnt/archgen/users/yilei/Data/hapCon_sim/SGDPcontaminant/B_French-3.chrX.bam

[hyl317]

In my very first post, when I say [path to my BAM file], I meant "/mnt/archgen/users/yilei/Data/hapCon_sim/SGDPcontaminant/B_French-3.chrX.bam", not the literal [path to my BAM file]. Sorry for the confusion. Still, by providing the path to BAM file as you have suggested, I got
Program ./deamSim terminated succesfully, wrote 0 sequences

[grenaud]

Ok just to be sure, these are not aligned fragments correct? deamSim is for unaligned sequences that need to be realigned.

[hyl317]

Hi Gabriel,

thanks for the quick reply. Yes my BAM file are aligned reads... Our goal is to throw down some aDNA-charactersitic damage to modern BAM files and investigate how different deamination levels affect certain downsteram analysis. What's your suggestion then? Do you recommend, first convert the aligned BAM file to fastq (samtools bam2fq), and then use deamSim, and then realign to the reference genome? Thanks in advance.

[grenaud]

Two ways to go about this:

1. create an unaligned bam file stripped of all mapping information
2. transferring them to fastq and fastq to bam. You can even do this on 1 line using samtools.

[hyl317]

sorry for multiple questions. I tried to transofrming my bam to fastq like the follwoing,

samtools bam2fq my.bam > my.fq
then I ran
./deamSim -mapdamage ../examplesMapDamage/results_Ust_Ishim/misincorporation.txt single my.fq
However, I got
Parse error#1 with line @HS2000-1259_195:2:2308:3592:28380/1

[grenaud]

reread my point 2, you missed converting back to bam.

[hyl317]

sorry for my mistake. Now I tried,samtools import -o test.bam my.fq
and then
./deamSim -mapdamage ../examplesMapDamage/results_Ust_Ishim/misincorporation.txt single test.bam
however, deamSim complains about the bam header,
Parse error#1 with line BAM5@CO Reverse with: samtools fastq -N -o paired.fastq

[grenaud]

samtools import does not really exist anymore, can you update samtools?

[hyl317]

I think samtools import is actually a rather new addition to samtools. see https://github.com/samtools/samtools/releases/tag/1.13
bullet point 4, "A new tool, samtools import, has been added. It reads one or more FASTQ files and converts them into unmapped SAM, BAM or CRAM."

But anyways, have you tried other samtools command for the purpose of converting fastq to unmapped BAM that works for deamSim?

[grenaud]

I see! I got:
samtools import
[main] "samtools import" has been removed. Please use "samtools view" instead.

This is weird, can you email me the first 100 records of your bam file?

[hyl317]

Hi Gabriel, See the attached file for the input (aligned) BAM I started with. Thanks again for your prompt reply! best, Yilei

…

On Fri, Jan 28, 2022 at 1:27 PM Gabriel Renaud ***@***.***> wrote: I see! I got: samtools import [main] "samtools import" has been removed. Please use "samtools view" instead. This is weird, can you email me the first 100 records of your bam file? — Reply to this email directly, view it on GitHub <#15 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ALNWTMG3SZG77D5SC3MES6TUYKDTNANCNFSM5M6VSM4A> . Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>. You are receiving this because you authored the thread.Message ID: ***@***.***>

[grenaud]

please send it via email, not github.

[grenaud]

I see now these are paired-end sequences, not single end. So I did the following, I just considered the forward reads. otherwise, the assumption of independence between fragments is gone.

samtools fastq -f 0x40 test.bam |samtools import -O BAM /dev/stdin > in.bam
[M::bam2fq_mainloop] discarded 0 singletons
[M::bam2fq_mainloop] processed 435 reads
gargammel/src/deamSim -b out.bam -matfile gargammel/src/matrices/double- in.bam

Program /home/gabriel/Software/gargammel/src/deamSim terminated succesfully, wrote 435 sequences