perl从多行文件中提取一段序列

2013-07-08 15:46 阅读(?)评论(0)

有序列文件:帮网友从一个文件中提取一度序列的脚本记录下来。

文件开头是:LOCUS       HQ912884                9667 bp    RNA     linear   VRL 18-OCT-2012一类,随后是蛋白序列,再后是核酸序列,需要把其中的一段核酸序列找出来。

功能:把原始序列提取出来到一个文件,把特定的这段长度也取出来。

perl代码:

#!/usr/bin/perl -w
use strict;

unless (@ARGV==3) {
    die"Usage: perl $0 <input.fa> <out.str1><out.str2>\n";
}

my ($infile,$outfile1,$outfile2) = @ARGV;
open IN,$infile || die"error: can't open infile: $infile";
open OUT1,">$outfile1" || die$!;
open OUT2,">$outfile2"||die$!;
$/="LOCUS";<IN>;

while ( my $seq = <IN>){
 my $id = $1 if($seq =~ /(\w+)/);  #id号赋值
 $seq =~ s/\n//g; #去掉读入每行的换行符
 $seq .= $seq; #用点号把读入的数据富集的一个标量$seq中
 $seq =~ s/[\d\D]*ORIGIN//g; #删除每一分隔符内容ORIGIN前的内容,获得sequence
 $seq =~ s/[\d]//g;#删除其中的数字
 $seq =~ s/\s//g;#删除其中的空白字符
 $seq =~ s/\/\/LOCUS//g;删除尾部的//LOCUS标记
 my$str = substr($seq,2405,1095); #取所需要的字符串
 print OUT1 "\n$id\n$seq\n";
 print OUT2 "\n$id\n$str\n";
}

$/="\n";
close IN;
close OUT1;
close OUT2;

perl $0 …………,命令运行后输出文件。搞定

作者已禁止网友对该文进行评论