Unconfigured Ad

Collapse
X
 
  • Filter
  • Time
  • Show
Clear All
new posts
  • eni
    Junior Member
    • Dec 2009
    • 7

    Issues generating sam file format in ssaha2?

    Hi
    I'm trying to map transcriptome data on a refference genome and generate the output in SAM format. The problem: when I run it, says:
    Code:
    SSAHA2: Reading hashtable reference into memory.
    SSAHA2: Hashtable size      : 0.001007 GB
    SSAHA2: Number of sequences : 1
    SSAHA2: K-mer length        : 13
    SSAHA2: Skip step           : 3
    SSAHA2: Max. sequence length: 755119
    position_table - Query probably contains unusual characters!
    : Success
    and the sam file is empty.
    What does this mean? and How can I fix this?

    Thanks

    P.D.: I'm running the 2.5 version.
  • eni
    Junior Member
    • Dec 2009
    • 7

    #2
    More data:
    When I run using normal output, it stops while analizing reads like this:
    Code:
    @F35ERS101ARI6V
    GTAATAATGATGATATTAGTAGTAGTAGTAGTAGTAGTAGTGCTAGTAGTAGTAGTAGTAGTAGTAGTAGTAGTGCGAGTAGTAGTAGCACTAATAGTAGAAGTAAAGGAAAAAGTAAAA
    +
    HIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHCCB;9---,,63,,,,
    I don't know to much about fastq format, but comparing with reads that don't cause problems, the last line looks quite suspicious (related to de quality data, rigth).

    And the line corresponding to that read in the .qual file that I used to build the .fastq file is:
    Code:
    >F35ERS101ARI6V
    39 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 39 34 34 33 26 24 12 12 12 11 11 21 18 11 11 11 11

    The script I used to generate the fastq file is this:
    Code:
    #!/usr/bin/perl -w
     
    use strict;
     
    use Bio::SeqIO;
    use Bio::Seq::Quality;
     
    use Getopt::Long;
     
    die "pass a fasta and a fasta-quality file\n"
      unless @ARGV;
     
     
    my ($seq_infile,$qual_infile)
      = (scalar @ARGV == 1) ?($ARGV[0], "$ARGV[0].qual") : @ARGV;
     
    ## Create input objects for both a seq (fasta) and qual file
     
    my $in_seq_obj =
      Bio::SeqIO->new( -file   => $seq_infile,
    		   -format => 'fasta',
    		 );
     
    my $in_qual_obj =
      Bio::SeqIO->new( -file   => $qual_infile,
    		   -format => 'qual',
    		 );
     
    my $out_fastq_obj =
      Bio::SeqIO->new( -format => 'fastq'
    		 );
     
    while (1){
      ## create objects for both a seq and its associated qual
      my $seq_obj  = $in_seq_obj->next_seq || last;
      my $qual_obj = $in_qual_obj->next_seq;
     
      die "foo!\n"
        unless
          $seq_obj->id eq
    	$qual_obj->id;
     
      ## Here we use seq and qual object methods feed info for new BSQ
      ## object.
      my $bsq_obj =
        Bio::Seq::Quality->
    	new( -id   => $seq_obj->id,
    	     -seq  => $seq_obj->seq,
    	     -qual => $qual_obj->qual,
    	   );
     
      ## and print it out.
      $out_fastq_obj->write_fastq($bsq_obj);
    }
    So, what's the problem: my reads? the bioperl script? both?
    If it's the bioperl script, is there another format converter arround? because I don't want to discard reads if there is another way to fix this...HELP!!!

    Comment

    • jkbonfield
      Senior Member
      • Jul 2008
      • 146

      #3
      The quality line in fastq looks fine - with maximum character of "I" indicating Q40, as the .qual file states.

      I can't see anything obviously wrong with it. Maybe Ssaha2 thinks it doesn't look like a standard DNA sequence due to the almost complete lack of C, but that would be harsh.

      Comment

      • eni
        Junior Member
        • Dec 2009
        • 7

        #4
        Thanks for the answer.
        After I saw your post, I used the reads in fasta format and the program ran smoothly so it seems the sequences themselves are not the problem...

        Another question: Does SSAHA2 consider quality data while calculating the alignments scores or anything that could alter the results? I ask this because if not, I colud use the results obtained with the reads in fasta format an put the quality data in the .sam file...

        Comment

        Latest Articles

        Collapse

        • SEQadmin2
          Cancer Drug Resistance: The Lingering Barrier to Rising Survival
          by SEQadmin2



          Cancer survival rates have significantly increased in the last few decades in the United States, reaching a combined 70% 5-year survival rate by 2021. Behind this number, there are years of research to find new therapies, drug targets, and early detection methods. But there is one core challenge that keeps slowing down these advances, and it’s about drug resistance.

          There is no single reason why many patients don’t respond to treatment as expected. Cancer is...
          Today, 05:17 AM
        • GATTACAT
          Reply to Nine Things a Sample Prep Scientist Thinks About Before Sequencing
          by GATTACAT
          Love this - good data definitely starts from good input, and poor input can only give relatively poor data. I particularly like the mention of Nanodrop/absorbance based methods for quantification. It's such a toss up if you'll get an accurate reading or what amounts to a randomly generated number, and a lot of library/sequencing related issues can be traced back to poor quant.
          07-01-2026, 11:43 AM
        • SEQadmin2
          Nine Things a Sample Prep Scientist Thinks About Before Sequencing
          by SEQadmin2


          I’m not a sequencing expert. I’m a purification scientist who uses NGS to evaluate workflows my group develops. With this perspective, we think about the sample first and the NGS workflow second. The sequencer is an exceptionally honest reporter, but it can only report on what you give it, so whether you get clean, interpretable data from an NGS workflow is largely determined before you begin.

          Here are nine questions we think about, in roughly the order they matter, before...
          06-18-2026, 07:11 AM

        ad_right_rmr

        Collapse

        News

        Collapse

        Topics Statistics Last Post
        Started by SEQadmin2, Today, 10:08 AM
        0 responses
        5 views
        0 reactions
        Last Post SEQadmin2  
        Started by SEQadmin2, Yesterday, 11:05 AM
        0 responses
        7 views
        0 reactions
        Last Post SEQadmin2  
        Started by SEQadmin2, 07-02-2026, 11:08 AM
        0 responses
        30 views
        0 reactions
        Last Post SEQadmin2  
        Started by SEQadmin2, 06-30-2026, 05:37 AM
        0 responses
        28 views
        0 reactions
        Last Post SEQadmin2  
        Working...