从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
  • 阅读 ( 1039 )
  • 分类:perl

0 条评论

请先 登录 后评论
omicsgene
omicsgene

生物信息

252 篇文章

作家榜 »

  1. omicsgene 252 文章
  2. 安生水 212 文章
  3. Daitoue 165 文章
  4. 生物女学霸 94 文章
  5. landy 37 文章
  6. 组学生物-王运斌 34 文章
  7. 红橙子 33 文章
  8. omics007 22 文章