Question: For bam and mutation signature I need coverage and number of observations
gravatar for kulvait
4.0 years ago by
Czech Republic
kulvait250 wrote:

Hi, I have variant calling pipeline for somatic mutations. I am trying to uncover mutations on low allele balance frequencies.

When I find mutation in particular sample, I would like to see if it is present in other samples. Unfortunately just comparing vcf files does not do the trick because of many possible filtering steps or lack of some calls. Thus for given signature of the mutation from vcf file I would like to obtain info regarding this particular mutation in number of bam files. In particular I am interested the most what is the coverage at the position of the mutation and what is the number of reads supporting the mutation.

For example I have 1.bam, 2.bam, 3.bam. Based on particular criteria (allele balance > 0.1, coverage > 500, allele observations > 10) i produce 1.vcf, 2.vcf, 3.vcf. Now when there is some variant in 1.bam based on these criteria, then in 2.bam there can be this mutation with allele balance 0.09, coverage 1000 and 9 allele observations (and it does not make it to 2.vcf) while in 3.bam there is 0 coverage and no allele observation. I would like to somehow extract this info about signatures of the mutations in 1.vcf in files 2.bam and 3.bam. What tool should I use and is there one?

Thanks, Vojtěch.

bam next-gen vcf • 1.5k views
ADD COMMENTlink modified 3.2 years ago by cpad011214k • written 4.0 years ago by kulvait250

Thanks Pierre, at least for SNP detection this is a good solution. Unfortunately it does work only for my unprocessed bam files. After processing them further by picard/GATK there is error

Nov 10, 2016 11:54:30 AM scan

SEVERE: htsjdk.samtools.SAMFormatException: Invalid GZIP header
htsjdk.samtools.SAMFormatException: Invalid GZIP header
    at htsjdk.samtools.util.BlockGunzipper.unzipBlock(
    at htsjdk.samtools.util.BlockGunzipper.unzipBlock(
    at htsjdk.samtools.util.BlockCompressedInputStream.inflateBlock(
    at htsjdk.samtools.util.BlockCompressedInputStream.readBlock(
    at htsjdk.samtools.BAMFileReader$BAMFileIndexIterator.getNextRecord(
    at htsjdk.samtools.BAMFileReader$BAMFileIterator.advance(
    at htsjdk.samtools.BAMFileReader$BAMFileIndexIterator.<init>(
    at htsjdk.samtools.BAMFileReader.createIndexIterator(
    at htsjdk.samtools.BAMFileReader.query(
    at htsjdk.samtools.SamReader$PrimitiveSamReaderToSamReaderAdapter.query(
    at htsjdk.samtools.SamReader$PrimitiveSamReaderToSamReaderAdapter.query(
    at com.github.lindenb.jvarkit.util.command.Command.instanceMainWithExceptions(
    at com.github.lindenb.jvarkit.util.command.Command.instanceMain(
    at com.github.lindenb.jvarkit.util.command.Command.instanceMainWithExit(

I can send you my bam file since I have no clue what is wrong. When i use samtools view it normally shows the reads.

Kind regards, Vojtech.

ADD REPLYlink written 4.0 years ago by kulvait250

Please use `ADD COMMENT to reply to earlier answers, as such this thread remains logically structured and easy to follow.

ADD REPLYlink written 4.0 years ago by WouterDeCoster44k

I suspect there is something wrong in your bam . Test your files with or / and

ADD REPLYlink written 4.0 years ago by Pierre Lindenbaum131k

Well, I have plenty of errors of the type

ERROR: Read name 700834F:16:HAYDCADXX:1:1206:13536:58893, Mate not found for paired read

Problem is that GATK tools did some kind of filtering and removed some reads without their mates. So I probably need to do that manually. Or is there some tool to do this?

ADD REPLYlink written 4.0 years ago by kulvait250

my tool would ignore this kind of problem, please run CheckTerminatorBlock too.

ADD REPLYlink written 4.0 years ago by Pierre Lindenbaum131k

Well, thank you I found the problem. There was files together with file.bam in directory


Index file file.bam.bai contained up to date index while the file.bai contained outdated index. I needed to remove incorrect file.bai and then the program works with correct index file.bam.bai,

thank you Pierre.

ADD REPLYlink modified 4.0 years ago • written 4.0 years ago by kulvait250
gravatar for Pierre Lindenbaum
4.0 years ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum131k wrote:

I've writtern a progam "FindAllCoverageAtPosition"

It scans a list of bam and produces a summary for a given position.

$ find ./  -type f -name "*.bam" |\
   java -jar dist/findallcoverageatposition.jar -p "chr2:1234"
ADD COMMENTlink written 4.0 years ago by Pierre Lindenbaum131k
gravatar for cpad0112
3.2 years ago by
cpad011214k wrote:
igvtools count -w 1 --bases --query chr1:100-1000 input.bam output.wig hg38
ADD COMMENTlink modified 3.2 years ago • written 3.2 years ago by cpad011214k
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 2298 users visited in the last hour