前段时间做了一个克隆文库，有点大，拿到上百条序列，有了序列后，想选择一种合适的酶进行酶切，进而将克隆的结果与T-RFLP的结果联系起来。试了几个软件和一些在线酶切分析工具，好像都是可以用多种酶去切一个序列，而不能同时切多个序列。要是序列不多，一个个切切还可以，可是这上百条序列就需要几百次甚至更多的重复操作，不但能把手累的抽筋，眼都能被累的抽筋，实在让人无法忍受。于是决定自己做个小工具。

首先用了几个小时速成了一下Perl语言及Bioperl，发现这玩意处理字符串真的是太方便了！只需要下面这几行代码就可以自动或半自动同时查找多个序列酶切位点了(由于T-RFLP结果只与第一个酶切位点有关，因此只找每个序列的第一个切点)。

#!/usr/bin/perl -w
use strict;
use Bio::SeqIO;
my $file = 'All_Sequences.fas';
my $site = 'CCGG';
my $n = 2;
my $inputfile = Bio::SeqIO->new(-format => 'fasta', -file => $file);
my $outputfile = "Result.txt";
unless ( open(RESULT, ">$outputfile") ) {
print "Cannot open \"$outputfile\" to write to!!\n\n";
exit;
}
while ( my $seq = $inputfile->next_seq) {
my $str = $seq->seq;
if ( $str =~ /$site/ ) {
print RESULT index ($str, $site, 0)+$n, "\n";
print $seq->id, " -> ", index ($str, $site, 0)+$n,"\n";
}
}
close(RESULT);

使用说明：

首先要先安装Perl，另外还需要安装Bioperl。然后将上面这段代码保存到一个文本文档中，将后缀名改为.pl，把待切的多条序列合并到一起(fasta格式)，命名为All_Sequences.fas，和这个.pl文件放在同一个文件夹，然后和运行其它perl程序一样运行就可以了。需要说明的是上面只是用一个酶(MSP I)去切放在All_Sequences.fas这个文件中的所有序列，结果保存打印在屏幕上同时保存在Result.txt这个文件中，如果需要用别的酶，my $site = 'CCGG';和my $n = 2;这两行需要进行相应的修改，具体怎么改我就不说了，做酶切的同学肯定知道怎么改。

显而易见，这个东西使用起来不是很方便，因此如果需要处理的序列不是很多，使用这个代码没有什么优势，不需要使用，序列多的话，可以考虑用一下。

有人知道用什么软件可以进行类似的操作吗？如果知道，麻烦留言告诉我一下:)

Tags: , , ,