SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
samtools pileup BioSlayer Bioinformatics 9 01-16-2012 03:10 AM
samtools mpileup and bcftools command lines to compare individual differences? aslihan Bioinformatics 1 11-17-2011 02:20 PM
samtools pileup raghu1990 Bioinformatics 0 06-28-2011 02:26 AM
samtools pileup wisosonic Bioinformatics 0 05-11-2011 04:57 AM
SAMtools command line ??? Pawan Noel Bioinformatics 6 11-16-2010 10:42 AM

Reply
 
Thread Tools
Old 08-19-2009, 06:45 AM   #1
kellyv
Junior Member
 
Location: Penn State

Join Date: Jul 2009
Posts: 3
Post SAMtools pileup command

I am having trouble with the samtools pileup command and would appreciate any advice. I'm probably doing something wrong, but I have read the documentation and can't figure out what. Originally I was starting with a sam file produced by BWA. But that wasn't working so I switched to a sam file generated by using bowtie2sam.pl on a bowtie file.

I used the following command to create the unsorted bam file:
samtools view -buSt chrM.fa -o file.bam file.sam

And then the following to sort the file:
samtools sort file.bam file.bam.sorted

Then I tried a variety of things to produce the pileup file.

samtools pileup -f chrM.fa file.bam.sorted.bam
--> no output, no error

samtools pileup -t chrM.fa.fai file.bam.sorted.bam
samtools pileup -f chrM.fa -t chrM.fa.fai file.bam.sorted.bam
-->
[sam_header_read2] 1 sequences loaded.
[sam_read1] reference '$(((((((((((((((((((((((((NMCX1CMDZ7N0N27' is recognized as '*'.
[sam_read1] reference '(((((((((((((((((((((((((((((NMCX1CMDZ7N0N27' is recognized as '*'.
Parse error at line 2: invalid CIGAR character
Abort trap

samtools pileup -f chrM.fa -S file.sam
-->
[samopen] no @SQ lines in the header.
[sam_read1] missing header? Abort!

samtools pileup -f chrM.fa file.sam
-->
[bam_pileup] fail to read the header: non-exisiting file or wrong format.


The only way I could get pileup output was as follows:
samtools pileup -t chrM.fa.fai file.sam
samtools pileup -t chrM.fa.fai -S file.sam
--> note there is no -S flag in the first one and yet it produced correct output

I want to be able to use the pileup command on the bam file because I am assuming it will be much faster on large files. Can I do that?

And I was told that it's better to use -f than -t--is that true, or does it not matter?
kellyv is offline   Reply With Quote
Old 08-19-2009, 07:19 AM   #2
dawe
Senior Member
 
Location: 4530'25.22"N / 915'53.00"E

Join Date: Apr 2009
Posts: 258
Default

I have same issues here when I use the conversion tools. I mean, if I perform the alignment with bwa (which produces SAM output) everything works fine, like this:

bwa aln genome file.fq > file.sai
bwa samse file.sai | samtools view -bt genome.fai -o results.bam -

samtools sort results.bam sorted && mv sorted.bam results.bam
samtools index results.bam

samtools pileup [-f genome] results.bam

I've tried export2sam.pl for old runs but it has some issues which I believe are related to the fact export files have the chromosome (or genome) file name in the 'chromosome' column. This is reported in SAM from the converter but there's no, say, chr5.fa chromosome in the reference used for samtools. The worst happens when I convert eland_export files aligned against a multifasta reference: the expected column for chromosome contains only, say, genome.fa... the export2sam.pl places all the reads as they were aligned to chromosome 'genome.fa':

11II-08F1-GA_1:6:1:0:1237 81 genome.fa 138895 275 36M = 138757 -174 CGAATTAGTAGATTGTTTGTAAAAATCGGTATTGTN <:<997<;<:<<::;9999;;;:;9<<:
<<;:;;/& MD:Z:C35 SM:i:133
11II-08F1-GA_1:6:1:0:1237 161 genome.fa 138757 275 36M = 138895 174 ACTTGTTTTAAATAAACAGATTTTCCACTCATGTTG @@@@?=BBBCBB@?@BAB=5>@B@B@A@
AA@@?<@@ MD:Z:36 SM:i:142
11II-08F1-GA_1:6:1:0:399 97 genome.fa 179127 275 36M = 179285 194 NATATCATAACCATAAAAAAAAAAGCACAATTCAAC &1<<<<<<<<<<<<<<<<<<<<<8579:
:<<:::98 MD:Z:A35 SM:i:133

I don't have bowtie ouputs here at the moment, you may check if the bowtie chromosome column is the name of the file or the chromosome...

HTH

Davide
dawe is offline   Reply With Quote
Old 08-20-2009, 04:48 PM   #3
ech
Junior Member
 
Location: USA

Join Date: Jul 2009
Posts: 7
Default * character in pileup output

Here is an example of line in pileup output which contains * characters
chrIII 8884091 T W 95 126 60 12 ,,*****,*aA^]A aab`a^aHa^aa

I do not believe this character is documented. Would be nice to hear back what it means.
ech is offline   Reply With Quote
Old 08-20-2009, 09:05 PM   #4
baohua100
Senior Member
 
Location: Canada

Join Date: Jun 2008
Posts: 103
Default

How can I pileup all the sites of reference sequence , but not only the sites covered by at least 1 ?
baohua100 is offline   Reply With Quote
Old 08-21-2009, 08:09 AM   #5
nilshomer
Nils Homer
 
nilshomer's Avatar
 
Location: Boston, MA, USA

Join Date: Nov 2008
Posts: 1,285
Default

Quote:
Originally Posted by baohua100 View Post
How can I pileup all the sites of reference sequence , but not only the sites covered by at least 1 ?
You can get a bioinformatician or computer scientist to write a quick perl script to do this for you. We need jobs!
nilshomer is offline   Reply With Quote
Reply

Tags
samtools

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:25 AM.


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