vsearch 去除重復序列和singleton 序列

在16S數據分析中，為了減少聚類的時間，提高准確度，需要去除重復序列，而singleton序列因為沒有其他的序列作為驗證，可信度不是很高，也需要去除，通常情況下使用usearch 完成這2項任務，但是usearch 64位是收費的，而32為的usearch 在64位的red hat 上測試時，去除重復序列時報錯了，libgomp: Thread creation failed: Resource temporarily unavailable

百度之后了解到是由於進程數達到上限，修改了上限后還是報錯，就放棄使用usearch 來去除重復序列，使用vsearch 來去除重復序列

源代碼：

　　https://github.com/torognes/vsearch/releases

安裝：

　　wget https://github.com/torognes/vsearch/archive/v1.11.1.tar.gz

      tar xzvf v1.11.1

      cd vsearch-1.11.1/

      ./autogen.sh

      ./configure

      make      

      make install

 

測試：

　　# 去除重復序列

      vsearch --derep_fulllength raw.reads.fasta  --output test.fasta --sizeout        

            

      # 去除singleton序列

      vsearch --sortbysize test.fasta  --output desingleton.fasta --minsize 2

 

vsearch 和 mothur 去除重復序列的結果完全一致，其實去除重復序列就是就將那些序列完全一致的序列只取一條就可以了，去除singleton 序列就是將那些只出現一次的序列去除，為了加深理解，我寫了個perl腳本

來完成這2個任務，經過測試和vsearch的結果一致，代碼如下：

#!/usr/bin/perl

use warnings;
use strict;

my ($fasta) = @ARGV;

my %unique = ();
local $/ = ">";
open FASTA, $fasta or die "Can't open $fasta\n";
while (<FASTA>) {
    chomp;
    next if /^\s*$/;
    my ($id, $seq) = split /\n/, $_, 2;
    $seq =~ s/\s+//;
    push @{$unique{$seq}}, $id;
}
close FASTA;


foreach my $x (keys %unique) {
    my $size = scalar @{$unique{$x}};
    unshift  @{$unique{$x}}, $size;
}

foreach my $x (sort {$unique{$b}->[0] <=>  $unique{$a}->[0] } keys %unique) {
    my $id   = $unique{$x}->[1];
    my $size = $unique{$x}->[0];
    next if $size = 1;
    print qq{>$id;size=$size;\n$x\n};

}