[thondeboer]

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

[lparsons]

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.

[lh3]

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.

[lparsons]

Yes, I sorted it just fine. In fact the indexing step will complain that the file isn't sorted.

One issue could be that I just realized that the ref_list file I gave during the import didn't have the reference size in it. I assume this means the length of the reference sequence? I'll have to give that a try when I first import the file (convert to bam).

[lh3]

Have you sorted the alignment first? Indexing in 0.1.2 has a bug, but should not cause segfault. Thanks.

[lparsons]

I imported an ELAND alignment and was able to convert into SAM, then to BAM, then sort it. However, at the indexing step I too ran into a segmentation fault. I'm using the 0.1.2 version from the download page, not from SVN. Any suggestions?

[webbrewer]

Have you sorted the alignment? Indexing only works for sorted alignment. Also remember to use the latest bwa. The old version may generate some funny alignments, though this happens very rarely.

I hadn't sorted it before. Now I ran "samtools sort", then "samtools index" on the sorted output. It resulted the same with seg fault.
I am using bwa version 0.4.5. Is there a newer svn version?
samtools index works without issue on converted MAQ alignments.

[lh3]

Have you sorted the alignment? Indexing only works for sorted alignment. Also remember to use the latest bwa. The old version may generate some funny alignments, though this happens very rarely.

[myrna]

samtools index

Have you tried samtools view foo.bam?

If you get the sam alignments back, then all should be well. I believe you get a warning if the .bam file is unsorted, but perhaps you should try this if you haven't already:

samtools sort foo.bam bar.sort

[webbrewer]

samtools import is for making a .bam file from a .sam file. Why are you attempting to run this command on a .bam file?

Oops. I meant to say that "samtools index" seg faults.

[myrna]

samtools import

samtools import is for making a .bam file from a .sam file. Why are you attempting to run this command on a .bam file?

[webbrewer]

samtools index seg fault

I am using the most current version of samtools from svn.
I successfully ran the "samtools import" command on my .sam file from bwa.
When I then run "samtools index" on the .bam file, it seg faults.
Let me know if you need more information to determine what is causing this.

[lh3]

maq2sam-short is for the .map files generated by maq-0.6.x, while maq2sam-long for files generated by maq-0.7.x. Sorry for the confusion, and one of the aims of SAM is to avoid such confusions in future.

[kon104]

What's the difference between maq2sam-short and -long?

Also, short seems to segfault on 64-bit versions of Red Hat and Ubuntu... Am I missing something?

[lh3]

I am planning to release samtools-0.1.2 which fixed some bugs in the old version and added new features. For now you can check out source codes from SVN. It should be quite close to 0.1.2.

The new version comes with a Bayesian indel caller, although it is just a prototype at present. The strength of the samtools' caller is that it makes use of reads mapped without indel. Using this information helps to reduce false negatives. In addition, the new caller gives genotype rather than just saying there is an indel. You cannot easily tell from maq's indelpe if the indel is a heterozygote or a homozygote. With the new caller, the filters could be: a) the indel score; b) two indels should not be too close to each other.