Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • sort

    can any one help me , how to sort a file to have only unique reads, the format is the following

    >HWI-EAS373:1:1:10:1119#0/1
    GGAGCATTGCGACATTTAGCTNATTCAGCGGCATGGAAGG
    >HWI-EAS373:1:1:10:1119#0/2
    GGAATAAGCAGAGTGAGCATAAGAGAAGATGNCTTCATGC
    >HWI-EAS373:1:1:11:1157#0/1
    CGGGAATCCGGGACTTGAAAANATTCAGAACACTTTTCAG
    >HWI-EAS373:1:1:11:1157#0/2
    GATCCAGTTAAGGCAGAGATTGAGAGTTCACNGGAGAGTA
    >HWI-EAS373:1:1:11:986#0/1
    GTTTCCTCATACTAGACAACANATTCATGGTGTTTGAAGC
    >HWI-EAS373:1:1:11:986#0/2
    CCGGATGAGAGAAAAAGGGCTGCTGCAGCCANGGTTGAGC
    >HWI-EAS373:1:1:11:1274#0/1
    CGTCCTCCTCTCCTCCTGTCGNCTCCAGCGCACCA

  • #3
    Note that there weren't any repeated sequences in your example dataset, but something like this should pull out unique sequences. It will be a bit memory heavy since it needs to keep all unique sequences in memory but there's no getting around that really.

    Code:
    #!/usr/bin/perl
    use warnings;
    use strict;
    
    my %kept;
    my %rejected;
    
    while(my $header = <>) {
      chomp $header;
      my $seq = <>;
      chomp $seq;
      next if (exists $rejected{$seq});
      if (exists $kept{$seq}) {
        delete $kept{$seq};
        $rejected{$seq} = 1;
      }
      else {
        $kept{$seq} = $header;
      }
    }
    
    foreach (keys %kept) {
      print $kept{$_},"\n",$_,"\n";
    }
    
    warn "Rejected ".(scalar keys %rejected)." sequences\n";

    Comment


    • #4
      grep -v '>' myfasta | sort -u

      is a short solution. But perhaps you need fasta? with or without original ids?

      Comment


      • #5
        Originally posted by simonandrews View Post
        Note that there weren't any repeated sequences in your example dataset, but something like this should pull out unique sequences. It will be a bit memory heavy since it needs to keep all unique sequences in memory but there's no getting around that really.
        Using a combination of Bio::SeqIO and Digest::MD5 you can generate signatures for all the sequences traversed and only keep the MD5 signatures in memory. Works best when the average sequence is much longer than the MD5 signature.

        The code is left as a homework exercise for bioenvisage.

        Comment


        • #6
          We've implemented an efficient two-pass algorithm in Goby (see my other post). Here's how to proceed once you have downloaded the tool:

          If you have a fasta/fastq format, first convert to compact format:
          java -Xmx3g -jar goby.jar -m fasta-to-compact with-redundancy.fasta
          (The file with-redundancy.compact-reads should have been created.)
          Use the tally-reads mode to calculate how many times each sequence appears in the input:
          java -Xmx3g -jar goby.jar -m tally-reads with-redundancy.compact-reads -f myfilter
          Tally reads leverages sequence digests and works in two passes to minimize memory usage. Input files can contain tens of millions of reads.
          Convert back to fasta, excluding sequences that appear more than once:
          java -Xmx3g -jar goby.jar -m compact-to-fasta -i with-redundancy.compact-reads -f myfilter-keep.filter -o unique-reads.fa

          See our wiki for nicer formatting and other options:
          Last edited by Fabien Campagne; 02-02-2010, 02:24 PM.

          Comment


          • #7
            awk will do a much faster job than sort | uniq -c (which I used a lot before I discovered awk). First, use grep to get the lines not starting with '>', then use awk. Same as with @simonandrews solution: it will store all unique reads in memory.

            grep -v '>' INFILE | awk '{r[$0]=1}END{for (x in r){print x}}'

            Comment

            Latest Articles

            Collapse

            • seqadmin
              Recent Advances in Sequencing Analysis Tools
              by seqadmin


              The sequencing world is rapidly changing due to declining costs, enhanced accuracies, and the advent of newer, cutting-edge instruments. Equally important to these developments are improvements in sequencing analysis, a process that converts vast amounts of raw data into a comprehensible and meaningful form. This complex task requires expertise and the right analysis tools. In this article, we highlight the progress and innovation in sequencing analysis by reviewing several of the...
              05-06-2024, 07:48 AM
            • 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...
              04-22-2024, 07:01 AM

            ad_right_rmr

            Collapse

            News

            Collapse

            Topics Statistics Last Post
            Started by seqadmin, 05-10-2024, 06:35 AM
            0 responses
            20 views
            0 likes
            Last Post seqadmin  
            Started by seqadmin, 05-09-2024, 02:46 PM
            0 responses
            25 views
            0 likes
            Last Post seqadmin  
            Started by seqadmin, 05-07-2024, 06:57 AM
            0 responses
            21 views
            0 likes
            Last Post seqadmin  
            Started by seqadmin, 05-06-2024, 07:17 AM
            0 responses
            21 views
            0 likes
            Last Post seqadmin  
            Working...
            X