SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Calculate depth of every base position Lv Ray Bioinformatics 5 10-17-2014 11:04 PM
base coverage per position on scaffolds boetsie Bioinformatics 5 01-23-2014 05:44 AM
base position cmccabe Bioinformatics 1 01-15-2014 06:08 AM
How to find out the base on a certain position shuang Bioinformatics 8 09-09-2011 06:45 AM
MAQ: finding the stop position of an alignment scirocco Bioinformatics 0 12-28-2008 11:33 AM

Reply
 
Thread Tools
Old 06-30-2015, 06:10 PM   #1
susanklein
Senior Member
 
Location: oceania

Join Date: Feb 2014
Posts: 115
Default finding the base at a position in a bam

Hi All,

I have seen some other threads on this subject, sorry I don't have the links. But none of them seem to resolve the problem. I am trying to do something very simple but no tools for bam files seem to have the facility. I have tried samtools and GATK to no avail.

All I need is to get the consensus read base at particular positions of the reference in a bam file. I am trying to check the state of SNPs from a particular bacterial SNP typing scheme in my NGS data. Not all of the positions will actually be a SNP, in which cae they have the same base as the reference, and I want that reported, so .vcf is no good to me.

Greatly appreciate any suggestions. I suppose I can just use samtools to make the vcf and then if a position is not listed then assume it is the reference base. I do find samtools almost impenetrable.. would this line be appropriate?

samtools mpileup -I -uBf ref.fasta 1_sorted.bam | bcftools view -bcg - > 1.bcf && bcftools view 1.bcf > 1.vcf

but I don't understand the output vcf file because it seems to have a '.' for the reference base and then the proper A T C or G for the 'alt' base(s).

Thanks for any help,

S.
susanklein is offline   Reply With Quote
Old 06-30-2015, 10:33 PM   #2
SNPsaurus
Registered Vendor
 
Location: Eugene, OR

Join Date: May 2013
Posts: 521
Default

I do think the mpileup command is what you want. I am not sure what your pipeline is making, but it sounds like the pileup format (. or , for the reference base with ATCG with alternates).
You can output a vcf with genotype compared to the reference for each position:
samtools mpileup -gu -f ref.fasta sorted.bam | bcftools call -c - > all.vcf
The genotype will be 0/0 for homozygous reference, 0/1 for heterozygous with the alternate listed allele, and 1/1 for homozygous alternate.

bcftools call -cv would output the variants only.

You might want to use bcftools consensus to create a modified fasta file with your variants:
http://www.htslib.org/doc/bcftools.html#consensus
__________________
Providing nextRAD genotyping and PacBio sequencing services. http://snpsaurus.com
SNPsaurus is offline   Reply With Quote
Old 06-30-2015, 10:39 PM   #3
susanklein
Senior Member
 
Location: oceania

Join Date: Feb 2014
Posts: 115
Default

Thanks SNPsaurus, I will try that. But I've now found that the published SNP locations I have been using have mistakes so have to go back a step..

S.
susanklein is offline   Reply With Quote
Old 07-01-2015, 11:54 AM   #4
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

samtools mpileup should also take as a commend line entry a position, or a bed file of positions, so you don't need to generate a line for every base in the genome.
swbarnes2 is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 12:13 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO