利用gff文件提取基因结构,结果文件是空的

这是我的GFF文件

attachments-2020-08-ImwhNhwR5f3cc82467dc4.png

attachments-2020-08-ImwhNhwR5f3cc82467dc4.png下图是运行程序,运行过程中没报错,但结果文件是空的

attachments-2020-08-36UPzaPA5f3cc87c0f2db.png这是我的ID

attachments-2020-08-YZNWRKtW5f3cc90e74661.png

请先 登录 后评论

2 个回答

omicsgene - 生物信息
擅长:重测序,遗传进化,转录组,GWAS

脚本用错了吧,你可以把perl脚本贴一下我看看问题

请先 登录 后评论
Sirius

use Getopt::Long;

my %opts;

use Data::Dumper;

GetOptions( \%opts, "in1=s", "in2=s", "out=s", "h" );

if (   !defined( $opts{in1} )

|| !defined( $opts{in2} )

|| !defined( $opts{out} )

|| defined( $opts{h} ) )

{

&USAGE;

}

open( IN1, "$opts{in1}" )  || die "open $opts{in1} failed\n";

open( IN2, "$opts{in2}" )  || die "open $opts{in2} failed\n";

open( OUT, ">$opts{out}" ) || die "open $opts{out} failed\n";

my %gffs;

while (<IN1>) {

chomp;

next if /^#/;

my @b = split/\st/, $_;

$gffs{$b[0]} = 1;

}


#print Dumper(\%gffs);

while (<IN2>) {

chomp;

next if (/^#/);

my @a = split /\t/, $_;

next if $a[2]=~/exon/i;

if ($a[2] =~/^mRNA$/i or $a[2] =~/^transcript$/i ) {

($id1) =  ($a[8] =~ m/ID=([^;]*)/);


}elsif ( $a[2] =~/^CDS$/i or $a[2] =~/utr/i ) {


($id1) =  ($a[8] =~ m/Parent=([^;]*)/);

}else{

next;

}


if ( exists $gffs{$id1} ) {

print OUT "$_\n";

}


}

close OUT;

close IN1;

close IN2;


sub USAGE {

print "usage: perl $0 -in1  mRNA_id.txt -in2  genome.gff3  -out gene_location.txt ";

exit;

}

老师,这是全部脚本内容
请先 登录 后评论
  • 3 关注
  • 0 收藏,2476 浏览
  • Sirius 提出于 2020-08-19 14:39

相似问题