vgteam / vg

tools for working with genome variation graphs
https://biostars.org/tag/vg/
Other
1.1k stars 194 forks source link

BAM file surjected from GAM includes unmapped reads #2581

Closed kubranarci closed 4 years ago

kubranarci commented 4 years ago

I want to use VG map alignment file to use with GATK4 haplotypecaller. In order to do that I used the instructions mentioned here (https://github.com/vgteam/vg/wiki/Working-with-a-whole-genome-variation-graph) to work through whole-genome. Then I converted (sorted) GAM file to BAM file using VG surject (to bam option). Yet, in sorted and indexed BAM file, ALL reads seems to be unmapped (mostly flag 4). Do you have any idea on what I am missing though the process?

Mostly like this: HISEQ1:22:H9UJNADXX:2:1106:20572:14228/1 4 0 0 0 0 AGTTGGCAGAGAAGGCCTCTTGTGCTCACGCTCCCCCGCAGTCCCCAGCCCTTCTGCCTTTCTCCCCCGACACTGCTGCACCAGAGTGAAAGGGCTATGGCAAGGGGGTGTCATCTGAGGAGTATTAAGAATGCAGATTCCTGGGCCT @@@FFFFFHHHHGJIIJIJFIJEHIGIJIJJJJIJJIHIJIIJJIJIEIJH@HGFFFFEFEEEDEDDDBDDBBCDCDDDDDDDDDDCCCACDDDD?BDDDBDDBBBDB5>BCCDDDC>CBA@9?CDDEDDDDDDACDCCDDDDDDBB@ HISEQ1:22:H9UJNADXX:2:2212:3639:78995/1 4 0 0 0 0 TTCCCAGCTCACCATCAGATCTATTCTGAGCTTGTGCATTTGGCCTAGGCAAGACTGGCTGGCTGCTTGGGGCAGAGGGAGGGGCCCTGGGGACAAAAGGAGCTGGGTACAGGGGTTGCTCTAGCAGGCACAGCACAGACCCCGGAGA CCCFFFFFHHHHHJJJJJJJJJJJJJJGJJJJJJIJIJJJJJJJJJIJJJJJJJJJJJJJJJJJJJJJJJHHFFFDDDDDDDDDDDDDDDDDDDDDDDDDDBBDDDD?CDDEDDDBBDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDB HISEQ1:22:H9UJNADXX:2:2212:3639:78995/1 4 0 0 0 0 TTCCCAGCTCACCATCAGATCTATTCTGAGCTTGTGCATTTGGCCTAGGCAAGACTGGCTGGCTGCTTGGGGCAGAGGGAGGGGCCCTGGGGACAAAAGGAGCTGGGTACAGGGGTTGCTCTAGCAGGCACAGCACAGACCCCGGAGA CCCFFFFFHHHHHJJJJJJJJJJJJJJGJJJJJJIJIJJJJJJJJJIJJJJJJJJJJJJJJJJJJJJJJJHHFFFDDDDDDDDDDDDDDDDDDDDDDDDDDBBDDDD?CDDEDDDBBDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDB HISEQ1:22:H9UJNADXX:1:2204:1672:85287/1 4 0 0 0 0 CAGATCTATTCTGAGCTTGTGCATTTGGCCTAGGCAAGACTGGCTGGCTGCTTGGGGCAGAGGGAGGGGCACTGGGGACAAAAGGAGCTGGGTACAGGGATTGCTCTAGCAGGGACAGCACAGACCCCGGAGAAAGCCCACAACAATA CCCFFFFFHHHHHJJJJJGGHIJJJJJJJJJJJJJJJJJJJJJIIJJJJJJJJJJJJJJJIJJHFHFDDD';?CDDDBDDDDDDDDDDDDDDCCDDEDD(+3<CDDDDDDDDB)2<BDDDDDDDDDDDDBDD@BDDDDDDDDDDDDDC HISEQ1:22:H9UJNADXX:1:2204:1672:85287/1 4 0 0 * 0 0 CAGATCTATTCTGAGCTTGTGCATTTGGCCTAGGCAAGACTGGCTGGCTGCTTGGGGCAGAGGGAGGGGCACTGGGGACAAAAGGAGCTGGGTACAGGGATTGCTCTAGCAGGGACAGCACAGACCCCGGAGAAAGCCCACAACAATA CCCFFFFFHHHHHJJJJJGGHIJJJJJJJJJJJJJJJJJJJJJIIJJJJJJJJJJJJJJJIJJHFHFDDD';?CDDDBDDDDDDDDDDDDDDCCDDEDD(+3<CDDDDDDDDB)2<BDDDDDDDDDDDDBDD@BDDDDDDDDDDDDDC

ekg commented 4 years ago

Note that vg map will write BAM directly. Try running a few reads that way and see if they are still unmapped.

On Fri, Dec 20, 2019, 09:41 Kübra NARCI notifications@github.com wrote:

I want to use VG map alignment file to use with GATK4 haplotypecaller. In order to do that I used the instructions mentioned here ( https://github.com/vgteam/vg/wiki/Working-with-a-whole-genome-variation-graph) to work through whole-genome. Then I converted (sorted) GAM file to BAM file using VG surject (to bam option). Yet, in sorted and indexed BAM file, ALL reads seems to be unmapped (mostly flag 4). Do you have any idea on what I am missing though the process?

Mostly like this: HISEQ1:22:H9UJNADXX:2:1106:20572:14228/1 4 0 0 0 0 AGTTGGCAGAGAAGGCCTCTTGTGCTCACGCTCCCCCGCAGTCCCCAGCCCTTCTGCCTTTCTCCCCCGACACTGCTGCACCAGAGTGAAAGGGCTATGGCAAGGGGGTGTCATCTGAGGAGTATTAAGAATGCAGATTCCTGGGCCT @@@FFFFFHHHHGJIIJIJFIJEHIGIJIJJJJIJJIHIJIIJJIJIEIJH@HGFFFFEFEEEDEDDDBDDBBCDCDDDDDDDDDDCCCACDDDD ?BDDDBDDBBBDB5>BCCDDDC>CBA@9?CDDEDDDDDDACDCCDDDDDDBB@ HISEQ1:22:H9UJNADXX:2:2212:3639:78995/1 4 0 0 0 0 TTCCCAGCTCACCATCAGATCTATTCTGAGCTTGTGCATTTGGCCTAGGCAAGACTGGCTGGCTGCTTGGGGCAGAGGGAGGGGCCCTGGGGACAAAAGGAGCTGGGTACAGGGGTTGCTCTAGCAGGCACAGCACAGACCCCGGAGA CCCFFFFFHHHHHJJJJJJJJJJJJJJGJJJJJJIJIJJJJJJJJJIJJJJJJJJJJJJJJJJJJJJJJJHHFFFDDDDDDDDDDDDDDDDDDDDDDDDDDBBDDDD?CDDEDDDBBDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDB HISEQ1:22:H9UJNADXX:2:2212:3639:78995/1 4 0 0 0 0 TTCCCAGCTCACCATCAGATCTATTCTGAGCTTGTGCATTTGGCCTAGGCAAGACTGGCTGGCTGCTTGGGGCAGAGGGAGGGGCCCTGGGGACAAAAGGAGCTGGGTACAGGGGTTGCTCTAGCAGGCACAGCACAGACCCCGGAGA CCCFFFFFHHHHHJJJJJJJJJJJJJJGJJJJJJIJIJJJJJJJJJIJJJJJJJJJJJJJJJJJJJJJJJHHFFFDDDDDDDDDDDDDDDDDDDDDDDDDDBBDDDD?CDDEDDDBBDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDB HISEQ1:22:H9UJNADXX:1:2204:1672:85287/1 4 0 0 0 0 CAGATCTATTCTGAGCTTGTGCATTTGGCCTAGGCAAGACTGGCTGGCTGCTTGGGGCAGAGGGAGGGGCACTGGGGACAAAAGGAGCTGGGTACAGGGATTGCTCTAGCAGGGACAGCACAGACCCCGGAGAAAGCCCACAACAATA CCCFFFFFHHHHHJJJJJGGHIJJJJJJJJJJJJJJJJJJJJJIIJJJJJJJJJJJJJJJIJJHFHFDDD';?CDDDBDDDDDDDDDDDDDDCCDDEDD(+3<CDDDDDDDDB)2<BDDDDDDDDDDDDBDD@BDDDDDDDDDDDDDC HISEQ1:22:H9UJNADXX:1:2204:1672:85287/1 4 0 0 * 0 0 CAGATCTATTCTGAGCTTGTGCATTTGGCCTAGGCAAGACTGGCTGGCTGCTTGGGGCAGAGGGAGGGGCACTGGGGACAAAAGGAGCTGGGTACAGGGATTGCTCTAGCAGGGACAGCACAGACCCCGGAGAAAGCCCACAACAATA CCCFFFFFHHHHHJJJJJGGHIJJJJJJJJJJJJJJJJJJJJJIIJJJJJJJJJJJJJJJIJJHFHFDDD';?CDDDBDDDDDDDDDDDDDDCCDDEDD(+3<CDDDDDDDDB)2<BDDDDDDDDDDDDBDD@BDDDDDDDDDDDDDC

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/vgteam/vg/issues/2581?email_source=notifications&email_token=AABDQEKWLCPRGSMG4B2EVNTQZSAKZA5CNFSM4J52KPTKYY3PNVWWK3TUL52HS4DFUVEXG43VMWVGG33NNVSW45C7NFSM4IB43ZBA, or unsubscribe https://github.com/notifications/unsubscribe-auth/AABDQEMM4W3B5PZH3MUZDRLQZSAKZANCNFSM4J52KPTA .

kubranarci commented 4 years ago

Hi, sorry for not mentioning, I done that too but it was not working also. Also, I just inspected on GAM file though VG stats and the alignments seem to be aligned in there;

The result: Total alignments: 10848688 Total primary: 10848688 Total secondary: 0 Total aligned: 10646043 Total perfect: 8105365 Total gapless (softclips allowed): 10300682 Insertions: 818749 bp in 242882 read events

ekg commented 4 years ago

Would you please provide a reproducible example of reads that should be mapped but don't map with vg map --surject-to bam?

On Fri, Dec 20, 2019 at 10:02 AM Kübra NARCI notifications@github.com wrote:

Hi, sorry for not mentioning, I done that too but it was not working also. Also, I just inspected on GAM file though VG stats and the alignments seem to be aligned in there;

The result: Total alignments: 10848688 Total primary: 10848688 Total secondary: 0 Total aligned: 10646043 Total perfect: 8105365 Total gapless (softclips allowed): 10300682 Insertions: 818749 bp in 242882 read events

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/vgteam/vg/issues/2581?email_source=notifications&email_token=AABDQEPMOPZ2EOQVM7IWL6LQZSC3RA5CNFSM4J52KPTKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEHMK47I#issuecomment-567848573, or unsubscribe https://github.com/notifications/unsubscribe-auth/AABDQEJNS2UUHS4LDEWOMY3QZSC3RANCNFSM4J52KPTA .

kubranarci commented 4 years ago

The number flagged as 4 is 21697376: samtools view ~/Desktop/softwares/vg/test.sam -f 4 |wc -l

ekg commented 4 years ago

@kubranarci we have continuous integration testing of vg that should catch a total failure of alignment. If you've run into something we're not testing, we need an example of data and how to build a graph that reproduces the problem you're seeing in order to resolve it.

Please start by sharing all the commands you used to generate this result.

kubranarci commented 4 years ago

Raw reads are subsampled from 50X HG002 (GIAB). 1KG.phase3.AF001 is a AF filtered phase3 1000G project variant. This is how the processes look like:

VG CONSTRUCT AND VG IDS

tabix -p vcf 1KG.phase3.AF001.vcf.gz && (seq 1 22; echo X; echo Y) | parallel -j 24 "vg construct -C -t 1 -m 32 -R {} -r hs37d5.fa -v 1KG.phase3.AF001.vcf.gz > 1KG.phase3.AF001.vcf{}.vg" && vg ids -j $(for i in $(seq 1 22; echo X; echo Y); do echo 1KG.phase3.AF001.vcf$i.vg; done)

VG PRUNE AND VG INDEX vg prune -r --prune 1KG.phase3.AF001.vcf1.vg > 1KG.phase3.AF001.vcf1.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf10.vg > 1KG.phase3.AF001.vcf10.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf11.vg > 1KG.phase3.AF001.vcf11.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf12.vg > 1KG.phase3.AF001.vcf12.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf13.vg > 1KG.phase3.AF001.vcf13.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf14.vg > 1KG.phase3.AF001.vcf14.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf15.vg > 1KG.phase3.AF001.vcf15.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf16.vg > 1KG.phase3.AF001.vcf16.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf17.vg > 1KG.phase3.AF001.vcf17.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf18.vg > 1KG.phase3.AF001.vcf18.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf19.vg > 1KG.phase3.AF001.vcf19.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf2.vg > 1KG.phase3.AF001.vcf2.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf20.vg > 1KG.phase3.AF001.vcf20.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf21.vg > 1KG.phase3.AF001.vcf21.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf22.vg > 1KG.phase3.AF001.vcf22.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf3.vg > 1KG.phase3.AF001.vcf3.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf4.vg > 1KG.phase3.AF001.vcf4.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf5.vg > 1KG.phase3.AF001.vcf5.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf6.vg > 1KG.phase3.AF001.vcf6.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf7.vg > 1KG.phase3.AF001.vcf7.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf8.vg > 1KG.phase3.AF001.vcf8.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf9.vg > 1KG.phase3.AF001.vcf9.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcfX.vg > 1KG.phase3.AF001.vcfX.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcfY.vg > 1KG.phase3.AF001.vcfY.pruned.vg && vg index -x 1KG.phase3.AF001.xg -g 1KG.phase3.AF001.gcsa 1KG.phase3.AF001.vcf1.pruned.vg 1KG.phase3.AF001.vcf10.pruned.vg 1KG.phase3.AF001.vcf11.pruned.vg 1KG.phase3.AF001.vcf12.pruned.vg 1KG.phase3.AF001.vcf13.pruned.vg 1KG.phase3.AF001.vcf14.pruned.vg 1KG.phase3.AF001.vcf15.pruned.vg 1KG.phase3.AF001.vcf16.pruned.vg 1KG.phase3.AF001.vcf17.pruned.vg 1KG.phase3.AF001.vcf18.pruned.vg 1KG.phase3.AF001.vcf19.pruned.vg 1KG.phase3.AF001.vcf2.pruned.vg 1KG.phase3.AF001.vcf20.pruned.vg 1KG.phase3.AF001.vcf21.pruned.vg 1KG.phase3.AF001.vcf22.pruned.vg 1KG.phase3.AF001.vcf3.pruned.vg 1KG.phase3.AF001.vcf4.pruned.vg 1KG.phase3.AF001.vcf5.pruned.vg 1KG.phase3.AF001.vcf6.pruned.vg 1KG.phase3.AF001.vcf7.pruned.vg 1KG.phase3.AF001.vcf8.pruned.vg 1KG.phase3.AF001.vcf9.pruned.vg 1KG.phase3.AF001.vcfX.pruned.vg 1KG.phase3.AF001.vcfY.pruned.vg

VG MAP vg map --gcsa-name 1KG.phase3.AF001.gcsa --xg-name 1KG.phase3.AF001.xg --fastq HG002-NA24385-0.5x.p1.fastq --fastq HG002-NA24385-0.5x.p2.fastq > HG002-NA24385-0.5x.p2.gam

VG SURJECT vg surject --xg-name 1KG.phase3.AF001.xg HG002-NA24385-0.5x.p2.gam --threads 48 --bam-output > HG002-NA24385-0.5x.p2.bam

ekg commented 4 years ago

Would you please share a small set of reads that map with e.g. bwa mem, but not with vg map?

Do they map with vg map, but not surject?

On Fri, Dec 20, 2019 at 12:13 PM Kübra NARCI notifications@github.com wrote:

Raw reads are subsampled from 50X HG002 (GIAB). 1KG.phase3.AF001 is a AF filtered phase3 1000G project variant. This is how the processes look like:

VG CONSTRUCT AND VG IDS

tabix -p vcf 1KG.phase3.AF001.vcf.gz && (seq 1 22; echo X; echo Y) | parallel -j 24 "vg construct -C -t 1 -m 32 -R {} -r hs37d5.fa -v 1KG.phase3.AF001.vcf.gz > 1KG.phase3.AF001.vcf{}.vg" && vg ids -j $(for i in $(seq 1 22; echo X; echo Y); do echo 1KG.phase3.AF001.vcf$i.vg; done)

VG PRUNE AND VG INDEX vg prune -r --prune 1KG.phase3.AF001.vcf1.vg > 1KG.phase3.AF001.vcf1.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf10.vg > 1KG.phase3.AF001.vcf10.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf11.vg > 1KG.phase3.AF001.vcf11.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf12.vg > 1KG.phase3.AF001.vcf12.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf13.vg > 1KG.phase3.AF001.vcf13.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf14.vg > 1KG.phase3.AF001.vcf14.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf15.vg > 1KG.phase3.AF001.vcf15.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf16.vg > 1KG.phase3.AF001.vcf16.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf17.vg > 1KG.phase3.AF001.vcf17.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf18.vg > 1KG.phase3.AF001.vcf18.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf19.vg > 1KG.phase3.AF001.vcf19.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf2.vg > 1KG.phase3.AF001.vcf2.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf20.vg > 1KG.phase3.AF001.vcf20.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf21.vg > 1KG.phase3.AF001.vcf21.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf22.vg > 1KG.phase3.AF001.vcf22.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf3.vg > 1KG.phase3.AF001.vcf3.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf4.vg > 1KG.phase3.AF001.vcf4.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf5.vg > 1KG.phase3.AF001.vcf5.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf6.vg > 1KG.phase3.AF001.vcf6.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf7.vg > 1KG.phase3.AF001.vcf7.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf8.vg > 1KG.phase3.AF001.vcf8.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcf9.vg > 1KG.phase3.AF001.vcf9.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcfX.vg > 1KG.phase3.AF001.vcfX.pruned.vg && vg prune -r --prune 1KG.phase3.AF001.vcfY.vg > 1KG.phase3.AF001.vcfY.pruned.vg && vg index -x 1KG.phase3.AF001.xg -g 1KG.phase3.AF001.gcsa 1KG.phase3.AF001.vcf1.pruned.vg 1KG.phase3.AF001.vcf10.pruned.vg 1KG.phase3.AF001.vcf11.pruned.vg 1KG.phase3.AF001.vcf12.pruned.vg 1KG.phase3.AF001.vcf13.pruned.vg 1KG.phase3.AF001.vcf14.pruned.vg 1KG.phase3.AF001.vcf15.pruned.vg 1KG.phase3.AF001.vcf16.pruned.vg 1KG.phase3.AF001.vcf17.pruned.vg 1KG.phase3.AF001.vcf18.pruned.vg 1KG.phase3.AF001.vcf19.pruned.vg 1KG.phase3.AF001.vcf2.pruned.vg 1KG.phase3.AF001.vcf20.pruned.vg 1KG.phase3.AF001.vcf21.pruned.vg 1KG.phase3.AF001.vcf22.pruned.vg 1KG.phase3.AF001.vcf3.pruned.vg 1KG.phase3.AF001.vcf4.pruned.vg 1KG.phase3.AF001.vcf5.pruned.vg 1KG.phase3.AF001.vcf6.pruned.vg 1KG.phase3.AF001.vcf7.pruned.vg 1KG.phase3.AF001.vcf8.pruned.vg 1KG.phase3.AF001.vcf9.pruned.vg 1KG.phase3.AF001.vcfX.pruned.vg 1KG.phase3.AF001.vcfY.pruned.vg

VG MAP vg map --gcsa-name 1KG.phase3.AF001.gcsa --xg-name 1KG.phase3.AF001.xg --fastq HG002-NA24385-0.5x.p1.fastq --fastq HG002-NA24385-0.5x.p2.fastq > HG002-NA24385-0.5x.p2.gam

VG SURJECT vg surject --xg-name 1KG.phase3.AF001.xg HG002-NA24385-0.5x.p2.gam --threads 48 --bam-output > HG002-NA24385-0.5x.p2.bam

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/vgteam/vg/issues/2581?email_source=notifications&email_token=AABDQEPBEV6VJSQ72YTBABDQZSSEFA5CNFSM4J52KPTKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEHMU5RA#issuecomment-567889604, or unsubscribe https://github.com/notifications/unsubscribe-auth/AABDQEOCEYMVI5QL6D54LHDQZSSEFANCNFSM4J52KPTA .

kubranarci commented 4 years ago

I need to align the reads using bwa mem, I will update you as soon as possible.

kubranarci commented 4 years ago

Do you think in vg index providing input vg files not sorted by chr just like in my example will cause a problem?

ekg commented 4 years ago

Is there something with the flags that is wrong?

Try just mapping the reads from one chromosome to a graph made only for that chromosome.

On Fri, Dec 20, 2019, 13:21 Kübra NARCI notifications@github.com wrote:

@ekg https://github.com/ekg, using bwa mem %100 of the reads were mapped according to GATK CollectAlignmentSummaryMetrics run.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/vgteam/vg/issues/2581?email_source=notifications&email_token=AABDQENRRLNI3RQWVFIBYXLQZS2DXA5CNFSM4J52KPTKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEHMZEOA#issuecomment-567906872, or unsubscribe https://github.com/notifications/unsubscribe-auth/AABDQEPN3EQTG55EIKMOLADQZS2DXANCNFSM4J52KPTA .

ekg commented 4 years ago

Also 100% mapping rate is very weird! Are you sure it's not actually ~99% or so?

On Fri, Dec 20, 2019, 15:03 Erik Garrison erik.garrison@gmail.com wrote:

Is there something with the flags that is wrong?

Try just mapping the reads from one chromosome to a graph made only for that chromosome.

On Fri, Dec 20, 2019, 13:21 Kübra NARCI notifications@github.com wrote:

@ekg https://github.com/ekg, using bwa mem %100 of the reads were mapped according to GATK CollectAlignmentSummaryMetrics run.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/vgteam/vg/issues/2581?email_source=notifications&email_token=AABDQENRRLNI3RQWVFIBYXLQZS2DXA5CNFSM4J52KPTKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEHMZEOA#issuecomment-567906872, or unsubscribe https://github.com/notifications/unsubscribe-auth/AABDQEPN3EQTG55EIKMOLADQZS2DXANCNFSM4J52KPTA .

kubranarci commented 4 years ago

I have just completed mapping with bwa mem, the result shows that the read number flagged as 4 is 58432 while it was 21697376 for vg map. As you see the difference is huge. Now,I am mapping one chr to the vg graph. Yet can be please overview my pipeline I believe something is missed effecting mapping.

kubranarci commented 4 years ago

Do you think in vg index providing input vg files not sorted by chr just like in my example will cause a problem?

can you also answer this question please?

ekg commented 4 years ago

Can you build and index the graph for one chromosome and map reads which you expect to map to it?

What happens when you run make test in the vg source repo test directory?

As for the order, all that matters is that you maintain the same order of inputs.

On Mon, Dec 23, 2019, 11:08 Kübra NARCI notifications@github.com wrote:

Do you think in vg index providing input vg files not sorted by chr just like in my example will cause a problem?

can you also answer this question please?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/vgteam/vg/issues/2581?email_source=notifications&email_token=AABDQEPJPQJFEDTIYJYRRUDQ2CEZLA5CNFSM4J52KPTKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEHQYR5I#issuecomment-568428789, or unsubscribe https://github.com/notifications/unsubscribe-auth/AABDQELZYXB6MIJAKC44RXLQ2CEZLANCNFSM4J52KPTA .

ekg commented 4 years ago

Perhaps you should try to provide the files in the same order that you compacted their IDs. I don't know why that would cause the error though.

It's key to verify that a simple example works first.

On Mon, Dec 23, 2019, 11:35 Erik Garrison erik.garrison@gmail.com wrote:

Can you build and index the graph for one chromosome and map reads which you expect to map to it?

What happens when you run make test in the vg source repo test directory?

As for the order, all that matters is that you maintain the same order of inputs.

On Mon, Dec 23, 2019, 11:08 Kübra NARCI notifications@github.com wrote:

Do you think in vg index providing input vg files not sorted by chr just like in my example will cause a problem?

can you also answer this question please?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/vgteam/vg/issues/2581?email_source=notifications&email_token=AABDQEPJPQJFEDTIYJYRRUDQ2CEZLA5CNFSM4J52KPTKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEHQYR5I#issuecomment-568428789, or unsubscribe https://github.com/notifications/unsubscribe-auth/AABDQELZYXB6MIJAKC44RXLQ2CEZLANCNFSM4J52KPTA .

kubranarci commented 4 years ago

Hi, I just tried to construct a graph for chr10 and mapped (without vg ids and prune steps ), it seems to be aligned well. Therefore, there should be something in the workflow, maybe in vg index (maybe difference in order as you said). Now I will try to make it for all chromosomes providing the same order. If you see anything wrong with the provided workflow, please enlight me :) Thanks for your help.

kubranarci commented 4 years ago

I do also realize that I am indexing xg and gcsa files in one step (after vg prune), do you think that can cause a problem?

kubranarci commented 4 years ago

Indexing XG and GCSA files in separate steps and pruning graph before GCSA files only solved my problem. FYI.

ekg commented 4 years ago

Thank you for bearing with me on this one.

It seems that the indexing won't work in one step. That's not expected but I guess we haven't been testing it properly. Sorry for the trouble and thanks for the help.

On Wed, Dec 25, 2019, 15:08 Kübra NARCI notifications@github.com wrote:

Closed #2581 https://github.com/vgteam/vg/issues/2581.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/vgteam/vg/issues/2581?email_source=notifications&email_token=AABDQEJ7V3IIYCSGQYY66BTQ2NSNRA5CNFSM4J52KPTKYY3PNVWWK3TUL52HS4DFWZEXG43VMVCXMZLOORHG65DJMZUWGYLUNFXW5KTDN5WW2ZLOORPWSZGOVVKNV3Y#event-2908019439, or unsubscribe https://github.com/notifications/unsubscribe-auth/AABDQELYPOARK7ZNAZID7Y3Q2NSNRANCNFSM4J52KPTA .

kubranarci commented 4 years ago

No problem, thanks for your help.