SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
How to obtain full length RNA transcript sequence Annibal RNA Sequencing 9 03-14-2012 06:51 AM
Geneid can predict differnet genes when the length of input sequence are different? bio_tt Bioinformatics 0 01-06-2012 04:47 PM
filter sequence by length NicoBxl Bioinformatics 6 09-09-2011 11:00 AM
Roche gsMapper output exon contigs rather than full-length sequence? sulicon Bioinformatics 0 02-28-2011 04:51 PM
samtools view error: CIGAR and sequence length are inconsistent (tophat/bowtie) glacierbird Bioinformatics 2 06-29-2010 01:58 AM

Reply
 
Thread Tools
Old 06-25-2012, 06:58 AM   #1
bnfoguy
Member
 
Location: New Jersey, USA

Join Date: May 2011
Posts: 17
Exclamation CIGAR and Sequence length incosistent

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.
bnfoguy is offline   Reply With Quote
Old 06-25-2012, 08:20 AM   #2
maubp
Peter (Biopython etc)
 
Location: Dundee, Scotland, UK

Join Date: Jul 2009
Posts: 1,542
Default

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?
maubp is offline   Reply With Quote
Old 06-25-2012, 08:51 AM   #3
Richard Finney
Senior Member
 
Location: bethesda

Join Date: Feb 2009
Posts: 700
Default

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.
Richard Finney is offline   Reply With Quote
Old 07-09-2012, 11:11 AM   #4
katussa10
Member
 
Location: Columbus, OH, USA

Join Date: Jun 2010
Posts: 11
Default

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?
katussa10 is offline   Reply With Quote
Old 07-09-2012, 05:19 PM   #5
Richard Finney
Senior Member
 
Location: bethesda

Join Date: Feb 2009
Posts: 700
Default

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.
Richard Finney is offline   Reply With Quote
Old 12-09-2013, 12:51 PM   #6
mxr1895
Junior Member
 
Location: new zealand

Join Date: Feb 2012
Posts: 6
Default

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
mxr1895 is offline   Reply With Quote
Old 12-09-2013, 04:21 PM   #7
EricHaugen
Member
 
Location: Seattle, WA

Join Date: Sep 2009
Posts: 13
Default

Quote:
Originally Posted by mxr1895 View Post
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:
Code:
#!/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;
    }
}
EricHaugen is offline   Reply With Quote
Old 12-10-2013, 05:50 AM   #8
maubp
Peter (Biopython etc)
 
Location: Dundee, Scotland, UK

Join Date: Jul 2009
Posts: 1,542
Default

Bastien hasn't yet fixed this bug as of MIRA 4.0 RC5,
http://www.freelists.org/post/mira_t...m-conversion,3
maubp is offline   Reply With Quote
Old 02-12-2014, 04:33 AM   #9
cankut
Junior Member
 
Location: Valencia

Join Date: Feb 2014
Posts: 1
Default

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.
cankut is offline   Reply With Quote
Reply

Tags
cigar, sam to bam, samtools, sequence length

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 04:01 AM.


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