SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
SOAP alignment format convert to SAM/BAM KevinLam Bioinformatics 31 01-10-2018 08:05 PM
SAM/BAM format to wiggle format pinki999 Bioinformatics 19 08-12-2015 12:35 AM
SAM to CUFFLINKS SAM format repinementer Bioinformatics 4 03-15-2012 08:53 AM
Looking process to convert gff3 format into ace format or sam format andylai Bioinformatics 1 05-17-2011 02:09 AM
anyone help me on bowtie format -> sam format! tninja Bioinformatics 2 04-25-2010 09:33 PM

Reply
 
Thread Tools
Old 03-10-2009, 02:36 PM   #21
lh3
Senior Member
 
Location: Boston

Join Date: Feb 2008
Posts: 693
Default

You can generate ref_list file by running "bam faidx" on your reference sequence. The index file can be used with import. Note that faidx allows you to quickly extract subsequence from the genome, which may be useful to you.
lh3 is offline   Reply With Quote
Old 03-11-2009, 03:19 PM   #22
lparsons
Member
 
Location: NJ

Join Date: Nov 2008
Posts: 28
Default

Thanks for the quick replies. I've tried with various values, etc. and the indexing step still seg faults. Any other ideas on debugging this? Perhaps using an older version? Thanks for any ideas.
lparsons is offline   Reply With Quote
Old 03-17-2009, 10:55 PM   #23
thondeboer
Member
 
Location: Bay Area

Join Date: Jan 2009
Posts: 24
Default

I have been looking at the SAM format to see if it is something we should consider for the output of the mapping for a genome assembly we are doing at Complete Genomics. As you may know, we have quite an unusual read structure that may be difficult to represent in SAM (5+10+10+10 times two, mate paired reads). The problem lies in the fact that the there is some overlap between the first 5 base read and the second 10 bases (we call them negative gaps).

The read could be

acgtc tcgattgcgg ...

which maps to the reference like this

acgTCgattgcc...

The capital TC show the negative gaps where there were actual overlaps in the read sequence.

Anyone now how we could represent this in SAM? Can the CIGAR standard deal with negative numbers? 5M-2M8M ?)

Should we map the 5 base read and the other 30 bp read as two separate reads? But we also have mate pairs, so how should we represent those? And would any other tool be able to deal with a read structure like this?

Thanks,

Thon
Complete Genomics
thondeboer is offline   Reply With Quote
Old 03-18-2009, 01:31 AM   #24
lh3
Senior Member
 
Location: Boston

Join Date: Feb 2008
Posts: 693
Default

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?
lh3 is offline   Reply With Quote
Old 03-18-2009, 09:13 AM   #25
thondeboer
Member
 
Location: Bay Area

Join Date: Jan 2009
Posts: 24
Default

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
thondeboer is offline   Reply With Quote
Old 03-18-2009, 02:20 PM   #26
lh3
Senior Member
 
Location: Boston

Join Date: Feb 2008
Posts: 693
Default

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?
lh3 is offline   Reply With Quote
Old 03-20-2009, 01:39 AM   #27
lh3
Senior Member
 
Location: Boston

Join Date: Feb 2008
Posts: 693
Default

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.

Quote:
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 is offline   Reply With Quote
Old 03-22-2009, 06:51 PM   #28
nilshomer
Nils Homer
 
nilshomer's Avatar
 
Location: Boston, MA, USA

Join Date: Nov 2008
Posts: 1,285
Default

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
nilshomer is offline   Reply With Quote
Old 03-25-2009, 11:43 AM   #29
qtrinh
Member
 
Location: Canada

Join Date: May 2008
Posts: 20
Default

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
qtrinh is offline   Reply With Quote
Old 03-25-2009, 11:48 AM   #30
lh3
Senior Member
 
Location: Boston

Join Date: Feb 2008
Posts: 693
Default

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.
lh3 is offline   Reply With Quote
Old 03-25-2009, 11:58 AM   #31
qtrinh
Member
 
Location: Canada

Join Date: May 2008
Posts: 20
Default

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
qtrinh is offline   Reply With Quote
Old 04-08-2009, 09:49 AM   #32
thondeboer
Member
 
Location: Bay Area

Join Date: Jan 2009
Posts: 24
Default

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...
thondeboer is offline   Reply With Quote
Old 04-08-2009, 10:23 AM   #33
apfejes
Senior Member
 
Location: Oakland, California

Join Date: Feb 2008
Posts: 236
Default

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?
__________________
The more you know, the more you know you don't know. —Aristotle
apfejes is offline   Reply With Quote
Old 04-09-2009, 03:28 AM   #34
lh3
Senior Member
 
Location: Boston

Join Date: Feb 2008
Posts: 693
Default

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.
lh3 is offline   Reply With Quote
Old 04-15-2009, 03:29 PM   #35
TylerBackman
Member
 
Location: Riverside, CA

Join Date: Oct 2008
Posts: 13
Default

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!
TylerBackman is offline   Reply With Quote
Old 04-17-2009, 02:05 PM   #36
gcrdb
Junior Member
 
Location: usa

Join Date: Jan 2009
Posts: 9
Smile 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
gcrdb is offline   Reply With Quote
Old 04-18-2009, 01:07 AM   #37
lh3
Senior Member
 
Location: Boston

Join Date: Feb 2008
Posts: 693
Default

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.
lh3 is offline   Reply With Quote
Old 04-20-2009, 10:45 AM   #38
gcrdb
Junior Member
 
Location: usa

Join Date: Jan 2009
Posts: 9
Default

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
gcrdb is offline   Reply With Quote
Old 04-20-2009, 11:03 AM   #39
lh3
Senior Member
 
Location: Boston

Join Date: Feb 2008
Posts: 693
Default

should be: ../samtools pileup -t ex1.fa.fai ex1.sam or ../samtools pileup ex1.bam. I have added a Makefile.

Note that there is a companion format called BAM which is the binary representation of SAM. Most of samtools commands work on BAM only. I know having two formats is a bit confusing, but this is necessary for faster parsing.
lh3 is offline   Reply With Quote
Old 04-20-2009, 03:46 PM   #40
gcrdb
Junior Member
 
Location: usa

Join Date: Jan 2009
Posts: 9
Default

lh3,
Thanks for quick response, pileup is working now!
Here is some questions about pileup format when I look at them at first time:
(1) what is a "*" in read bases , which is not documented in "http://samtools.sourceforge.net/pileup.shtml".
(2) Is it okay for a base in the same position pile-up twice ? (chr1 1949878 occur twice in my first pileup output)
thanks,
Below is the piece pile-up output I found the problem:
chr1 1949878 A A 142 0 60 55 C$....,,,...,,.C..,.,,..,...,,,.,.,.+1C,,,,..,.,........,^F
,^], &5,2IIII5I<II+%5=I+II(II8@CII3*I0I+IIII,I$@IIAAIIDI@I*5
chr1 1949878 * */+C 38 38 * +C 13 4 30 8
chr1 1949879 A A 150 0 60 54 ....,G,...,,.-1G..-1G.,.,,..,...,,,.,.,.,,,,..,.,........,,
, II*IIII;I.II&(&1I3II%9III:II7&I&@&IIII5I"6IIIIG.I33I$6
chr1 1949879 * -G/* 481 481 -G * 6 9 30 9
chr1 1949880 G G 25 0 60 55 .$A$..,A,..A,,*.*A,A,,A.,A..,,,.,A,.,,,,.+1A.,.,........A,,
^], +(,.III)8-II8%D.I0II#,I@5III.$I+I,IIII$I2EIIIIIIIIIE&?/
chr1 1949880 * */+A 350 350 * +A 24 3 19 9
chr1 1949881 A A 162 0 60 53 ..,,,...,,....,.,,..,...,,,.,.,.,,,,..,.,........G,,, 6DI
II3I$II81D)I%II+.III'III$I2I+IIII+II<II46IIAHI.2IB
gcrdb 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:14 PM.


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