[gcrdb]

lh3,
Actually, I tried SAMtools before but somehow "pileup" it's not outputing anything so I am thinking go back to use MAQ.
Did I run the program correctly? (use the example files which come with samtools package)
examples> ../samtools pileup -f ex1.fa ex1.sam
--> return nothing
examples> ../samtools pileup ex1.sam
--> return nothing

thanks again!
g

[lh3]

SAMtools has a SNP caller, based on the same code of MAQ. See this page for more information: http://samtools.sourceforge.net/cns0.shtml. What is missing in SAMtools is a SNP filtration script like "maq.pl SNPfilter", but it is easy to write your own at the moment.

SAMtools' indel caller uses a different algorithm. It outperforms MAQ.

[gcrdb]

Can SAMtools convert SAM back to MAQ?

Hi lh3,
I am glad that SAMtools can do maq2sam , but will it be easy to do sam2maq?
The reason I ask is that I want to MAQ to generate SNP and INDELs. BWA can not do that (yet).

thanks
g

[TylerBackman]

This is an excellent, and very exciting idea. With a standard alignment format (SAM), and a standard raw read format (Phred/Sanger fastq) we can drastically reduce the time most of us spend writing our own file format parsers and converters, and eliminate a common source of error in data analysis (incorrect parsing).

It's great to see the bioinformatics community coming together in this way.

To anyone developing alignment tools: Please include support for this format in future versions of your software!

[lh3]

To thondeboer: currently BWA natively generates alignments in the SAM format. BFAST also generates SAM. We also provide converters for SOAP, Bowtie, Export, novoalign and even blast. However, most of these converters are incomplete in that sometimes they cannot convert every information due to the lack of documentation especially for short indels. So far as I know, all aligners generate its own format. SAM is probably the first effort in unifying the alignment format, in particular for alignment for the new sequencing data.

To apfejes: I am not able to comment much on the Java implemention. I know the I/O part is complete and actually does more nice things than the C version of samtools. you may send an email to the mailing list to ask for the documentations.

[apfejes]

On the subject of tools that use SAMtools, I'm very interested in adding in support for my project, based in java. I'm aware that there exists java based tools for reading/writing in this format, but I'm unable to find any documentation on the software. Has anyone come across any information on how to use the Java SAM tools code?

[thondeboer]

Is there a running list somewhere that shows which programs produce (or take as input) SAM/BAM format? It would help us make a final decision in what format we should produce our mapping results for the mapping we do here at Complete Genomics.

The negative gap structure will probably be something that is dealt with, with special tags (GS, GQ and GC) so software that wants to make optimal use of all the data should use those tags, but most of the software should work "out of the box" if they support SAM/BAM...

[qtrinh]

Hi Heng,
I did do "samtools sort" first. Here is what I did:

samtools import homo_sapiens.fasta.fai S1.sam S1.bam

samtools sort S1.bam S1_sorted

samtools pileup -f homo_sapiens.fasta -c S1_sorted.bam

Q

[lh3]

You should run "samtools sort" first. Pileup works on a stream without loading the whole alignment into the memory and therefore the alignments must be sorted on the chromosomal positions.

[qtrinh]

I got seg fault when doing "samtools pileup". The error message is "[bam_pileup_core] the input is not sorted. Abort!", however I do use the sorted input file. Anyone else seen this before?

Q

[nilshomer]

Just wanted to let everyone know that BFAST (http://genome.ucla.edu/bfast) now supports SAM output. Congratulations Heng (I recognize your coding style from reading MAQ's source!), and others for creating a great tool (samtools) and a good start at an open format (SAM).

Nils

[lh3]

I copied my reply to the samtools-devel mailing list here. I think the following strategy is the best for data generated by Complete Genomics.

Now I prefer to store "acgtctcgattgcgg" as:

SEQ=acgtctcgattgcgg
CIGAR=5M2S8M (or 3M2S10M)

where S stands for "soft clipping/skip on the read". I have checked the source code and think the current "samtools pileup" supports internal soft clipping, which means pileup, consensus/indel calling and glf can be applied. Even if samtools does not support internal soft clipping now, I think it should not be hard to change the code.

[lh3]

Hello Thon,

You can join samtools' mailing lists which can be found here:

http://sourceforge.net/mail/?group_id=246254

You may send around your suggestion and see what others respond. For the moment, I think you may save the read as "acgtggattgcc" and add an optional field to indicate that the first tg is actually an overlap between the first 5bp and the rest of read. In this way, you can use samtools without losing any information. Note that samtools will not look into the optional fields, but the information is kept anyway. Would this be enough for your application?

[thondeboer]

Well...The reads are independent and sometimes don't agree and this is something we need to capture so we could not just remove that information.

Is there any other way in SAM (now or future) that would allow us to capture our read structure? We are going to be producing thousands of genomes very soon and I'm sure many of our customers would want to use a format that is comparable to the one used in the 1000 genome project, but we also would want our read structure to be supported in that format...

Is there a way for us to get involved in the design of the SAM standard?

Thon
Complete Genomics

[lh3]

Negtive length in CIGAR would be good, but that is not supported at the moment. Alternatively, you can save the read as acgtggattgcc and write the CIGAR as 11M. In this way, you cannot get the original read sequence, but I guess this is not so important in most cases. What do you think?