Question: Bam to nucleotide frequencies
1
gravatar for auser_27
5.1 years ago by
auser_2720
Canada
auser_2720 wrote:

I have whole genome sequences of several clinical isolates and I want to  get the nucleotide frequencies (% A,C,T,G) at a given position.  Among novel SNPs we could discover, there are some very specific positions that I want to look at. There are enough clinical isolates that I don't wish to this manually using IGV.

I can do this, but it's extremely inelegant and I think that there is some more elegant way that I am missing. I have summarize that there are probably two approaches

Approach 1:

1 . Run the basic samtools SNP pipeline, outputting all positions without a VCF ( remove -V) and outputting all possible reference alleles. 

2. Write a script that parses the info column in the VCF file, get the depth info for each nucleotide, calculate the frequencies (for my variants of interest only); luckily I already have this script.

 

Approach 2: 

1. Use bedtools and convert that bam to fastq

2. Use seqtk and convert fastq to fasta

3. Use bedtools nuc command, specifying position of interest

 

I feel like this is too common a task to require the use of so many tools and file format conversions and even custom scripts, so I ask : what's a better way to do this?

snp sequence next-gen • 2.9k views
ADD COMMENTlink modified 4.4 years ago by Biostar ♦♦ 20 • written 5.1 years ago by auser_2720
3
gravatar for swbarnes2
5.1 years ago by
swbarnes26.5k
United States
swbarnes26.5k wrote:

What about parsing the relevant line from samtools mpileup?

ADD COMMENTlink written 5.1 years ago by swbarnes26.5k

I think using mpileup with a given genomic coordinate is the simplest solution. It supports filters by MAPQ and you can specify genomic coordinates (just be careful about the off-by-one error if youre zero-counting vs one-counting bases). 

ADD REPLYlink written 5.1 years ago by karl.stamm3.5k

This was the fastest way to do it.

ADD REPLYlink written 5.1 years ago by auser_2720
2
gravatar for Adrian Pelin
5.1 years ago by
Adrian Pelin2.3k
Canada
Adrian Pelin2.3k wrote:

Install popoolation, and run snp-diff script. You will get a table, a column for chromosome, a column for nucleotide position in that chromosome, and then columns for each of your samples that contain allele frequencies.

ADD COMMENTlink written 5.1 years ago by Adrian Pelin2.3k
1
gravatar for Sean Davis
5.1 years ago by
Sean Davis25k
National Institutes of Health, Bethesda, MD
Sean Davis25k wrote:

In R, the bam2R function in the deepSNV package is pretty useful for this.  

ADD COMMENTlink written 5.1 years ago by Sean Davis25k
0
gravatar for poisonAlien
5.1 years ago by
poisonAlien2.8k
Asgard
poisonAlien2.8k wrote:

bam-readcount is one option.

ADD COMMENTlink written 5.1 years ago by poisonAlien2.8k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 531 users visited in the last hour