This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Extract consensus gene sequence (from gff file) from bam file

Hi,

I have a bam file from Illumina DNA-seq that I obtained by mapping reads to the reference genome using Bowtie. I want to extract the consensus sequence of a specific gene from the bam file irrespective of the reference. The gene coordinates are present in the gff file but I did not use this when mapping to the reference.

Thanks very much!

bam dna-seq gff
• 4,469 views • link
updated by Ram Ram • written by AW AW

You can use SAMtools to specify the chromosome/interval of the gene, and its variant-calling workflow to identify the differences b/t the reference and your data.

• 0 views • link
written by harold.smith.tarheel harold.smith.tarheel

Thanks! How do I go about doing this? I've only seen examples specifying one region but I want to extract the multiple exons and not intron?

• 0 views • link
written by AW AW

1 answer

I have tried the following but get an error

samtools view -b input.bam "chr1:14382-15521" > subset.bam
samtools sort subset.bam subset_sort
samtools mpileup -d 10000000 -uf ref.fna subset_sort.bam | bcftools view -cg - | vcfutils.pl vcf2fq > cns.fq

Error: Could not parse --min-ac g
[mpileup] 1 samples in 1 input files
Use of uninitialized value $l in numeric lt (<) at /bcftools/bin/vcfutils.pl line 566.
Use of uninitialized value $l in numeric lt (<) at /bcftools/bin/vcfutils.pl line 566.

Any ideas?

• 0 views • link
updated by Ram Ram • written by AW AW

Log in to answer this question.

More posts like this