Seqanswers Leaderboard Ad

Collapse

Announcement

Collapse
No announcement yet.
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • 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.

  • #2
    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?

    Comment


    • #3
      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.

      Comment


      • #4
        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?

        Comment


        • #5
          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.

          Comment


          • #6
            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]:



            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

            Comment


            • #7
              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;
                  }
              }

              Comment


              • #8
                Bastien hasn't yet fixed this bug as of MIRA 4.0 RC5,

                Comment


                • #9
                  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.

                  Comment

                  Latest Articles

                  Collapse

                  • seqadmin
                    Essential Discoveries and Tools in Epitranscriptomics
                    by seqadmin




                    The field of epigenetics has traditionally concentrated more on DNA and how changes like methylation and phosphorylation of histones impact gene expression and regulation. However, our increased understanding of RNA modifications and their importance in cellular processes has led to a rise in epitranscriptomics research. “Epitranscriptomics brings together the concepts of epigenetics and gene expression,” explained Adrien Leger, PhD, Principal Research Scientist...
                    Yesterday, 07:01 AM
                  • seqadmin
                    Current Approaches to Protein Sequencing
                    by seqadmin


                    Proteins are often described as the workhorses of the cell, and identifying their sequences is key to understanding their role in biological processes and disease. Currently, the most common technique used to determine protein sequences is mass spectrometry. While still a valuable tool, mass spectrometry faces several limitations and requires a highly experienced scientist familiar with the equipment to operate it. Additionally, other proteomic methods, like affinity assays, are constrained...
                    04-04-2024, 04:25 PM

                  ad_right_rmr

                  Collapse

                  News

                  Collapse

                  Topics Statistics Last Post
                  Started by seqadmin, 04-11-2024, 12:08 PM
                  0 responses
                  55 views
                  0 likes
                  Last Post seqadmin  
                  Started by seqadmin, 04-10-2024, 10:19 PM
                  0 responses
                  52 views
                  0 likes
                  Last Post seqadmin  
                  Started by seqadmin, 04-10-2024, 09:21 AM
                  0 responses
                  45 views
                  0 likes
                  Last Post seqadmin  
                  Started by seqadmin, 04-04-2024, 09:00 AM
                  0 responses
                  55 views
                  0 likes
                  Last Post seqadmin  
                  Working...
                  X