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

}