sort

[bioenvisage]

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

[alex_kozik]

I am using combination of Unix 'sort' and 'uniq' commands to generate non-redundant fasta file, details are here:

http://code.google.com/p/atgc-tools/wiki/regexp_and_unix_tricks

http://cadmium-gcat.blogspot.com/2009/12/non-redundant-fasta-dataset-using-sort.html

[simonandrews]

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.

#!/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";

[cariaso]

grep -v '>' myfasta | sort -u

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

[Dave S.]

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.

[Fabien Campagne]

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:

http://icb.med.cornell.edu/wiki/index.php/Goby/Non-redundant#Filter_redundant_reads

[flxlex]

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}}'