How to extract sequences from bam file using variant position
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?
0 answers
No answers yet.
Log in to answer this question.
More posts like this
-
Separate vcf file creation for matched tumor-normal samples
written by Gene_MMP8
-
How to extract protein coding sequence from fasta file?
written by doubleA •
-
why do bed files not contain information about the reference?
written by Marvin
-
Correct statistical test to determine the significance of nucleotides present
written by Gene_MMP8
-
Difference in chromosome numbering
written by Gene_MMP8
-
Similarity between reference files
written by Gene_MMP8
-
Finding cancer gene sequences
written by Gene_MMP8
-
Is there any location encoding available for the unwrapped DNA?
written by alireza.mdi •
-
Merging/Concatenating Vcf Files
written by bioinfo
-
Vcf Locations In Consensus Sequence
written by Nupur G •
You can do that using the VCF and ref fasta files. Why involve the BAM?
If s/he does want to just extract reads from the BAM, one can simply do this:
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.
Just use
bedtools getfasta: http://bedtools.readthedocs.io/en/latest/content/tools/getfasta.htmlCreate a custom BED file with start and end co-ordinates 5 bp from your variant.
Thanks. This is exactly what I needed. Just to be clear, the fasta file is the consensus sequence I generated, right?
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 faidxprior 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
headcommand in BASH) in order to learn the true behaviour ofbedtools getfastaThe 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.
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 getfastacommand on that. If there are indels, it may upset the coordinates.You may want to check the behaviour of vcf2fq in relation to indels.
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?
I want to extract the said bases around the variants from the bam file. Can I use the ref.fa here? @Ram
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.
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
How do you know there's a 10bp indel in your BAM file without any downstream processing?
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?
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.
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
Second BAM? You already have the VCF - that is the "another file" I was referring to, not another BAM.
Got it. Can you have a look at my edited answer above? Just to be clear. I included an exmaple
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.
I have posted there. Can you please reply? This is the link
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.