[juan]

I was curious if anyone has found a workaround to this error? It is coming up for me when viewing novoalign alignments. I am using samtools 0.1.6 (r453), and get the following error when loading a samtools converted bam file from novoalign's sam file

samtools: bam_lpileup.c:116: tview_func: Assertion `l == tv->n_pre' failed.


thanks

I am also having this problem on linux. Did you find a solution?

[dgacquer]

I am working with the Blast aligner and therefore I have to convert my blast alignment files to sam format. I think I have found a bug in the third party perl script "blast2sam.pl" which is available with samtools on Sourceforge.

The problem comes from the "aln2cm" subroutine used to build the CIGAR for a query. Apparently, the script switches 'D' with 'I'. For instance, if the CIGAR must be "103M1I2D10M", the script outputs the following "103M1D2I10M".

If you look at the blast2sam.pl script, line 88:

"push(@$cigar, $cmaux->[1] . substr("MDI", $cmaux->[0], 1));"

and replace this line with:

"push(@$cigar, $cmaux->[1] . substr("MID", $cmaux->[0], 1));" (to swith 'I' and 'D' in the "MDI" string)

the problem seems to be fixed.

Note that this problem is only visible when you are working with gapped queries, which leads to a "CIGAR length and query length are inconsistent" error message, and that you have to use the -s option with blast2sam.pl so that the query string can be included in the resulting sam file.

I think a user named corthay made a post some time ago about a similar problem. If other people have faced the same problem, please let me know.

I hope this helps

Best regards

David

[Alexey]

eference '163285701' is recognized as '*', sequence and quality are inconsistent

>

@luisczul

>Samtools indicates that the error happens to line 164507. What does >that line look like?As for the second problem, it seems like a bug. You >are using very short reference. Could you send me an example? >Thanks.

>@henry
>Are you generating results with "samse -n"? With -n, the output is NOT >sam. You can tell this from the bwa manual page.

>These three questions are less relevant to samtools. They are mostly related to bwa.

I am having similar problem with converting sam to bam:

[root@master maq_bwa_run]# samtools view -bt hg19.fa.fai -o aln.bam aln.sam
[samopen] SAM header is present: 93 sequences.
[sam_read1] reference '163285701' is recognized as '*'.
Parse warning at line 29354: mapped sequence without CIGAR
Parse error at line 29354: sequence and quality are inconsistent
Aborted

What could happen?

the line 29354 in aln.sam looks like this

0 chr5 163285701 0 0M * 0 0 * XT:A:R NM:i:0 X0:i:100307042 XM:i:0 XO:i:0 XG:i:0 MD:Z:0
~


the lines 29352, 29353, 29354 and 29355 look like this:


ran_melanocytes:853_1092_1366 4 * 0 0 * * 0 0 TTCGCGGGGGATAGACTGTACAATCTAGGCCGGTGGTGTATGGCCCTGT AA>>AA@=>A>A;><>>>;>?<:<:7?<='6<(,))<6/&38&&(531)
ran_melanocytes:853_1093_11 4 * 0 0 * * 0 0 NANTNGNAGGCGNNANNNTNGNCNGACNNNNNNGNGGANANTNTNNNGA -"1-"=-"/-"):8,:-"-"5-"-"-"<-"1-"1-";,,-"-"-"-"-"
0 chr5 163285701 0 0M * 0 0 * XT:A:R NM:i:0 X0:i:100307042 XM:i:0 XO:i:0 XG:i:0 MD:Z:0
ran_melanocytes:853_1093_36 4 * 0 0 * * 0 0 TATAAATTTCGAGGAAATTAGACATATCTCCGTCGTAGGGATCCCCTGG =3:=;=<<<=@/=?7:7;8?><?15A57:?=;;6:=95?:8,,,,,/
~



Thanks
Alexey

[bgulko]

Parse error at line 1: invalid CIGAR operation

RE: Posts 72 & 74 .

hi, I have trouble conveting sam to bam.. I tried both:

samtools import ref .fai in.sam out.bam
got error:
[sam_header_read2] 22 sequences loaded.
[sam_read1] reference '-143963499' is recognized as '*'.
Parse error at line 1: invalid CIGAR operation
Aborted

samtools view -bt ref .fai -o in.sam out.bam
and got similar error:
[sam_header_read2] 22 sequences loaded.
[sam_read1] reference '' is recognized as '*'.
[main_samview] truncated file.

thanks,

While this may have been answered in the interim, I didn't see one, so I thought I'd post my own.

This happens when you try to use the

bwa samse

command with the -n argument to generate a .sam file. With the -n argument, this command lists multiple matches per read. The resulting output file is a list of matches and NOT a .sam file. If that output file is treated as a .SAM file and fed into

samtools import

to create a BAM file, then samtools will generate the sort of error output you described. If you remove the "-n" argument to "bwa samse" then the output will be a proper SAM file, but with only one alignment per read.

I ran into this myself, and was only able to figure it out when I ran bwa in a debugger and traced the behavior.

--bg

[Fred13]

Am I doing something wrong ??



I used MAQGene to align sequence fragments of a mutant against the reference sequence. I want to use SAMTools to see the alignment of all the fragments against the ref. But I have a problem with the sam file ...

[root@localhost MaqToSam]# /home/.../SAMtools_svn/samtools/misc/maq2sam-long
Version: r439
Usage: maq2sam <in.map> [<readGroup>]
[root@localhost MaqToSam]# /home/.../SAMtools_svn/samtools/misc/maq2sam-long fr6_34.map > fr6_34.sam
[root@localhost MaqToSam]# ../samtools faidx c_elegans.WS200.dna.fa
[root@localhost MaqToSam]# ../samtools view -bS -T c_elegans.WS200.dna.fa -o fr6_34.bam fr6_34.sam
[COLOR="Red"][main_samview] fail to open file for reading.[/COLOR]
[root@localhost MaqToSam]# ../samtools view fr6_34.sam
[COLOR="Red"][bam_header_read] EOF marker is absent.
[main_samview] fail to read the header.[/COLOR]
[root@localhost MaqToSam]# ../samtools import c_elegans.WS200.dna.fa fr6_34.sam fr6_34.bam
[COLOR="Red"][main_samview] fail to open file for reading.[/COLOR]

The problem is that I do it once time and I could see the alignment with tview (I would like to see it in an other format : http://seqanswers.com/forums/showthread.php?t=3249) and now I cannot reiterate this ... what could be the problem ??

[bgulko]

BWA/SAMTOOLS SAM-&gt; BAM - invalid CIGAR operation.

I am encountering exactly the same problem as the post from around #72 in this thread,

... APIs need to start from BAM files (-bam) , not SAM files(no "-sam" at all). I only have SAM files which from bwa, all I need is to convert SAM to BAM.

got error:
[sam_header_read2] 22 sequences loaded.
[sam_read1] reference '-143963499' is recognized as '*'.
Parse error at line 1: invalid CIGAR operation
Aborted

thanks,

Here is a sample of my .SAM file, produced by bwa

>SRR003966.14 1 1
chr14   -37395245       0
>SRR003966.49 1 1
chr12   -25942795       0
>SRR003966.85 5 14
chr14   +102382726      1
chr3    -24756032       2
chr2    -47795311       2
chr5    -176639213      2
chr3    -149933203      2

This SAM file seems a little thin (according to the SAM file format specification I think there should be more fields).

Here is the subsequent command line I used to try to convert it to a BAM file, and its accompanying error

-sh-3.1$ ~/bin/samtools import hg18.fa.gz.fai SRR003966-05.sam SRR003966-05.bam
[sam_header_read2] 49 sequences loaded.
[sam_read1] reference '-37395245' is recognized as '*'.
Parse error at line 1: invalid CIGAR operation
Aborted

I am running BWA 0.5.5 and SAMTOOLS v 0.1.7.

I apologize in advance if this is a newbie question, but I am generally following Heng Li's process as outlined in Short-read Alignment with MAQ and BWA

Any feedback or pointer would be most welcome.

-- bg

[lh3]

SAM/SAMtools FAQ page.

[Zigster]

is there a package available to directly analyze BAM files in R?

[zlu]

I'm wondering for a genome assembly project, will duplicate removal (with samtolls rmdup) and flitering out low mapping quality (e.g mapQ <10) improve my assembly? What do people usually do after mapping with e.g bwa for QC purpose?

Another slightly different issue. Does it matter if 2 fatsq files have the exact identical headers (from 2 solexa runs)? How does samtools sort the bam file? Does it take the header IDs into consideration?

Thank you.

[lh3]

http://samtools.sourceforge.net/cns0.shtml

and

http://samtools.sourceforge.net/samtools.shtml

the pileup command.

[suseq]

varFilter out put

Hi,

I have used samtools to analyse variations using varFilter. So I have imported an alignment file from BWA in sam format, have sorted and run:
1. samtools pileup -vcf ...
2. samtools.pl varFilter...| awk '$6>=20' ...

It did run but I have problems to interpret all the columns. What I think is:
column 1: chromosome
column 2: first base coordinate from the ref.
column 3: ref. base
column 4: consensus base
column 5: ???
column 6: mapping quality
column 7: ???
column 8: read depth
column 9: read base column
column 10: ???

Does somebody know which values are in column 5, 7 and 10? I could not find this information.

[nilshomer]

Perhaps I have misunderstood it but isn't right that properly mapped flag (P string) are only used when read pairs are mapped to the same chromosome with correct insert size? I have about 6% of the properly mapped reads with the P string flag that have mate mapped to different chromosome with 0 insert size as shown below. Has anyone seen this before?

EBRI093151:1:90:555:299#0 pPR1 Chr11 107308221 23 36M Chr12 49 0 TATCCTATTCGAAAGTCGCCATGACCGTGGACATGA BCCBCBBACCB?CBA@BBACCBCAAB<6<?BABBBB XT:A:U NM:i:0 SM:i:23 AM:i:23 X0:i:1 X1:i:1 XM:i:0 XO:i:0 XG:i:0 MD:Z:36

EBRI093151:1:90:555:299#0 pPr2 Chr12 49 23 11S7M1D12M6S Chr11 107308221 0 CTACCGCTTGGGTGGTCATGAATGATTAGCACGCCC AB99@B=BBA>ACCBCCBBCBBCBBBBCBCBB@B@A XT:A:M NM:i:4 SM:i:23 AM:i:23 XM:i:3 XO:i:1 XG:i:1 MD:Z:3A3^T2T5C3

I don't see anywhere in the specification how to set the "properly paired" bit. I would guess this is aligner dependent.

[zlu]

Try the latest version of bwa.

Heng,

This was done with BWA 0.5.4.

For resequencing project, does it really matter if the mates are not properly mapped? Can I instead just filter out those reads with low mapping qualities?

Thank you.

[lh3]

Try the latest version of bwa.

[zlu]

properly mapped Flag

Perhaps I have misunderstood it but isn't right that properly mapped flag (P string) are only used when read pairs are mapped to the same chromosome with correct insert size? I have about 6% of the properly mapped reads with the P string flag that have mate mapped to different chromosome with 0 insert size as shown below. Has anyone seen this before?

EBRI093151:1:90:555:299#0 pPR1 Chr11 107308221 23 36M Chr12 49 0 TATCCTATTCGAAAGTCGCCATGACCGTGGACATGA BCCBCBBACCB?CBA@BBACCBCAAB<6<?BABBBB XT:A:U NM:i:0 SM:i:23 AM:i:23 X0:i:1 X1:i:1 XM:i:0 XO:i:0 XG:i:0 MD:Z:36

EBRI093151:1:90:555:299#0 pPr2 Chr12 49 23 11S7M1D12M6S Chr11 107308221 0 CTACCGCTTGGGTGGTCATGAATGATTAGCACGCCC AB99@B=BBA>ACCBCCBBCBBCBBBBCBCBB@B@A XT:A:M NM:i:4 SM:i:23 AM:i:23 XM:i:3 XO:i:1 XG:i:1 MD:Z:3A3^T2T5C3