CIGAR and Sequence length incosistent

[bnfoguy]

Hello,

I am trying to convert a .sam file into .bam file and I get the following error:

CIGAR and Sequence length are inconsistent. Below is the offending line:

ERR035488.960541 147 chr1 27205782 99 72M = 27205550 -304 TATGCCTCCTTGAGTGTCAGTGGCGTGATCTTGGCCCGGCTCACACCGGCCGGCAGGAAGTCTAGTAGGCAG 8<=21<?7A<<DA8>>@CEDE.D<DDB?@;BAC6;48EB@AE4>4@&FEB6FEFGFD1GFAC@DFAFBCBB> PQ:i:46 SM:i:96 UQ:i:0 MQ:i:96 XQ:i:15 NM:i:0

I need to convert this into a .bam file and further into a .flt.vcf format for indel calling. So far, I have cut the first 1 million alignments from the original sam output and have gotten as far as creating an {prefix}_sorted.bam file.

[maubp]

I don't see what's wrong with that - the sequence and the quality string appear to both be 72 characters, and the CIGAR string of 72M is consistent with that (72 matches or mismatches). Are you sure you are looking at the right line number?

[Richard Finney]

What software produced the sam file?
What version of samtools are you using ?

Please do this:
check with this command:

cat yoursample.sam | awk '{if (length($10)!=length($11)) print p+1" "$0;p++}'

It will print out the line number and offending line.

Line number includes sam header, not just sequence read lines.

[katussa10]

I have the same problem:
Line 18632882, sequence length 48 vs 74 from CIGAR
Parse error at line 18632882: CIGAR and sequence length are inconsistent



I used:
"cat yoursample.sam | awk '{if (length($10)!=length($11)) print p+1" "$0;p++}'"

and got the following:

76327 @PG ID:Bowtie VN:0.12.7 CL:"bowtie -p 8 --sam -e 120 Trinity_name -1 ../PAIR_END_1_Trim_the_reads_by_quality.fastqsanger -2 ../PAIR_END_2_Trim_the_reads_by_quality.fastqsanger Mysample.sam"
18632882 SOLEXA_0019_FC:1:27:11595:10519#AACTAC 163 c0_seq1 4639 255 74M = 4646 81 ATGAAATCGAACGGTTTTCTACACGATCGTGCAAGTGAAACACATGGG

However, when I used:
head -n 18632882 Mysample.sam |tail -n 1



SOLEXA_0019_FC:1:27:11595:10519#AACTAC 163 1089-1104_All_comp1980_c0_seq1 4639 255 74M = 4646 81 ATGAAATCGAACGGTTTTCTACACGATCGTGCAAGTGAAACACATGGGTTATTGAGAGCATTGGGAGTCATTAG IIIIIIIIIIIIIIIIIIIIIIIIIIIIDFDGGGGEHIIIGIIIIIIIHIIIIEIIIIIHIIFGIGIGIIIIIH XA:i:0 MD:Z:74 NM:i:0


What should I do next?

[Richard Finney]

Well, there's "od" : "octal dump".

Try

head -n 18632882 Mysample.sam |tail -n 1 | od -c

(or try "od -x" , whichever you can read easier )

See if there's some strange control character snuck in by some piece of software.

I'd just cut out the read and move on, like this ..

cat yoursample.sam | awk '{if (p!=18632881) print $0;p++}'" > fixed.sam

# note the zerobased couting using 18632881, not 18632882

Nobody's going to miss one read.

[mxr1895]

I'm dealing with the same issue. More specifically, I am trying to convert a sam file generated using 'miraconvert' from a MIRA *.caf file.
This issue has been raised in [mira_talk]:

http://www.freelists.org/post/mira_t...sam-conversion

Note that the above code: "cat yoursample.sam | awk '{if (length($10)!=length($11)) print p+1" "$0;p++}'"

will not necessarily work for this issue, since in the example above (and in my case) both fields are of the same length.
anybody know of a quick way of reading/interpreting the sum of a CIGAR field in order to remove the offending lines?
Cheers,
Miguel

[EricHaugen]

anybody know of a quick way of reading/interpreting the sum of a CIGAR field in order to remove the offending lines?

Here's a little Perl script I've used for that purpose:

#!/usr/bin/perl -w
use strict;

# Pass through valid-looking SAM records
# Suppress records with invalid CIGAR strings, sending warning to STDERR

while (<>) {
    my ($readName,$flag,$refName,$refOffset,$mapQuality,$cigar,$ignore1,$ignore2,$ignore3,$seq,$qual,$tag,@other) = split /\t/;
    unless (defined($seq)) {
        warn "No sequence found in $_\n";
        next;
    }
    my $orig = $cigar;
    my $readbases = 0;
    while (length($cigar) > 0) {
        if ($cigar =~ /^([\d]+)([MISP\=X])(.*)/) {
            $readbases += int($1);
            $cigar = $3;
        } elsif ($cigar =~ /^([\d]+)([DNH])(.*)/) {
            # ignore
            $cigar = $3;
        } else {
            warn "Failed to parse ($orig) in $_\n";
            last;
        }
    }
    my $seqlen = length($seq); 
    if ($seqlen != $readbases) {
        warn "$seqlen bases in $seq but $readbases bases parsed from CIGAR string $orig, $_";
    } else {
        print;
    }
}

[maubp]

Bastien hasn't yet fixed this bug as of MIRA 4.0 RC5,
http://www.freelists.org/post/mira_t...m-conversion,3

[cankut]

Hi,

Here you can find information about why I got the same problem and how I solved it.

Today I had the same problem... when I was trying convert sam file to bam file with following command:
# samtools view -bS myfile.sam -o myfile.bam

To solve my problem first I tried
# cat yoursample.sam | awk '{if (p!=18632881) print $0;p++}' > fixed.sam
as Richard Finney mentioned.
But it is not the best solution when the error continues within different lines.

So I checked my .sam file.
# head -n error_line myfile.sam | tail -n1
Here everything was ok.
Then I copied the mapped sequence and search it on my readfile.fa
I noticed that the following line had "--" instead of seq data.
example:
>seq1
GATATTGGCGCGGCTCAATCA
--
>seq2
GATATTGGCGCGGCT
>seq3
GATATTGGCGCGGC

And it was exist several times in different lines.
So I removed these undesirable symbols with command:
# sed -i '/--/d' readfile.fa

Then I started again from aligning the curated readfile with using bwa
....
then I was succesful to convert sam to bam and all following steps done without any error.