去掉特定长度的短序列&取motif上下游一定长度的序列

/3 评论/在： /通过：

假如有大量的Fasta格式的多序列，假如要做以下处理：
问题1 去掉特定长度的短序列
若要去除其中小于五十碱基的序列，请问如何操作；
问题2 取motif上下游一定长度的序列（含motif）
假如我要取每个序列中motif为“AAAA”及其上下游10个碱基的序列片段，并输出位置信息，请问如何编程处理

第一个问题，除去 < 50 bp 的序列，可以直接用 bioperl

 

use Bio::SeqIO;

my $o_seqi = Bio::SeqIO->new(
   -file => $infile,
   -format => ‘fasta’,
);

my $o_seqo = Bio::SeqIO->new(
   -file => “>$outfile”,
   -format => ‘fasta’,
);

while (my $o_seq = $o_seqi->next_seq) {
   next if ($o_seq->length < 50);
  
   $o_seqo->write_seq($o_seq);
}

 

 

第二个，同样用 bioperl

 

use Bio::SeqIO;

my $o_seqi = Bio::SeqIO->new(
-file => $infile,
-format => ‘fasta’,
);

my $pattern = ‘AAAA’;

while (my $o_seq = $o_seqi->next_seq) {
   my $seq = $o_seq->seq;   # 提取序列，成为一个字符串

   if ($seq =~ /(.{10}$pattern.{10})/) {
     print $1;
   }
}

 

有点相关的文章

3 回复
  3. tempo8 says:

    motif为“AAAA”及其上下游10个碱基code有问题！
    借用http://liucheng.name/1285/，代码如下
    01 $/ = “>”;
    02 print “ID\tMotif\n”;
    03 while(){
    04 if($_ =~ /(.*?)\n(.*)/ms){
    05 $id = $1;
    06 $seq = $2;
    07 $seq =~ s/\n//g;
    08 while ($seq =~ /(.{0,10}NNNNN.{0,10})/g) {
    09 my $motif = $1;
    10 print “$id\t$motif\n”;
    11 }
    12 }
    13 }

评论已关闭。