Get Allele Frequency of each SNP in RNASeq

Question: Get Allele Frequency of each SNP in RNASeq

4.9 years ago by

a1249m • 0

European Union

a1249m • 0 wrote:

For a set of RNASeq bam files aligned to hg19, I would like to calculate the number of reads mapping to all possible (high-confidence) alleles for each position. What is the best way to do this? I have tried using samtools mpileup and bcftools as was recommended in a few solutions Extract Allele At Particular Positions From Bam Files, but the output either gives me 'Error: Could not parse --min-ac g', or gives a vcf file which is deprecated and shows me a vcf file whose output I am not sure is correct. Here is a sample:

chr2 100002198 . T . 32.9952 . DP=1;MQ0F=0;AF1=0;AC1=0;DP4=1,0,0,0;MQ=60;FQ=-29.991 GT:PL:DV 0/0:0:0

chr2 100002199 . A . 32.9952 . DP=1;MQ0F=0;AF1=0;AC1=0;DP4=1,0,0,0;MQ=60;FQ=-29.991 GT:PL:DV 0/0:0:0

chr2 100002200 . T . 32.9955 . DP=1;MQ0F=0;AF1=0;AC1=0;DP4=1,0,0,0;MQ=60;FQ=-29.9906 GT:PL:DV 0/0:0:0

chr2 100002201 . C . 32.9955 . DP=1;MQ0F=0;AF1=0;AC1=0;DP4=1,0,0,0;MQ=60;FQ=-29.9906 GT:PL:DV 0/0:0:0

What I need is a file that shows me 1) each base position, 2) the reference nucleotide at that position, 3) the range of alleles shown by my sample at that position, and 4) the frequency of each allele.

sequencing snp rna-seq • 1.9k views

ADD COMMENT • link •

modified 4.9 years ago by Chris Miller ♦ 21k • written 4.9 years ago by a1249m • 0

4.9 years ago by

mark.ziemann • 1.2k

Australia/Mebourne/Geelong/Deakin

mark.ziemann • 1.2k wrote:

I can't help you with the bcftools error, but I have done something similar using VarScan which worked nicely.

ADD COMMENT • link written 4.9 years ago by mark.ziemann • 1.2k

4.9 years ago by

Chris Miller ♦ 21k

Washington University in St. Louis, MO

Chris Miller ♦ 21k wrote:

If you know the positions that you want to query, bam-readcount is probably the tool you're looking for.(https://github.com/genome/bam-readcount)

ADD COMMENT • link written 4.9 years ago by Chris Miller ♦ 21k

Thanks Chris, this what I was looking but I got confused from the documentation about the format of regions file. Can you please explain the right format it accepts?

ADD REPLY • link written 4.0 years ago by fizer • 30

Tab-delimited, one based coordinates

1    1234   1234
17   9876  9876

Column 4+ will be ignored

ADD REPLY • link written 4.0 years ago by Chris Miller ♦ 21k

Please log in to add an answer.

Similar posts • Search »