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

}

 

时间: 2024-10-10 09:09:55

vsearch 去除重复序列和singleton 序列的相关文章

基因组注释

基因组注释主要包括四个研究方向:重复序列的识别:非编码RNA的预测:基因结构预测和基因功能注释.我们将分别对这四个领域进行阐述. 1 重复序列的识别. 1.1  重复序列的研究背景和意义:重复序列可分为串联重复序列(Tendam repeat)和散在重复序列(Interpersed repeat)两大类.其中串联重复序列包括有微卫星序列,小卫星序列等等:散在重复序列又称转座子元件,包括以DNA-DNA方式转座的DNA转座子和反转录转座子(retrotransposon).常见的反转录转座子类别有

DNA甲基化及其测量方法(转)

转自声明的奥秘 www.lifeomics.com DNA甲基化与肿瘤发生:         DNA甲基化水平和模式的改变是肿瘤发生的一个重要因素.这些变化包括CpG岛局部的高甲基化和基因组DNA低甲基化状态.如图1左所示,在正常细胞中,位于抑癌基因启动子区域的CpG岛处于低水平或未甲基化状态,此时抑癌基因处于正常的开放状态,抑癌基因不断表达抑制肿瘤的发生.而在肿瘤细胞中,该区域的CpG岛被高度甲基化,染色质构象发生改变,抑癌基因的表达被关闭,从而导致细胞进入细胞周期,凋亡丧失,DNA修复缺陷,

扩增子分析QIIME2分析实战Moving Pictures

本示例的的数据来自文章<Moving pictures of the human microbiome>,Genome Biology 2011,取样来自两个人身体四个部位五个时间点 进入环境 source activate qiime2-2017.6 退出环境 source deactivate 准备数据 # 创建并进入工作目录 mkdir -p qiime2-moving-pictures-tutorialcd qiime2-moving-pictures-tutorial # 下载实验设

Git中的merge命令实现和工作方式

想象一下有如下情形:代码库中存在两个分支,并且每个分支都进行了修改,最后你想要将其中的一个分支合并到其他的分支中.个人博客网址 http://swinghu.github.com/ 那么要问合并的处理过程是怎么样的呢?Git是对每个分支,依据分支的历史数据按照序列化操作,还是它只是合并每个分支里文件的最后版本?这是一个问题,我想对git的merge操作有必要进行分析一下. 回忆一下,我们知道Git的版本库内部结构是以有向无环图(directed acyclic graph)组织起来的:每一次co

浅谈凸包之Andrew 与 Graham

前言 脑补知识点: 1.向量的内积(数量积,点乘): 公式:a· b = |a| * |b| cos<a, b>=a.x* b.y + b.x * a.y 2.向量的外积(向量积,差乘): 公式:|c|= |a|*|b|*sin<a, b> = a.x * b.y - b.x * a.y 点在多边形内判定 多边形: 就是二维平面上被一系列首尾相接.闭合的折线段围成的区域 在程序中一般用定点数组表示 其中各个定点按照逆时针顺序排序 问: 给你一个点 如何判断它是在多边形内 呢? 1.

《BI那点儿事》数据流转换——排序

原文:<BI那点儿事>数据流转换--排序 排序转换允许对数据流中的数据按照某一列进行排序.这是五个常用的转换之一.连接数据源打开编辑界面,编辑这种任务.不想设置为排序列的字段不要选中,默认情况下所有列都会选中.如图所示,按照TotalSugar_Cnt排序,并将所有列输出. 在底部的表格中,可以设置输出列的别名,是否按照列来排序.Sort Order列显示列将会第一排序,第二排序还是第三排序.双击列去除重复的排序列.

Linq专题之查询操作

前面我们主要讲解的是Linq的查询表达式,Linq不但提供了一些基本的查询表达式,还提供了数十个查询操作.比如筛选操作.聚合操作.投影操作等等.通过这些查询操作可以更方便的对数据源进行处理. Linq提供了数十个查询操作,大多数的操作都是针对实现了IQueryable<T>和IEnumerbale<T>接口的序列. 序号     查询操作           对应的查询表达式                 说明                                   

转录组分析的正确姿势

转录组分析的正确姿势 转录组分析是目前应用最广的高通量测序分析技术之一.常见设计是不同样品之间比较,寻找差异基因.标志基因.协同变化基因.差异剪接和新转录本,并进行结果可视化.功能注释和网络分析等. 转录组的测序分析也相对成熟,从RNA提取.构建文库.上机测序再到结果解析既可以自己完成,又可以在专业公司进行. 概括来看转录组的分析流程比较简单,序列比对-转录本拼接 (可选)-表达定量-差异基因-功能富集-定制分析.整个环节清晰流畅,可以作为最开始接触高通量测序学习最合适的技术之一. 但重点和难点

解读生命密码的基本手段 ——DNA测序技术的前世今生

解读生命密码的基本手段 ——DNA测序技术的前世今生 任鲁风  于军 (中国科学院基因组科学及信息重点实验室,北京基因组研究所) DNA(脱氧核糖核酸)和RNA(核糖核酸)是生命体的两种最基本组成物质,其序列的组成和变化造就了形形色色的生命世界.这两种承担了生命体遗传信息载体功能的物质,一方面在生命的不断繁衍中保持了各个物种的独特面目,另一方面又通过不断的演变改变着自身性状,同时又影响着与之相关的物种,这一规律在生命科学领域被归纳为“中心法则”.笼统而言,几乎全部的生命现象均来源于A.T.C.G