本篇内容主要讲解“perl如何提取GFF中所有转录本的位置信息”,感兴趣的朋友不妨来看看。本文介绍的方法操作简单快捷,实用性强。下面就让小编来带大家学习“perl如何提取GFF中所有转录本的位置信息”吧!
提取基因组注释文件GFF中所有基因转录本的位置信息,以及转录本对应的基因的ID:
perl代码如下:
#!/usr/bin/perl -w use strict; use Cwd qw(abs_path getcwd); use Getopt::Long; use Data::Dumper; die "perl $0 <gff> <outfile>" unless(@ARGV==2); my$gff=$ARGV[0]; my%gene=(); my%gene_region=(); my%mRNA2Gene=(); open IN,"$gff" or die "$!"; open OUT ,">$ARGV[1]" or die "$!"; print OUT "#mRNA_ID\tgene_ID\tchr\tstart\tend\tstrand\n"; while(<IN>){ chomp; next if (/^#/); my@tmp=split(/\t/); if($tmp[2] =~/^gene/){ my($id)=($tmp[8]=~/ID=([^;]+)/); $gene{$id}=1; $gene_region{$id}=[$tmp[0],$tmp[3],$tmp[4],$tmp[6]]; #print "gene:$id\n"; #my$gene_chr->{$id}=$tmp[0]; } if($tmp[2] =~/mRNA|transcript/i){ my($id)=($tmp[8]=~/ID=([^;]+)/); my($pid)=($tmp[8]=~/Parent=([^;]+)/); if(exists $gene{$pid}){ print OUT "$id\t$pid\t$tmp[0]\t$tmp[3]\t$tmp[4]\t$tmp[6]\n"; } #print "mRNA:$id\n"; } } close(IN); close(OUT);
到此,相信大家对“perl如何提取GFF中所有转录本的位置信息”有了更深的了解,不妨来实际操作一番吧!这里是亿速云网站,更多相关内容可以进入相关频道进行查询,关注我们,继续学习!
免责声明:本站发布的内容(图片、视频和文字)以原创、转载和分享为主,文章观点不代表本网站立场,如果涉及侵权请联系站长邮箱:is@yisu.com进行举报,并提供相关证据,一经查实,将立刻删除涉嫌侵权内容。