Genes of interest from SAM/BAM files

[ksw9]

Hi,
I am trying to compare paired-end Illumina data with standard MLST approaches for genotyping bacteria. What is the best way to extract known marker genes (using coordinates from a reference) from SAM/BAM files created by mapping samples to this reference?

Ultimately, I just want to compare the genotyping capacity of this shotgun data with traditional, MLST methods and ensure that we are getting good coverage of the MLST markers of choice. I would like to generate consensus sequences from short read data of the makers of interest and am unsure the best way to do this.

Thank you!

[crazyhottommy]

I am not sure what you want to do, but look at bedtools. it may help you.

[GenoMax]

1. http://seqanswers.com/forums/showthread.php?t=39766

Assuming that you are referring to extracting reads that map to certain marker genes: samtools view should allow you to pull out reads from specified gene regions.

2. http://seqanswers.com/forums/showthread.php?t=38969

Samtools mpileup should generate the consensus sequence.

[ksw9]

Thank you for your help!

[ksw9]

I am trying to generate a consensus sequence using the following command:

samtools mpileup -uf B31.fna Sample_Bbcap10_006_004_MLST.bam | bcftools view -cg - | /home/bioinfo/software/samtools/samtools-0.1.19/bcftools/vcfutils.pl vcf2fq > cns.fq

and struggling with the error:

Use of uninitialized value in length at /home/bioinfo/software/samtools/samtools-0.1.19/bcftools/vcfutils.pl line 544, <> line 57.


Does anyone have any insight on this? I'm just beginning with SAMtools and really appreciate the support.

[GenoMax]

Have you indexed your reference genome file (B31.fna) with samtools faidx? That is a fasta format file correct?

$ samtools faidx B31.fna

[ksw9]

Yes, it is a fasta and I did index with samtools faidx.

[shuoguo]

looks to me the error occurs at last step. does "vcfutils.pl" takes stdin?

[GenoMax]

@ksw9: This may be an obvious question but was the reference file used for generation of the indexes and the alignments to generate the BAM's? Based on the name it looks like that may be a 454 sequence file. What aligner did you use?

[ksw9]

I used BWA for mapping paired-end Illumina reads to the reference genome.

[ksw9]

Okay, I think my problem was in my original .bed file which now has contigs appropriately labeled.

I guess my question is this: if I have extracted reads mapping only to loci of interest using:
samtools view -h -L MLST.bed Sample_006_004.bam > Sample_006_004_MLST.sam

then should I create a new reference sequence including only the loci of interest in order to generate the consensus sequence only for the loci of interest (and not generate the consensus for the entire genome)?

Thank you for all of your help!

[shuoguo]

Okay, I think my problem was in my original .bed file which now has contigs appropriately labeled.

I guess my question is this: if I have extracted reads mapping only to loci of interest using:
samtools view -h -L MLST.bed Sample_006_004.bam > Sample_006_004_MLST.sam

then should I create a new reference sequence including only the loci of interest in order to generate the consensus sequence only for the loci of interest (and not generate the consensus for the entire genome)?

Thank you for all of your help!

Good to know, thanks!