This is a test version of Biostars. For the public version, visit https://www.biostars.org.
BCFTOOLS: One-liner to get individuals with more than one heterozygote call in a region

Hello!

I have a regular multi-sample VCF from which I am looking to find individuals that have >1 HET call in the same gene. The vcf is chunked so relatively small in size (with the gene of interest in the chunk)/

Does anyone have a one liner for such a task with bcftools or even R?

Many thanks

bcftools
• 2,612 views • link
written by tacrolimus tacrolimus

multiple individuals merged into one across a single gene

VCF format does not work that way. How did you get to this from a real valid VCF?

• 0 views • link
written by Ram Ram

Yup - its a multi sample vcf chunked by region of the genome.

• 0 views • link
written by tacrolimus tacrolimus

Yup

Yes to what, exactly?

it's a multi-sample VCF

That is not the same thing as saying "multiple individuals merged into one". Please be clear - is it a multi-sample VCF or some sort of mutated file with multiple GTs per sample?

• 0 views • link
written by Ram Ram

Dear RamRS,

It is a multi-sample vcf. Whilst it may contain mutations I can assure you it has not been mutated.

• 0 views • link
written by tacrolimus tacrolimus

Thank you for clarifying that. Your initial statement "I have a vcf with multiple individuals merged into one across a single gene" is quite confusing as individuals are not "merged", they are merely part of the VCF, and I'm also assuming loci are not merged (the "across a single gene" part). What you have is a regular VCF from which you are looking to find individuals that have >1 HET call in the same gene.

• 0 views • link
written by Ram Ram

Correct! Apologies if that was too vague. I have amended by question accordingly.

• 0 views • link
written by tacrolimus tacrolimus

1 answer

using http://lindenb.github.io/jvarkit/BioAlcidaeJdk.html to create a BED tools + bcftools to get the variants. (you could get some false positives if an indel overlap a position )

bcftools view --targets-file <(java -jar dist/bioalcidaejdk.jar -e 'final Map<String,List<Interval>> sample2ctx = new HashMap<>();while(iter.hasNext()) { final VariantContext v= iter.next(); for(final Genotype g:v.getGenotypes()) { if(!g.isHet()) continue; List<Interval> list = sample2ctx.get(g.getSampleName()); if(list==null) { list = new ArrayList<>(); sample2ctx.put(g.getSampleName(),list); } list.add(new Interval(v)); } }for(final String sn: sample2ctx.keySet()) { List<Interval> list = sample2ctx.get(sn); if(list.size()<2) continue; for(Interval r:list) { print(r.getContig()); print("\t"); print(r.getStart()-1); print("\t"); print(r.getEnd()); print("\t"); println(sn); }}'  src/test/resources/rotavirus_rf.vcf.gz) src/test/resources/rotavirus_rf.vcf.gz
• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

That's quite the one liner! Thanks I'll give it a go. All the best.

• 0 views • link
written by tacrolimus tacrolimus

Pierre, I think the code you've written here is complicated enough to warrant an explainer. How is your java code picking regions to use with the --targets-file parameter?

• 0 views • link
written by Ram Ram

Right, I'm clear on the process sub part - just wondering how the bed file is created. In other words, can you describe the code in English please?

• 0 views • link
written by Ram Ram 
create an associative array [sample,list[interval]]
for each variant
   for each genotype
       if the genotype is het add the interval for this sample

for each sample
   print each interval as bed if the number is > 1
• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Log in to answer this question.

More posts like this