SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Samtools SNP Caller versus MAQ SNP caller BertieWooster Bioinformatics 5 05-18-2012 01:41 PM
NGS whole genome sequence versus sequence capture for quality control houkto Bioinformatics 0 02-02-2012 04:16 AM
Velvet Assembler: expected coverage versus estimated coverage versus effective covera DMCH Bioinformatics 1 11-30-2011 04:21 AM
Transcriptome assembler versus genome assemblers samanta General 9 08-02-2011 09:21 AM
All chromosomes and maq Layla Bioinformatics 1 02-13-2009 08:36 AM

Reply
 
Thread Tools
Old 04-10-2012, 04:17 PM   #1
ramirob
Member
 
Location: Vermont

Join Date: Apr 2012
Posts: 14
Default Samtools on Chromosomes Versus Genome

Hi everyone,

Does anybody have experience running samtools on CASAVA alignments?

I am trying to call SNPs using samtools on a CASAVA/illumina alignment.

CASAVA outputs an alignment of the type name_export.txt.gz

and it automatically converts to BAM format when using their SNP caller.

Now I want to call SNPs using Samtools.

The alignment was done using each chromosome individually:

genomes/Homo_sapiens/UCSC/hg18/Sequence/Chromosomes/Chr...

However, the BAM alignment provided by CASAVA is for the whole genome, so then should I use the whole reference genome?:

samtools mpileup -f ../genomes/Homo_sapiens/UCSC/hg18/Sequence/WholeGenomeFasta/genome.fa sorted.bam > sorted.mpileupOutput

the problem is that, in that case the mpileup output doesn't seem to recognize it as I get:

bash-3.2$ head -n 100 sorted.mpileupOutput
chr1.fa 3309 N 1 ^NT C
chr1.fa 3310 N 1 G C
chr1.fa 3311 N 1 C F
chr1.fa 3312 N 1 C F
chr1.fa 3313 N 1 C F


Any ideas??

Thanks in advance,

Ramiro
ramirob is offline   Reply With Quote
Old 04-10-2012, 04:33 PM   #2
sdriscoll
I like code
 
Location: San Diego, CA, USA

Join Date: Sep 2009
Posts: 438
Default

try this. mpileup doesn't do too much. bcftools is the program that actually does the SNP calling.

Code:
samtools mpileup -uf ../genomes/Homo_sapiens/UCSC/hg18/Sequence/WholeGenomeFasta/genome.fa sorted.bam | bcftools view -bvcg - > sorted.bcf
bcftools view sorted.bcf > sorted.vcf
head sorted.vcf
sdriscoll is offline   Reply With Quote
Old 04-10-2012, 06:38 PM   #3
bioinfosm
Senior Member
 
Location: USA

Join Date: Jan 2008
Posts: 482
Default

or you can try splitting the bam file by chromosome and then calling SNP on them separately.. you could also try starting from the fastq or qseq files, run alignment using something like bwa and then call variants
__________________
--
bioinfosm
bioinfosm is offline   Reply With Quote
Old 04-10-2012, 08:11 PM   #4
sdriscoll
I like code
 
Location: San Diego, CA, USA

Join Date: Sep 2009
Posts: 438
Default

indeed...i align data using BWA when i'm looking for SNPs.
sdriscoll is offline   Reply With Quote
Old 04-11-2012, 05:22 PM   #5
AllanLindh
Junior Member
 
Location: Santa Cruz, CA, USA, Earth

Join Date: Apr 2012
Posts: 5
Default How to combine .vcf files into whole genome

This discussion is very useful, but when I run mpileup plus bcftools to produce .vcf files for individual Chrm, I have found no way to combine them into whole genome .vcf file for further processing. vcf-concat seems like it might do it, but it never runs due to diverse errors. Most recent is:
****************************
The use of -1 for unknown number of values is deprecated, please use '.' instead.
FORMAT=<ID=PL,Number=-1,Type=Integer,Description="List of Phred-scaled genotype likelihoods, number of values is (#ALT+1)*(#ALT+2)/2">
****************************
I am running most recent version of samtools (0.1.18) -- has anyone else noticed a mismatch in the FORMAT specs between samtools and vcftools??

Any suggestions at all would be very welcom, this is driving me nuts.

Thanks
Allan Lindh
AllanLindh is offline   Reply With Quote
Old 04-13-2012, 11:01 AM   #6
ramirob
Member
 
Location: Vermont

Join Date: Apr 2012
Posts: 14
Default

sdriscoll,

It seems that mpileup looking like:

chr1.fa 3309 N 1 ^NT C
chr1.fa 3310 N 1 G C
chr1.fa 3311 N 1 C F
chr1.fa 3312 N 1 C F

is a bad sign, it looks to me that it is not mapping to the ref. genome properly. I did as you suggested and everything was a SNP.
ramirob is offline   Reply With Quote
Old 04-13-2012, 11:14 AM   #7
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

Double check that the chromosome names in your .bam file match the chromosome names in your fasta. And I'd remake the fa.fai with samtools faidx, see if that completes without an error.
swbarnes2 is offline   Reply With Quote
Old 05-19-2012, 02:19 PM   #8
frc1230
Member
 
Location: United States

Join Date: Jul 2011
Posts: 10
Default

I am currently trying to deal with exactly this problem: CASAVA alignments that choke mpileup. We use mpileup output directly for analysis with our own SNP finders, and it works very well with output from a BWA-based aligner (stampy). But with CASAVA output, as stated in another response, mpileup thinks everything is a SNP, all reads with phred scores of 0. This is not the case, fortunately, as the same read data run through stampy produces an alignment file that behaves well with mpileup.
frc1230 is offline   Reply With Quote
Old 05-22-2012, 03:38 AM   #9
ramirob
Member
 
Location: Vermont

Join Date: Apr 2012
Posts: 14
Default

I must admit I just gave up on Casava alignments, just use BWA and I am exploring other aligners.
ramirob is offline   Reply With Quote
Old 05-22-2012, 08:19 AM   #10
frc1230
Member
 
Location: United States

Join Date: Jul 2011
Posts: 10
Default

I like stampy-BWA, but it does take a lot of cpu time
frc1230 is offline   Reply With Quote
Reply

Tags
samtools casava mpileup

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 05:00 PM.


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