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


假如有大量的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;
   }
}

 


《 “去掉特定长度的短序列&取motif上下游一定长度的序列” 》 有 3 条评论

  1. 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 }