Question: (Closed) How to extract sequences from bam file using variant position
0
gravatar for banerjeeshayantan
17 months ago by
banerjeeshayantan140 wrote:

I have a bam file and I have called variants(only SNPs) by aligning it with the reference hg38. Now for each variant location x, I want to extract DNA sequence from x-5 to x+5.
I created a consensus sequence and it starts with a lot of Ns. So how should I count the coordinates for the start of the genome, or location 1.
Say I have a variant at location 1233 of chromosome 2. So starting from what location of the consensus sequence, can I reach 1233?

sequencing next-gen • 1.3k views
ADD COMMENTlink written 17 months ago by banerjeeshayantan140

You can do that using the VCF and ref fasta files. Why involve the BAM?

ADD REPLYlink written 17 months ago by RamRS24k

If s/he does want to just extract reads from the BAM, one can simply do this:

samtools view -h MyAlignedReads.bam chr1:1000-2000
ADD REPLYlink written 17 months ago by Kevin Blighe50k

I want to extract a string of 11 nucleotides x-5:x+5 where x is the location of the variant. Can I use the above command? I don't want the entire read as such. Just a part of it.

ADD REPLYlink written 17 months ago by banerjeeshayantan140

Just use bedtools getfasta: http://bedtools.readthedocs.io/en/latest/content/tools/getfasta.html

Create a custom BED file with start and end co-ordinates 5 bp from your variant.

ADD REPLYlink written 17 months ago by Kevin Blighe50k

Thanks. This is exactly what I needed. Just to be clear, the fasta file is the consensus sequence I generated, right?

ADD REPLYlink written 17 months ago by banerjeeshayantan140

If that's the FASTA sequence against which you aligned your reads to produce your BAM and variant calls, then, yes, and you should index that consensus FASTA with samtools faidx prior to using bedtools.

Just test it with a single variant first in order to confirm that it does indeed work. You should even test it by just extracting the first few bases in the reference genome (and then simple look at those bases in the consensus FASTA via the head command in BASH) in order to learn the true behaviour of bedtools getfasta

ADD REPLYlink modified 17 months ago • written 17 months ago by Kevin Blighe50k

The consensus sequence was created using
samtools mpileup -d8000 -uf hg38.fa fixed_Realigned_tumor.bam | bcftools call -c | vcfutils.pl vcf2fq > cnsnew.fq
cnsnew.fq is my consensus file. I am generating this as without a consensus sequence I might be left with multiple strings (x-5:x+5) for each variants. I want a single string for any given mutation.

ADD REPLYlink written 17 months ago by banerjeeshayantan140

Based on your series of commands, you're just 'imputing' the hg38.fa reference with your called variants - is that right? In that case, you should be able to just index the cnsnew.fq file and then run the bedtools getfasta command on that. If there are indels, it may upset the coordinates.

You may want to check the behaviour of vcf2fq in relation to indels.

ADD REPLYlink modified 17 months ago • written 17 months ago by Kevin Blighe50k

Can you kindly go through the comments below posted by me today? I was wondering what would happen if I use only the ref.fa and there are some large INDELs in my tumour dataset. Can I create a ref.fa based on the variant information?

ADD REPLYlink written 17 months ago by banerjeeshayantan140

I want to extract the said bases around the variants from the bam file. Can I use the ref.fa here? @Ram

ADD REPLYlink modified 17 months ago • written 17 months ago by banerjeeshayantan140

BAM has more information than you need. You already have the positions, the actual change as well as the reference. What is the information you need that is in the BAM file but not in the VCF or the ref fa? Plus, it is more expensive to query BAMs, so avoid it unless you have really good reason to do that.

ADD REPLYlink modified 17 months ago • written 17 months ago by RamRS24k

Apologies for a late question. So I understand the point you are trying to make about using ref.fa instead of the bam file. But what if there are some large INDELS that might mess up my extracted sequence. Say, a 10bp INDEL is present in my tumour bam. I am trying to extract some sequence, a part of which falls within this INDEL. So what then? Can I change my ref.fa to account for INDELS in my bam file? Is there any tool available? Thanks

ADD REPLYlink written 17 months ago by banerjeeshayantan140

How do you know there's a 10bp indel in your BAM file without any downstream processing?

ADD REPLYlink written 17 months ago by RamRS24k

I am just entertaining the possibility. If there are INDELS of more than some x bp length, then I would have to change my ref.fa accordingly, right?

ADD REPLYlink written 17 months ago by banerjeeshayantan140

By the time you're done finding those indels, you'll have another file that will be used in place of the BAM file. BAM does not let you access the information you want without considerable computational cost for either use case. In case you have significant structural variants, you should be able to use the SV-calling results (the VCF) to produce a new reference. In no case can I imagine a BAM file being the only data source for your problem.

ADD REPLYlink written 17 months ago by RamRS24k

Thanks for your explanation and patience. I understand your point regarding computational cost and that querying a FASTA file is easier, but Can you please elaborate on "By the time you're done finding those indels, you'll have another file that will be used in place of the BAM file." Just to be clear, I start with an aligned normal/tumor pair sample. I find the variants by comparing normal vs tumor. This produces a VCF file which contains both SNP and INDEL information. So where is this second bam file coming from?
Eg:
Ref: ATGGCGCA
Tumor bam: ACGTTACTCA
Now I am using the first ref file to extract two bases to the left and right of the 6th variant. So in the ref file I will get GC and CA, but in the tumour file it is AC and CA. So how to proceed here?

Thanks again

ADD REPLYlink modified 17 months ago • written 17 months ago by banerjeeshayantan140

Second BAM? You already have the VCF - that is the "another file" I was referring to, not another BAM.

ADD REPLYlink written 17 months ago by RamRS24k

Got it. Can you have a look at my edited answer above? Just to be clear. I included an exmaple

ADD REPLYlink modified 17 months ago • written 17 months ago by banerjeeshayantan140

Is the neighborhood discussion not happening in a different thread already? Please discuss the neighborhood topic in that thread. Please explain (over there) why you wish to use the BAM and not the VCF.

ADD REPLYlink written 17 months ago by RamRS24k

I have posted there. Can you please reply? This is the link

ADD REPLYlink modified 17 months ago • written 17 months ago by banerjeeshayantan140

Hello banerjeeshayantan!

Questions similar to yours can already be found at:

We have closed your question to allow us to keep similar content in the same thread.

If you disagree with this please tell us why in a reply below. We'll be happy to talk about it.

Cheers!

PS: Do not open multiple posts for thee same question.

ADD REPLYlink modified 17 months ago • written 17 months ago by RamRS24k
Please log in to add an answer.
The thread is closed. No new answers may be added.

Help
Access

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