To speed up variant calling (SNP / Indel) on a set of large BAM files I want to do the variant calling in paralell on a Sun Grid Engine cluster.
I already managed to split the BAMs by reference (chromosome) using https://github.com/pezmaster31/bamtools/ . Afterward I use vcf tools to concatenate the chromosome VCF files to one sample VCF file.
These bamfiles are still relatively big so I would like to split the BAMS further on large N polymer sequences in the reference. I already have a file with the location of these large N polymer sequences in our reference.
Is there a tool that already has the functionality to split BAM files by N polymer sequences?
Or what would be the easiest way to use this region information to split BAM files?
2 answers
You can use samtools: http://samtools.sourceforge.net/
In samtools view, there is a -L parameter that'll allow you to extract reads that overlap with a specified BED file:
-L FILE output alignments overlapping the input BED FILE
Convert your N-polymer coordinates into a bed file and extract reads that only overlap those coordinates with the -L parameter. You'll have to make a bed file that only has regions without the N-polymer coordinates to extract everything else. I don't think there is an option to exclude regions in samtools.
Hi Damian,
Your post is almost 3 years old but I am currently thinking about this problem.
If you split your chromosome into 1Kb for example, what about the reads that start in the 0-1Kb chunk and finish in the 1Kb-2Kb chunks? There are in both chunks right? Isn't that a problem? Doesn't it induce errors?
I am not sure how samtools handles the -L parameter or if it even still implements it. If samtools will grab any reads that intersects (totally or partially) the region defined the bed file, then it should be fine. The only problem is that you might get duplicated calls among consecutive region's vcfs which you can check for and take out.
Thanks for the answer.
Yes, then I want to merge my vcfs.
I was just worried that 2 vcfs could disagree about an event (like 0-1Kb chunk vcf says that position 1008 is a sequencing error and 1Kb-2Kb chunk vcf says that it is a SNP).
Log in to answer this question.
-
How to subset large BAM files specifically/ extract specific subsets?
written by ella •
-
Find reference genome regions spanned by only mapping quality 0 reads in multiple WGS samples
written by William
-
Changing output directory for bamtools split
written by phosphorus •
-
how to split a Bam by chromosome and unplaced scaffolds.
written by Maxine •
-
Best way to find intervals for parallelization of joint variant calling
written by William
-
variant calling with platypus
written by evelyn
-
unsuccessful split of bam file by samtools and bamtools
written by grobelsek.v •
-
Creating a cell by SNP matrix
written by dominicdhall •
-
Cannot install ATHLATES, issues with bamtools
written by darkmatter1996 •
-
Downsampling Bam Files
written by Allpowerde
Is I/O your limiting factor here? if not, just place a single bam on a network-accessible disk, and use samtools to extract the reads you need. (samtools view file.bam 2:1234-5678 | variantcaller)
I can't pipe the samtools output, I need to pass a bam file path to an 3rd party variant caller.
Not a direct answer but GATK 2 has support for reduced bams now so they might reduce your bam sizes for multi sample calling.
but of course if your issue is getting more chunks to parallelise then I think the samtools -L bedfile seems like a good idea