从gff文件当中提取对应转录本ID的基因结构信息用于GSDS绘制结构图脚本更新

get_gene_exon_from_gff.pl 脚本更新,

脚本名称:get_gene_exon_from_gff.pl

用法:perl ../script/get_gene_exon_from_gff.pl -in1 WRKY_domain_new_out_emoved_redundant.txt -in2 Arabidopsis_thaliana.TAIR10.31.gff3 -out gene_exon_info.gff

输出结果(部分):

attachments-2019-03-MpLqEWYv5c8b386f199eb.jpg


脚本源代码:

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/\t/, $_;
	$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;
}


脚本下载地址:get_gene_exon_from_gff.pl  可下载替换以前的版本。
  • 发表于 2019-03-15 13:34
  • 阅读 ( 601 )
  • 分类:perl

0 条评论

请先 登录 后评论
omicsgene
omicsgene

生物信息

223 篇文章

作家榜 »

  1. omicsgene 223 文章
  2. 安生水 187 文章
  3. Daitoue 154 文章
  4. 生物女学霸 42 文章
  5. landy 37 文章
  6. 生信老顽童 33 文章
  7. orange 32 文章
  8. omics007 18 文章