This is a test version of Biostars. For the public version, visit https://www.biostars.org.
extracting Allele Read Counts

Dear all,

please could you advise : given a (tumor) BAM file and a (germline) VCF file, what tool shall i use in order to extract the Allele Read Counts for each heterozygous SNP (from the vcf file) ?

many thanks,

-- bogdan

snp titan cnv
• 4,489 views • link
updated by Pierre Lindenbaum Pierre Lindenbaum • written by Bogdan Bogdan

Depending on the variant caller you are using. For some of them you will see the following two fields under FORMAT column:

##FORMAT=<ID=RD,Number=1,Type=Integer,Description="Depth of reference-supporting bases (reads1)">
##FORMAT=<ID=AD,Number=1,Type=Integer,Description="Depth of variant-supporting bases (reads2)">
• 0 views • link
written by iraun iraun

The easiest way to do this is to extract the sites using bcftools view -v snps -g het. You can add -m2 -M2 if you are only interested in biallelic site. Then you can use allelecounter.py (https://github.com/secastel/allelecounter) to extract the read counts at ref and alt allele.

• 0 views • link
written by prasundutta87 prasundutta87

thank you for your suggestions ;)

• 0 views • link
written by Bogdan Bogdan

1 answer

a combination of some programs I've written.

  • Extract the heterozygous sites with vcffilterjs.
  • use awk to extract the 'CHROM:POS' of those site.
  • for each position , call FindAllCoverageAtPosition with the specified bam
  • sort + uniq to get one uniq header.
    java -jar jvarkit-git/dist/vcffilterjs.jar -e 'variant.getGenotype(0).isHet()' input.vcf | grep -v "^#" |\
awk '{printf("%s:%d\n",$1,$2);}' | while read  S; do echo "input.bam" | java -jar /jvarkit-git/dist/findallcoverageatposition.jar -p "${S}" ; done | \
LC_ALL=C sort | uniq | column -t


#File   CHROM      POS  SAMPLE  DEPTH  M    I   D  N  S   H  P  EQ  X  Base(A)  Base(C)  Base(G)  Base(T)  Base(N)  Base(^)  Base(-)
S1.bam  rotavirus  130  S1      328    328  0   0  0  52  0  0  0   0  33       0        0        295      0        0        0
S1.bam  rotavirus  232  S1      363    363  35  0  0  56  0  0  0   0  25       0        0        338      0        35       0
S1.bam  rotavirus  267  S1      315    315  1   1  0  44  0  0  0   0  0        296      19       0        0        1        1
S1.bam  rotavirus  961  S1      327    327  0   0  0  45  0  0  0   0  0        0        29       298      0        0        0

(later) pushed a fresh version on github a few minutes ago. there was a problem in the old findallcoverageatposition (all files were 'cram')

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

thank you Pierre for your suggestions.

• 0 views • link
written by Bogdan Bogdan

Log in to answer this question.

More posts like this