Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • Samtools fasta line length segfault, redux

    There's been at least one thread on this already, so I'm wondering if anyone has a handy shell script they could share for fixing .fa files with bad line lengths. In my particular case, I get this error from samtools faidx:

    [fai_build_core] different line length in sequence 'chr13'.

    ...when using our standard reference file, hsap_36.1_hg18.fa.

    I want to reformat this file rather than downloading some known-correct one, since I don't know exactly where our file came from and what differences there might be between it and another.

    The thread linked above has a snippet of BioPerl, but alas, I don't have that available on the server where I'm working. The last poster outlines an approach, but doesn't provide actual code.

    Hopefully I'm not being too lazy here by asking if someone's got such a script in their back pocket...I'm just a bit overwhelmed at the moment. I would love not to have to solve this particular problem from scratch, and I'm sure a solution posted here would be useful to others in the future. Thanks!

  • #2
    Turns out that in my case it's actually a blank line that's causing the problem. I started wondering why the issue was coming up at all, given that awk says the longest line in the file is 70 characters long:

    awk '{ if (x < length()) x = length() } END { print x }' hsap_36.1_hg18.fa

    I thought surely it couldn't be complaining about short lines...so I looked for sub-70-character lines:

    awk 'length($0) < 70' hsap_36.1_hg18.fa

    I got the expected set of header lines, and one trailing incomplete line per header block. But the last line of the chr13 block was blank due to the actual final line of sequence being exactly 70 characters.

    So the brute-force solution for me is just to eliminate blank lines from the file:

    grep . hsap_36.1_hg18.fa > tmp; mv tmp hsap_36.1_hg18.fa

    This still leaves the questions of (1) why samtools can't handle a simple silly blank line and (2) how to fold a FASTA file that does have anomalous line lengths...but my immediate crisis is over. Thanks in advance, though, for any insights anyone has on either of those questions!

    Comment


    • #3
      Thanks.

      I did this:
      Code:
      awk '(NF!=0){print}' sample.fa>tmp.fa
      awk '(NF==0){print "\blank"}' tmp.fa # check
      mv tmp.fa sample.fa

      Comment


      • #4
        That works too. But there are canned answers to this sort of question nowadays, fortunately, that didn't exist in 2009 when I originally posted. Picard, for example, has a FASTA-fixer:



        I suspect BioPerl and similar packages were doing this even in 2009, but I was pretty new to the field then and didn't know about them.

        Comment


        • #5
          In my case, the actual problem was different (I found out later).

          I extracted the genomic sequences from UCSC table browser in fasta format and saved it as my reference file. Invoking samtools faidx was giving the error 'different line length...'. However, fixing empty lines did not stop the error, and I discovered that my fasta file was truncated due to timeout while downloading from UCSC; the error message was the last lines of my fasta file and it was causing samtools emit this message.

          Comment


          • #6
            I found another case where this was happening I thought I'd just add to the list for future reference.

            My input fasta file was line wrapped, and there was a blank line between sequences. The problem sequence was where coincidentally the seq ended exactly at the line wrap position, ie:

            Code:
            >not_a_problem_seq
            AATCGACGTACGTAGCTGATC
            ATCGATGCTAGCTATAATCGT
            ACTAGCTACTGG
            
            >problem_seq
            ACTGCTAGCTGATCGATCGTAG
            ACGTACGTAGCTAGCTGACTGA
            ACTGATCGTAGCTAGCTGATCG <- here!
            
            >next_seq
            ATCGACGTAGCTAGCTAGCTAT
            ACTGATGCTACTAGCTGATGCT
            ACTGA
            Issue was resolved by actually adding an extra blank line between the end of problem_seq and the header for next_seq.

            Comment


            • #7
              Originally posted by reubennowell View Post
              My input fasta file was line wrapped, and there was a blank line between sequences. The problem sequence was where coincidentally the seq ended exactly at the line wrap position
              Thanks for reporting this. I've now fixed it and the fix will appear in samtools (and htslib) 1.3.

              Comment

              Latest Articles

              Collapse

              • seqadmin
                Strategies for Sequencing Challenging Samples
                by seqadmin


                Despite advancements in sequencing platforms and related sample preparation technologies, certain sample types continue to present significant challenges that can compromise sequencing results. Pedro Echave, Senior Manager of the Global Business Segment at Revvity, explained that the success of a sequencing experiment ultimately depends on the amount and integrity of the nucleic acid template (RNA or DNA) obtained from a sample. “The better the quality of the nucleic acid isolated...
                03-22-2024, 06:39 AM
              • seqadmin
                Techniques and Challenges in Conservation Genomics
                by seqadmin



                The field of conservation genomics centers on applying genomics technologies in support of conservation efforts and the preservation of biodiversity. This article features interviews with two researchers who showcase their innovative work and highlight the current state and future of conservation genomics.

                Avian Conservation
                Matthew DeSaix, a recent doctoral graduate from Kristen Ruegg’s lab at The University of Colorado, shared that most of his research...
                03-08-2024, 10:41 AM

              ad_right_rmr

              Collapse

              News

              Collapse

              Topics Statistics Last Post
              Started by seqadmin, Yesterday, 06:37 PM
              0 responses
              10 views
              0 likes
              Last Post seqadmin  
              Started by seqadmin, Yesterday, 06:07 PM
              0 responses
              9 views
              0 likes
              Last Post seqadmin  
              Started by seqadmin, 03-22-2024, 10:03 AM
              0 responses
              49 views
              0 likes
              Last Post seqadmin  
              Started by seqadmin, 03-21-2024, 07:32 AM
              0 responses
              67 views
              0 likes
              Last Post seqadmin  
              Working...
              X