SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
calling SNPs in haploid genomes d17 Genomic Resequencing 27 10-01-2012 06:06 AM
minimum depth variant calling samtools/gatk m_elena_bioinfo Bioinformatics 1 12-06-2011 08:31 AM
variant calling kjaja Bioinformatics 1 11-04-2011 07:16 AM
Samtools variant calling questions Chiel Bioinformatics 2 06-07-2011 09:10 AM
variant calling using samtools -v- bcftools ksc Bioinformatics 2 04-13-2011 06:44 AM

Reply
 
Thread Tools
Old 06-25-2012, 01:24 AM   #1
HTrewby
Junior Member
 
Location: Glasgow UK

Join Date: Jun 2012
Posts: 1
Default haploid variant calling in SAMtools

We’re using the standard SAMtools SNP calling protocol (samtools mpileup -uf ref.fa aln1.bam aln2.bam | bcftools view -bvcg - > var.raw.bcf) to look for sequence variants in bacterial genomes.

I’m really new to this, but it sounds like the protocol is geared towards diploid organisms by default. I hoping someone might know if there’s anything that should be changed to use it on haploid genomes? The –s or –p options in SAMtools mpileup look like they might be useful, but I’ve no feel for how to apply them correctly for this. Any help would be really appreciated!
HTrewby is offline   Reply With Quote
Old 06-25-2012, 08:18 AM   #2
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

I use samtools on bacteria all the time, and it's not really a problem. Remember that just because your organism is haploid doesn't mean that you have been giving a clonal culture to analyze.
swbarnes2 is offline   Reply With Quote
Old 12-17-2012, 02:32 PM   #3
dictseq
Junior Member
 
Location: Midwestern US

Join Date: Dec 2012
Posts: 1
Default

We found that the following workflow works well for variant calling in haploids using samtools. To generate a consensus sequence from the realigned BAM file in FASTQ format, we used:

samtools mpileup -uf ref.fasta realigned.bam | bcftools view -cg -s sample.txt - | perl vcfutils.pl vcf2fq > consensus.fq

where sample.txt is a plain text file with the sample name in the first column and the ploidy level in the second column. For example,

sample_name 1

The columns should be tab-delimited.
dictseq is offline   Reply With Quote
Old 12-18-2012, 08:23 AM   #4
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

Keep in mind that vcfutils vcf2fq will not handle indels. It just creates a small window of lowercase letters around the putative indel.
swbarnes2 is offline   Reply With Quote
Old 12-19-2012, 04:11 AM   #5
ekg
Member
 
Location: Boston, MA

Join Date: Apr 2010
Posts: 36
Default

Why not use freebayes? It's designed to support haploids and polyploids (use the --ploidy option).
ekg is offline   Reply With Quote
Old 12-19-2012, 04:15 AM   #6
ekg
Member
 
Location: Boston, MA

Join Date: Apr 2010
Posts: 36
Default

Also, you can use vcfgeno2haplo (https://github.com/ekg/vcflib/blob/m...geno2haplo.cpp) with an arbitrarily large window size to generate a consensus fasta sequence from the output.

If you input a file with one sample, the recorded consensus sequence will be the ALT in the single VCF record in the output.

It works for indels as well as SNPs. It does not generate quality estimates.
ekg is offline   Reply With Quote
Old 12-19-2012, 04:04 PM   #7
Zam
Member
 
Location: Oxford

Join Date: Apr 2010
Posts: 51
Default

This paper (apologies for self-publicity) might be of interest to you:

High-throughput microbial population genomics using the Cortex variation assembler. Z Iqbal, I Turner, G McVean, Bioinformatics 2012
http://bioinformatics.oxfordjournals....full.pdf+html

Fast highly specific and sensitive variant discovery and genotyping for microbial genomes by multi-sample de novo assembly, allowing direct genome comparison without using a reference. The paper contains examples with short and long reads, reproducing results from published papers (drug-resistant/susceptible S.aureus, in-host evolution) with a fraction of the time and effort.
Zam is offline   Reply With Quote
Old 12-19-2012, 11:44 PM   #8
ekg
Member
 
Location: Boston, MA

Join Date: Apr 2010
Posts: 36
Default

@Zam, this looks great! I'm looking forward to reading it.
ekg is offline   Reply With Quote
Old 02-27-2014, 01:09 PM   #9
kjm
Junior Member
 
Location: North Carolina

Join Date: Jul 2013
Posts: 4
Default

Hi,

I've been reading through a number of threads and manuals about this and still don't understand how to identify if SNPs are homozygous or heterozygous when viewing the vcf file. Has anyone seen any issues with using the -s parameter in the bcftools view for setting the ploidy as mentioned above? Not sure if it matters or not, but my samples aren't bacteria, but are haploid males. Thanks!

Oh and yes, we are looking into Freebayes too.
kjm is offline   Reply With Quote
Old 02-28-2014, 12:32 PM   #10
kjm
Junior Member
 
Location: North Carolina

Join Date: Jul 2013
Posts: 4
Default

Hi again. So I decided to do a test a merged sample (merged from three runs of the same sample) and got an error I'm not sure about.

Here is the command I used: samtools mpileup –Dsuf ref.fa Sample.bam | bcftools view –bcgv –s Sample.txt - > Sample.raw.bcf

Then got this:
[mpileup] 1 samples in 1 input file
<mpileup> set max per-file depth to 8000
<bcf_hdr_subsam> 1 samples in the list but not in BCF.

So it's something wrong with the text file I used for the -s parameter in bcftools for ploidy(I think?). The text file seems like it would be pretty straight forward, but I don't know what could cause the issue for the "samples in list but not in BCF" error. Any ideas? Thanks.
kjm is offline   Reply With Quote
Old 05-11-2014, 01:26 PM   #11
ZLounsberry
Junior Member
 
Location: Davis, CA

Join Date: Mar 2014
Posts: 2
Default

Hey kjm,
Assuming you haven't answered this on your own yet (or assuming someone else is having a similar error and needs the thread to not die at your question), I think I have an answer.

If you are piping it the way you are in your post, making the "Sample.txt" file simply a dash, tab, and 1 should cover it. That way it reads your file name as the piped filename and continues as usual. Got rid of my "<bcf_hdr_subsam> 1 samples in the list but not in BCF." error anyway... Hope that helps someone.

Sample.txt should read:
- 1
ZLounsberry is offline   Reply With Quote
Old 01-13-2015, 12:28 AM   #12
polpol
Junior Member
 
Location: Israel

Join Date: Jan 2015
Posts: 3
Default variant calling in SAMtools

Dear ZLounsberry,

I saw your answer , but I didn't understand how the text file should be .

I am running the following command:
samtools mpileup –uf ucsc.hg19.fasta child.bam father.bam mother.bam | bcftools view -bvcgT trioauto -s family.txt - > femvar.vcf &

My txt file is:

child.bam
father.bam
mother.bam


I get an error:
[mpileup] 3 samples in 3 input files
<mpileup> Set max per-file depth to 2666
<bcf_hdr_subsam> 3 samples in the list but not in BCF.


Do you have any suggestion how to fix the problem?
Thank you for your help,
polpol is offline   Reply With Quote
Old 01-13-2015, 07:14 AM   #13
ZLounsberry
Junior Member
 
Location: Davis, CA

Join Date: Mar 2014
Posts: 2
Default

Hello polpol,
Try changing your "family.txt" file to read:

- 1

instead of:

child.bam
father.bam
mother.bam

If that does not work, try:

- 3

(so a dash, a tab character, and a 1 or, if that doesn't work, maybe a 3). It's a bit off topic for this thread, so feel free to email me (my username AT gmail DOT com) and let me know if you need a hand.
ZLounsberry 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 11:38 PM.


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