Hi, Suppose I have a bam file and a vcf file containing variant calling result. I want to extract only reads with their mate that support variant allele in the vcf. It would be nice to get those reads in bam format. I tried googling such tools to do this and found like VariantBAM but it reports both reads that supporting and not supporting variant.
Thanks
•
written
by
pitithat •
I quickly wrote something : http://lindenb.github.io/jvarkit/Biostar322664.html
caveat: BAM files must be sorted with picard SortSam/queryname and variants are loaded in memory and only SNP are considered.
$ java -jar picard.jar SortSam I=src/test/resources/S1.bam O=query.bam SO=queryname
$ java -jar dist/biostar322664.jar -V src/test/resources/S1.vcf.gz query.bam
(...)
RF02_358_926_2:0:0_2:1:0_83 83 RF02 857 60 70M = 358 -569 GACGTGAACTATATAATTAAAATGGACAGAAATCTGCCATCAACAGCTAGATATATAAGACCTAATTTAC 2222222222222222222222222222222222222222222222222222222222222222222222 RG:Z:S1 NM:i:3 AS:i:55 XS:i:0
RF02_362_917_2:0:0_2:1:0_6f 147 RF02 848 60 70M = 362 -556 ATAAGGAATCACGTTAACTATATACTTAAAATGGACTGAAATCTGCCATCAACAGCTAGATATATAAGAC 2222222222222222222222222222222222222222222222222222222222222222222222 RG:Z:S1 NM:i:3 AS:i:55 XS:i:0
(...)
it's not fully tested, tell me if something is looks wrong.
Thank you much. But I still find that your program also report reads that don't support the variant allele.
•
written
by
pitithat •
can you please send a minimal VCF (one variant) and a minimal bam. Thanks. https://github.com/lindenb/jvarkit/issues/new
Just as a warning, strictly speaking, it is not possible to extract reads supporting a variant in a VCF since you don't know what reads the underlying variant caller has used. Mutect2, for example, discard reads failing some configurable options about base quality, mapping quality, flags etc plus reads are locally re-aligned. So the number of reads supporting a variant in the VCF output of Mutect2 is typically lower than the one you would "see" in a genome viewer. This gets even more complicated for structural variants. In fact, what my program does is simply to extract reads in a region if the read bases in that region do not match the reference.
Hello, I wanna make sure how your program work. Is this the command "filterVariantReads -r 1000 vars.*vcf" to do so? If it is, what will it give me and in what format?
Thanks
Log in to answer this question.
Hello,
have a look at this thread. It might be useful for you.
fin swimmer
Hello, Yeah, I found this thread however they don't give me the mate reads.
If it gives you the read names that support the variant, you can take this list as an input for FilterSamReads or a simple
grep. Depending on your aligner the mate have the same read name and should be returned as well.fin swimmer
This would be a nice tool indeed. For visualization, the program I've written, ASCIIGenome, has the option filterVariantReads that gets close to what you need if combined with the print command. It could be scripted and automated but it's not quite what you ask for.
I am going to tag Pierre Lindenbaum
He may already have something written to do this.
hi Pierre, Excellent!!! This is something I was trying to do a couple of weeks ago and settled for VariantBam. But now I see this discussion, I would like to ask if this is applicable for structural variants as well. As we all know, the vcf for structural variants is a bit different compared to SNVs. Here are how the vcf records look like
Getting the supporting reads for these variants is a little confusing. And at the moment, I am taking the reads from both breakpoints separately and then combining. This is a little laborious way but as I was on a tight deadline, I just did it manually. It would be nice to have a way to extract just supporting reads for such structural variants (helps visualization and validation).
Thanks in advance, Venkatesh Chellappa (Venki)
it's too complicated for now, but anyway, why would you need to extract the variant of a SV ? Isn't visualization enough to validate a variant ?
I need the "reads supporting the variant" to visualize in IGV. I am currently using an "evidence bam" that contains reads that support breakpoints on left and right ends of the structural variants.
so yeah, me manually extracting the reads based on the loci of breakpoints is too laborious and I want to know if anyone has working solution for this.
I am thinking of opening a new discussion!