5 提取蛋白编号信息

cp

@omicsgene 老师,按照视频的脚本仔细核对没有错误,编译执行后也有没有报错,但pep.fa文件一直是空,文件名字没有空格,是不是文件读取权限向题


如果脚本有问题,编译起来就会报错,还得请教老师来解决,谢谢




use Bio::SeqIO;
use Bio::Seq;


my$in = Bio::SeqIO->new(-file => "$ARGV[1]" , -format => 'Fasta');


my$out = Bio::SeqIO->new(-file => "$ARGV[2]" , -format => 'Fasta');


my%keep=();


open IN, "$ARGV[0]" or die "$!";
while(<IN>)
{
chomp;
next if /^#/;
my@tmp=split(/\s+/);
$keep{"$tmp[1].1"}=1;
}
close(IN);
while (my $seq = $in->next_seq() )
{
my($id,$sequence,$desc)=($seq->id,$seq->seq,$seq->desc);
if(exists $keep{$id})
{
$out->write_seq($seq);
}


}
$in->close();
$out->close();
请先 登录 后评论

2 个回答

cp

不好意思,第一次提问,下次一定写的详细来问您

我是这样运行的:

先创建了pep.fa文件

然后执行

[share] perl get_fa_by_id.pl AT.gff  TAIR10.pep.all.fa pep.fa

attachments-2018-10-TxGlb2hJ5bc4274e341b9.jpg

AT.gff是之前提取的基因位置信息,没有问题

运行完没有报错,但

pep.fa文件是0字节

请先 登录 后评论
omicsgene - 生物信息
擅长:重测序,遗传进化,转录组,GWAS

get_fa_by_id.pl 这个脚本是根据第一个文件信息的第一列ID信息,提取第二个文件fasta文件的中对应的序列信息,存储到第三个文件当中;所有两个文件当中的ID对应才能得到想要的序列;

分析代码:

代码中这个地方 给ID都加了个".1"  你看看你的gff 文件里面的第一列ID都加 .1 后的基因ID,与你提供的所有蛋白质序列里面的序列ID能不能对应?如果不能对应请改正。


attachments-2018-10-qks5F1eQ5bc4561458370.jpg

请先 登录 后评论