Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • Extract unaligned reads (Tophat) from FastQ

    I would like to further examine reads that Tophat didn't manage to align in a first run and i wonder, if there is any easy way to get these reads. With Bowtie this would be easy using the "--un" argument, but Tophat doesn't seem to have smth like this. I am so far able to extract read-ids of the reads that do align by:

    Code:
    cut --fields=1 accepted_hits.sam | sort --unique > accepted_hits_readsIds.txt
    From that point i'd need to extract fastq-entries that don't match any of the lines in the readIds file. Since FastQ-entries do not consist of single lines i got stuck here - any ideas/help would be appreciated!

    Thanks in advance & Cheers
    Uwe
    Last edited by Uwe Appelt; 09-15-2010, 03:10 AM.

  • #2
    you can try
    available at no cost under a non-open-source license by requesting from the web-site; Binary: available for direct download at no cost. For-Profit: Submit request for for-profit license from the web-site.

    Bioinformatics. 2010 Jul 6. [Epub ahead of print]
    G-SQZ: Compact Encoding of Genomic Sequence and Quality Data.

    Tembe W, Lowey J, Suh E.

    Translational Genomics Research Institute, 445 N 5th Street, Phoenix, AZ 85004, USA.
    Abstract

    SUMMARY: Large volumes of data generated by high-throughput sequencing instruments present non-trivial challenges in data storage, content access, and transfer. We present G-SQZ, a Huffman coding-based sequencing-reads specific representation scheme that compresses data without altering the relative order. G-SQZ has achieved from 65% to 81% compression on benchmark datasets, and it allows selective access without scanning and decoding from start. This paper focuses on describing the underlying encoding scheme and its software implementation, and a more theoretical problem of optimal compression is out of scope. The immediate practical benefits include reduced infrastructure and informatics costs in managing and analyzing large sequencing data. AVAILABILITY: http://public.tgen.org/sqz. Academic/non-profit: Source: available at no cost under a non-open-source license by requesting from the web-site; Binary: available for direct download at no cost. For-Profit: Submit request for for-profit license from the web-site. CONTACT: Waibhav Tembe ([email protected]).

    or maybe use bioperl

    if the number of reads are not a lot
    http://kevin-gattaca.blogspot.com/

    Comment


    • #3
      You probably need to do this in two passes. It's also a pain that tophat seems to alter the sequence id (but maybe this is because I was using paired end data?), so you have to ajdust the ids a bit.

      The code below seems to work on the tophat files I just ran it against.

      Code:
      #!/usr/bin/perl
      use warnings;
      use strict;
      
      my ($fastq,$sam,$outfile) = @ARGV;
      
      unless ($outfile) {
        die "Usage is filter_unmapped_reads.pl [FastQ file] [SAM File] [File for unmapped reads]\n";
      }
      
      if (-e $outfile) {
        die "Won't overwrite an existing file, delete it first!";
      }
      
      open (FASTQ,$fastq) or die "Can't open fastq file: $!";
      open (SAM,$sam) or die "Can't open SAM file: $!";
      open (OUT,'>',$outfile) or die "Can't write to $outfile: $!";
      
      my $ids = read_ids();
      
      filter_fastq($ids);
      
      close OUT or die "Can't write to $outfile: $!";
      
      
      sub filter_fastq {
      
        warn "Filtering FastQ file\n";
      
        my ($ids) = @_;
      
        while (<FASTQ>) {
      
          if (/^@(\S+)/) {
            my $id = $1;
      
            # Remove the end designator from paired end reads
            $id =~ s/\/\d+$//;
      
            my $seq = <FASTQ>;
            my $id2 = <FASTQ>;
            my $qual = <FASTQ>;
      
      
            unless (exists $ids->{$id}) {
      	print OUT $_,$seq,$id2,$qual;
            }
          }
          else {
            warn "Line '$_' should have been an id line, but wasn't\n";
          }
      
        }
      
      }
      
      
      sub read_ids {
      
        warn "Collecting mapped ids\n";
      
        my $ids;
      
        while (<SAM>) {
      
          next if (/^@/);
          my ($id) = split(/\t/);
          $ids->{$id} = 1;
        }
      
        return $ids;
      }

      Comment


      • #4
        Hi Simon,

        thank you so much for that chunk of code, it works like a charme! I worked out a solution of my own as well, but besides the obvious (e.g. poor parsing) drawbacks, fgrep appears to consume 40Gb of RAM in order to filter for ~18e6 read ids (in accepted_readIds.txt).

        Code:
        cut --fields=1 ./tophat_out/accepted_hits.sam | sort --unique > ./tophat_out/accepted_readIds.txt
        paste - - - - < ./reads_1.fq | fgrep --invert-match --file ./tophat_out/accepted_readsIds.txt | tr "\t" "\n" > ./readsFilt_1.fq
        So thanks again and Cheers,
        Uwe

        Comment


        • #5
          Simon, that script works wonderfully. Thanks.

          Comment


          • #6
            Excellent Simon. Thanks!

            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
            12 views
            0 likes
            Last Post seqadmin  
            Started by seqadmin, Yesterday, 06:07 PM
            0 responses
            10 views
            0 likes
            Last Post seqadmin  
            Started by seqadmin, 03-22-2024, 10:03 AM
            0 responses
            51 views
            0 likes
            Last Post seqadmin  
            Started by seqadmin, 03-21-2024, 07:32 AM
            0 responses
            68 views
            0 likes
            Last Post seqadmin  
            Working...
            X