guoweilong / cgmaptools

toolbox for analysing BS-seq data, advance features in SNV, ASM and DMR
https://cgmaptools.github.io
61 stars 26 forks source link

cgmaptools convert bam2cgmap exiting without error #70

Open lukesarre opened 2 weeks ago

lukesarre commented 2 weeks ago

Hello, thank you for this wonderful tool.

I'm running into a problem where cgmaptools is exiting early, without processing the aligned reads into the CGmap.gz file

This is happening only with some inputs, and not with others. For example, this one is working fine:

cgmaptools convert bam2cgmap --bam SRR042643.sorted.dedup.bam --genome ./genomes/Phycomyces_blakesleeanus_v2_total.fna

As an example of the inputs:

[btx832@frontend5 fungalInvestigations]$ samtools view SRR042643.sorted.dedup.bam | head -n 5
SRR042643.2218034       16      Pblake071009.mitochondria.fasta 1       255     5M1I59M *       0      0
        ATTACCTCAAATAAACTCAATATTATTATCAAAATAACTACAAATACTTATTTAATCTTAAAATA       *       XO:Z:-FW       XS:i:0   NM:i:4  XM:Z:---z-z----------------yx--------zyx-----z--z-y----zzz-zy---------  XG:Z:TA_TATCTCAAGATTAAATAAGTATCCGTAGTTATCCCGATAACAACACTGAGCCCACCTGAGGCAA_NN
SRR042643.1247985       0       Pblake071009.mitochondria.fasta 2       255     1M3I72M *       0      0
        TGTTGTTTTAGGTGGGTTTAGTGTTGTTATTGGGATAATTATGGATATTTATTTAATTTTGAGATATATATGAAAT    *       XO:Z:+FW
        XS:i:0  NM:i:3  XM:Z:-----zz-y-------z-y-----------x-------z--x-----z---------z----------z-------       XG:Z:NT_TGCCTCAGGTGGGCTCAGTGTTGTTATCGGGATAACTACGGATACTTATTTAATCTTGAGATATACATGAAAT_GG
SRR042643.8908027       16      Pblake071009.mitochondria.fasta 2       255     64M     *       0      0
        TACCTCAAATAAACTCAATATTANTATCAAAATAACTACAAATACTTATTTAATCTTAAAATAT        *       XO:Z:-FW       XS:i:0   NM:i:0  XM:Z:----z-z----------------yx--------zyx-----z--z-y----zzz-zy-----z-   XG:Z:GT_ATATCTCAAGATTAAATAAGTATCCGTAGTTATCCCGATAACAACACTGAGCCCACCTGAGGCA_AN
SRR042643.9549560       16      Pblake071009.mitochondria.fasta 2       255     76M     *       0      0
        TACCTCAAATAAACTCAATATTATTATCAAAATAACTACAAATACTTATTTAATCTTAAAATATACATAAAATAAA    *       XO:Z:-FW
        XS:i:0  NM:i:0  XM:Z:zzz----z--------z-z----------------yx--------zyx-----z--z-y----zzz-zy-----z-       XG:Z:TC_CCCATTTCATGTATATCTCAAGATTAAATAAGTATCCGTAGTTATCCCGATAACAACACTGAGCCCACCTGAGGCA_AN
SRR042643.10640227      0       Pblake071009.mitochondria.fasta 3       255     4M3I69M *       0      0
        GTTGTTTTAGGTGGGTTTAGTGTTGTTATTGGGATAATTATGGATATTTATTTAATTTTGAGATATATATGAAATG    *       XO:Z:+FW
        XS:i:0  NM:i:4  XM:Z:-zz----y-------z-y-----------x-------z--x-----z---------z----------z--------       XG:Z:TT_GCCTCAGGTGGGCTCAGTGTTGTTATCGGGATAACTACGGATACTTATTTAATCTTGAGATATACATGAAATG_GG
[btx832@frontend5 fungalInvestigations]$ head -n 5 ./genomes/Phycomyces_blakesleeanus_v2_total.fna
>Pblake071009.mitochondria.fasta
TTGCCTCAGGTGGGCTCAGTGTTGTTATCGGGATAACTACGGATACTTAT
TTAATCTTGAGATATACATGAAATGGGGATATCACAAACTTGTATCTACT
ATTTATCCAGGAAGAAACCCTCCTGGATAATGTAATGAGGAGATATCTCG
GACCCCCGGAGGAATCATCGGACTCTGGATCCGATTCCGGGGACGAAGAT

Whereas this one is not working (there are reads aligned to many more contigs):

[btx832@frontend5 fungalInvestigations]$ cgmaptools convert bam2cgmap --bam Lacbi2.sorted.dedup.bam --genome ./genomes/Lacbi2_AssemblyScaffolds.fna
# input bam file: Lacbi2.sorted.dedup.bam
# source genome file: ./genomes/Lacbi2_AssemblyScaffolds.fna
# Processing reads from: LG_1
[btx832@frontend5 fungalInvestigations]$

As an example of the inputs:

[btx832@frontend5 fungalInvestigations]$ samtools view Lacbi2.sorted.dedup.bam | head -n 5
SRR042633.9506604       16      LG_1    11      255     65M     *       0       0       AACCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCTAACCCTA       *       XO:Z:-FW        XS:i:0  NM:i:0  XM:Z:-----------------------------------------------------------------  XG:Z:GT_TAGGGTTAGGTTAGGGTTAGGGTTAGGGTTAGGGTTAGGGTTAGGGTTAGGGTTAGGGTTAGGTT_AG
SRR042633.8978555       16      LG_1    11      255     9M1D42M1I5M1I18M        *       0       0      AACCTAACCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCCTAACCCTAACCCTAACCCTAACCT     *       XO:Z:-FW       XS:i:0   NM:i:3  XM:Z:-----------------------------------------------------------------------------     XG:Z:TT_AGGTTAGGGTTAGGGTTAGGTTAGGGTTAGGGTTAGGGTTAGGGTTAGGGTTAGGGTTAGGGTTAGGGTTAGGTT_AG
SRR042633.4951499       16      LG_1    25      255     54M1D5M1I16M    *       0       0       CCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCTAACCCTAACCTAACCCTAACCCTAACCCTAAC    *       XO:Z:-FW        XS:i:0 NM:i:2   XM:Z:-----------------------------------------------------------------------------      XG:Z:GG_GTTAGGGTTAGGGTTAGGTTAGGGTTAGGGTTAGGTTAGGGTTAGGGTTAGGGTTAGGGTTAGGGTTAGGGTTAGG_GT
SRR042633.3437571       16      LG_1    52      255     10M1D5M1I35M1D6M1I10M   *       0       0      AACCCTAACCTAACCCTAACCCTAACCCTAACCTAACCCTAACCCTAACCCTAACCTTAACCCTAACC     *       XO:Z:-FW        XS:i:0 NM:i:4   XM:Z:----------------------------------------------------------------------     XG:Z:TA_GGTTAGGGTTAGGTTAGGGGTTAGGGTTAGGGTTAGGTTAGGGTTAGGGTTAGGTTAGGGTTAGGGTT_AG
SRR042633.4618143       16      LG_1    356     255     41M     *       0       0       TACTCTTACAATTCAATCAAAAACTAAACTTATCATATAAA       *       XO:Z:-FW        XS:i:0  NM:i:1  XM:Z:---------z---z-y--z-z-----z------z-----z-  XG:Z:TT_TTTATATGACAAGCTCAGCTCTTGATCAAATTGCAAGAGCA_TC
[btx832@frontend5 fungalInvestigations]$ head -n 5 ./genomes/Lacbi2_AssemblyScaffolds.fna
>LG_1
CCCTAACCCTAACCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCTAA
CCCTAACCCTAACCTAACCCTAACCCTAACCCCTAACCTAACCCTAACCTAACCCTAACCCAAAACCTAA
CCTAACCCTAACCAACAAGACACCCGAACCATTAACCGAGACACCGAGACATACCCGTTTCCATCCTAAC
CCCTTAACCGAGACACCCGAACCCTTAACCGTGTTTCCAAAGGACATTCCAAATCCATGGCACAATAACC

I can't tell the difference, as they've both been processed in the same way (aligned with bs_seeker2). Have you experienced this bug before?

guoweilong commented 2 weeks ago

As "cgmaptools convert" may not work well if there are too many contigs, you may try to split bam files by chromosomes, and then run the command for each chr individually.

Weilong

At 2024-09-23 21:59:23, "lukesarre" @.***> wrote:

Hello, thank you for this wonderful tool.

I'm running into a problem where cgmaptools is exiting early, without processing the aligned reads into the CGmap.gz file

This is happening only with some inputs, and not with others. For example, this one is working fine:

cgmaptools convert bam2cgmap --bam SRR042643.sorted.dedup.bam --genome ./genomes/Phycomyces_blakesleeanus_v2_total.fna

As an example of the inputs:

@.** fungalInvestigations]$ samtools view SRR042643.sorted.dedup.bam | head -n 5 SRR042643.2218034 16 Pblake071009.mitochondria.fasta 1 255 5M1I59M 0 0 ATTACCTCAAATAAACTCAATATTATTATCAAAATAACTACAAATACTTATTTAATCTTAAAATA XO:Z:-FW XS:i:0 NM:i:4 XM:Z:---z-z----------------yx--------zyx-----z--z-y----zzz-zy--------- XG:Z:TA_TATCTCAAGATTAAATAAGTATCCGTAGTTATCCCGATAACAACACTGAGCCCACCTGAGGCAA_NN SRR042643.1247985 0 Pblake071009.mitochondria.fasta 2 255 1M3I72M 0 0 TGTTGTTTTAGGTGGGTTTAGTGTTGTTATTGGGATAATTATGGATATTTATTTAATTTTGAGATATATATGAAAT XO:Z:+FW XS:i:0 NM:i:3 XM:Z:-----zz-y-------z-y-----------x-------z--x-----z---------z----------z------- XG:Z:NT_TGCCTCAGGTGGGCTCAGTGTTGTTATCGGGATAACTACGGATACTTATTTAATCTTGAGATATACATGAAAT_GG SRR042643.8908027 16 Pblake071009.mitochondria.fasta 2 255 64M 0 0 TACCTCAAATAAACTCAATATTANTATCAAAATAACTACAAATACTTATTTAATCTTAAAATAT XO:Z:-FW XS:i:0 NM:i:0 XM:Z:----z-z----------------yx--------zyx-----z--z-y----zzz-zy-----z- XG:Z:GT_ATATCTCAAGATTAAATAAGTATCCGTAGTTATCCCGATAACAACACTGAGCCCACCTGAGGCA_AN SRR042643.9549560 16 Pblake071009.mitochondria.fasta 2 255 76M 0 0 TACCTCAAATAAACTCAATATTATTATCAAAATAACTACAAATACTTATTTAATCTTAAAATATACATAAAATAAA XO:Z:-FW XS:i:0 NM:i:0 XM:Z:zzz----z--------z-z----------------yx--------zyx-----z--z-y----zzz-zy-----z- XG:Z:TC_CCCATTTCATGTATATCTCAAGATTAAATAAGTATCCGTAGTTATCCCGATAACAACACTGAGCCCACCTGAGGCA_AN SRR042643.10640227 0 Pblake071009.mitochondria.fasta 3 255 4M3I69M 0 0 GTTGTTTTAGGTGGGTTTAGTGTTGTTATTGGGATAATTATGGATATTTATTTAATTTTGAGATATATATGAAATG * XO:Z:+FW XS:i:0 NM:i:4 XM:Z:-zz----y-------z-y-----------x-------z--x-----z---------z----------z-------- XG:Z:TT_GCCTCAGGTGGGCTCAGTGTTGTTATCGGGATAACTACGGATACTTATTTAATCTTGAGATATACATGAAATG_GG

@.*** fungalInvestigations]$ head -n 5 ./genomes/Phycomyces_blakesleeanus_v2_total.fna

Pblake071009.mitochondria.fasta TTGCCTCAGGTGGGCTCAGTGTTGTTATCGGGATAACTACGGATACTTAT TTAATCTTGAGATATACATGAAATGGGGATATCACAAACTTGTATCTACT ATTTATCCAGGAAGAAACCCTCCTGGATAATGTAATGAGGAGATATCTCG GACCCCCGGAGGAATCATCGGACTCTGGATCCGATTCCGGGGACGAAGAT

Whereas this one is not working (there are reads aligned to many more contigs):

@.*** fungalInvestigations]$ cgmaptools convert bam2cgmap --bam Lacbi2.sorted.dedup.bam --genome ./genomes/Lacbi2_AssemblyScaffolds.fna

input bam file: Lacbi2.sorted.dedup.bam

source genome file: ./genomes/Lacbi2_AssemblyScaffolds.fna

Processing reads from: LG_1

@.*** fungalInvestigations]$

As an example of the inputs:

@.** fungalInvestigations]$ samtools view Lacbi2.sorted.dedup.bam | head -n 5 SRR042633.9506604 16 LG_1 11 255 65M 0 0 AACCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCTAACCCTA XO:Z:-FW XS:i:0 NM:i:0 XM:Z:----------------------------------------------------------------- XG:Z:GT_TAGGGTTAGGTTAGGGTTAGGGTTAGGGTTAGGGTTAGGGTTAGGGTTAGGGTTAGGGTTAGGTT_AG SRR042633.8978555 16 LG_1 11 255 9M1D42M1I5M1I18M 0 0 AACCTAACCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCCTAACCCTAACCCTAACCCTAACCT XO:Z:-FW XS:i:0 NM:i:3 XM:Z:----------------------------------------------------------------------------- XG:Z:TT_AGGTTAGGGTTAGGGTTAGGTTAGGGTTAGGGTTAGGGTTAGGGTTAGGGTTAGGGTTAGGGTTAGGGTTAGGTT_AG SRR042633.4951499 16 LG_1 25 255 54M1D5M1I16M 0 0 CCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCTAACCCTAACCTAACCCTAACCCTAACCCTAAC XO:Z:-FW XS:i:0 NM:i:2 XM:Z:----------------------------------------------------------------------------- XG:Z:GG_GTTAGGGTTAGGGTTAGGTTAGGGTTAGGGTTAGGTTAGGGTTAGGGTTAGGGTTAGGGTTAGGGTTAGGGTTAGG_GT SRR042633.3437571 16 LG_1 52 255 10M1D5M1I35M1D6M1I10M 0 0 AACCCTAACCTAACCCTAACCCTAACCCTAACCTAACCCTAACCCTAACCCTAACCTTAACCCTAACC XO:Z:-FW XS:i:0 NM:i:4 XM:Z:---------------------------------------------------------------------- XG:Z:TA_GGTTAGGGTTAGGTTAGGGGTTAGGGTTAGGGTTAGGTTAGGGTTAGGGTTAGGTTAGGGTTAGGGTT_AG SRR042633.4618143 16 LG_1 356 255 41M 0 0 TACTCTTACAATTCAATCAAAAACTAAACTTATCATATAAA * XO:Z:-FW XS:i:0 NM:i:1 XM:Z:---------z---z-y--z-z-----z------z-----z- XG:Z:TT_TTTATATGACAAGCTCAGCTCTTGATCAAATTGCAAGAGCA_TC

@.*** fungalInvestigations]$ head -n 5 ./genomes/Lacbi2_AssemblyScaffolds.fna

LG_1 CCCTAACCCTAACCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCTAA CCCTAACCCTAACCTAACCCTAACCCTAACCCCTAACCTAACCCTAACCTAACCCTAACCCAAAACCTAA CCTAACCCTAACCAACAAGACACCCGAACCATTAACCGAGACACCGAGACATACCCGTTTCCATCCTAAC CCCTTAACCGAGACACCCGAACCCTTAACCGTGTTTCCAAAGGACATTCCAAATCCATGGCACAATAACC

I can't tell the difference, as they've both been processed in the same way (aligned with bs_seeker2). Have you experienced this bug before?

— Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you are subscribed to this thread.Message ID: @.***>