Fixes for regional BAM re-alignment with bgzipped inputs: glia: out_of_range and tabix errors

[chapmanb]

Eric;
This fixes two issues I ran into while using glia for realignment included
called variants.

• It avoid erroring out when tabix region retrieval has no variants. The current
  logic would fail both on tabix retrieval due to a return problem with vcflib
  (fixed at Complementary fix for glia regional realignment around variants vcflib/vcflib#28) and in regions which don't
  have any variants where setRegion returns false. The fix checks the logic but
  avoids erroring out, instead only iterating over variants if they are present.
• Avoids std::out_of_range errors from referenceSequence.at(sp) when an
  alignment falls over the end of reference sequence. The fix counts these as
  mismatches.

Here's a small reproducible test case that triggers these errors.

wget https://github.com/chapmanb/bcbio.variation.plus/raw/master/test/data/ensemble/NA12878-10.bam
wget https://raw.github.com/chapmanb/bcbio.variation.plus/master/test/data/ensemble/chr10-start.fa
wget https://s3.amazonaws.com/chapmanb/union-10_250000_400000.vcf.gz
tabix -p vcf union-10_250000_400000.vcf.gz
samtools view -bu NA12878-10.bam | glia -d -Rr -w 1000 -S 200 -Q 200 -G 4 -f chr10-start.fa -v union-10_250000_400000.vcf.gz

Thanks again for this tool. I'm excited to have a minimal test case running with
it and looking forward to exploring more.

[ekg]

Hi Brad,

In trying to run the test I hit a snag. The BAM file (https://github.com/chapmanb/bcbio.variation.plus/raw/master/test/data/ensemble/NA12878-10.bam) does not appear to exist, or the URL has changed. Could you please update this and I'll run the tests.

[chapmanb]

Erik;
Sorry about this. I moved these over to another repository. The new wgets are:

wget https://github.com/chapmanb/bcbio.variation.recall/raw/master/test/data/ensemble/NA12878-10.bam
wget https://raw.github.com/chapmanb/bcbio.variation.recall/master/test/data/ensemble/chr10-start.fa

Please let me know if you run into anything else. Thanks much for looking at these.

[ekg]

Great, I've just downloaded them.

On Fri, Jan 31, 2014 at 7:22 AM, Brad Chapman [email protected]:

Erik;
Sorry about this. I moved these over to another repository. The new wgets
are:

wget https://github.com/chapmanb/bcbio.variation.recall/raw/master/test/data/ensemble/NA12878-10.bam
wget https://raw.github.com/chapmanb/bcbio.variation.recall/master/test/data/ensemble/chr10-start.fa

Please let me know if you run into anything else. Thanks much for looking
at these.

Reply to this email directly or view it on GitHubhttps://github.com//pull/1#issuecomment-33789855
.

[yifangt]

I was wondering where the "fastq.h" header is or, it may not be needed?

[jollymrt]

Hi Erik,

I am trying to run glia to realign a bam file around indels but getting the following error "Command not found : glia Did you mean ...."
how to resolve this issue. I am currently in glia directory and have run make and cmake on both bam and glia. I dont see any executable or dll files after running make. Can you please help.

Regards
Jolly

[ekg]

Glia needs to be in a directory in your PATH environment variable. This is a Unix level issue that you can resolve by putting the executable in an expected place or adding the current directory with it to the path. I should reiterate, use vg rather than glia. It has had hundreds of times the use and developmental effort. https://github.com/vgteam/vg

…

On Fri, Aug 10, 2018, 00:07 jollymrt ***@***.***> wrote: Hi Erik, I am trying to run glia to realign a bam file around indels but getting the following error "Command not found : glia Did you mean ...." how to resolve this issue. I am currently in glia directory and have run make and cmake on both bam and glia. I dont see any executable or dll files after running make. Can you please help. Regards Jolly — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#1 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAI4EYBaeKE2Js46IJUchaartmnwpkvZks5uPLKugaJpZM4BVx3Z> .

[jollymrt]

Hi Erik, Thanks for your response. Basically I need to run a realignment on the bam files that will then be used for our in-house variant caller. I looked up vg but it does not have a separate realignment tool or may be I did not looked up thoroughly. Can you please point if this exists in vg ? Also in glia, I did not see any executable file in the package after running make. Will the file be made in a different folder? Again thanks for your help. Regards On Fri, Aug 10, 2018 at 12:43 AM, Erik Garrison <[email protected]> wrote:

…

Glia needs to be in a directory in your PATH environment variable. This is a Unix level issue that you can resolve by putting the executable in an expected place or adding the current directory with it to the path. I should reiterate, use vg rather than glia. It has had hundreds of times the use and developmental effort. https://github.com/vgteam/vg On Fri, Aug 10, 2018, 00:07 jollymrt ***@***.***> wrote: > Hi Erik, > > I am trying to run glia to realign a bam file around indels but getting > the following error "Command not found : glia Did you mean ...." > how to resolve this issue. I am currently in glia directory and have run > make and cmake on both bam and glia. I dont see any executable or dll files > after running make. Can you please help. > > Regards > Jolly > > — > You are receiving this because you commented. > Reply to this email directly, view it on GitHub > <#1 (comment)>, or mute the > thread > <https://github.com/notifications/unsubscribe-auth/ AAI4EYBaeKE2Js46IJUchaartmnwpkvZks5uPLKugaJpZM4BVx3Z> > . > — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#1 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ACaFZjwHlZdqh2CX1ZExNnsxnQhpP5vzks5uPTmRgaJpZM4BVx3Z> .

-- Jolly Shrivastava Ph.D (Genetics, Genomics and Bioinformatics) University of California, Riverside

[ekg]

@jollymrt let's continue this conversation here: vgteam/vg#1823