perl提取genebank序列，程序写好了，但有点问题，望各位大哥帮忙看看

chenhao392

1. 请不要截图用序列A, 底下贴序列用序列B..
    test data 和error要用一份data...还要我自己去NCBI找你这个柑橘病毒.

2. 贴代码请加CODE的tag.

3.  猜测error来自 &seq(), 明显靠前的地方可以substr，后面找不到了. 就是说你的正则没拿到全序列.

4. 你就不该用此格式拿序列。就下fasta，另外想办法拿CDS的注释，找个GFF之类的，或者再parse这个文件找CDS位置就好.

b114213903

是不是再看看Pubmed网站的使用方式？应该可以直接下载对就的CDS区！

yzkcxt

谢谢，程序没有问题，是NCBI别人登录的序列存在问题，除了atcg，还有简并碱基，所有读取不出来。换成\w就好了。回复 2# chenhao392


   

yzkcxt

#!/usr/bin/perl
#2014/8/22-10:31:54
open(FILE1,"<CTV.gb") or die;
open(FILE3,">CTV.txt") or die;
while(<FILE1>){
        push @all,$_;
}
        $all=join "",@all;
        $all =~s/ +//mg;
$all =~ s/\s+//mg;
$all =~s/\/\//}/mg;
@private = split "}",$all;
foreach $a(@private){
$a =~/ACCESSION(?<first>(\w{1,4}_?\d{3,6}))/;
$name1 = $+{first};
$i = 0;
while($a =~ /CDS(\d+\.\.\d+)/mg){
$pos = $1;
$i += 1;
open(FILE2,">>CTV_$i.txt") or die;
$pos =~ s/\.\./;/;
($start,$end) = split ";",$pos;
print "$start,$end\n";
$sequence = &seq($a);
print FILE3 "\n>$name1\n$sequence";
my $sequ = substr($sequence,$start-1,$end-$start+1);
print FILE2 "\n>$name1($pos)\n$sequ";
}
}
sub seq{
$seq = shift(@_);
$seq =~ s/\d+//mg;
$seq =~ s/\s+//mg;
$seq =~ /ORIGIN([\w]+)/xsmg;
$sequen = ($1);
return $sequen;
}
close FILE1;
close FILE2;


回复 3# b114213903