基因组组装RNA-Seq 分析流程生物信息学

根据gff文件提取所有基因的位置信息

2021-06-17  本文已影响0人  MLD_TRNA

gff文件当中存储了基因组当中所有基因的注释信息,如果想得到基因组当中所有基因的位置信息可以利用awk命令批量的提取,命令如下:

$ grep -v '#' Arabidopsis_thaliana.TAIR10.41.gff3|awk -F "[\t=:;]" 'BEGIN{OFS="\t"}$3=="gene"{print $1,$4,$5,$10}' |head
#结果如下
1       3631    5899    AT1G01010
1       6788    9130    AT1G01020
1       11649   13714   AT1G01030
1       23121   31227   AT1G01040
1       31170   33171   AT1G01050
1       33365   37871   AT1G01060

GFF中提取所有转录本的位置信息,以及对应的基因ID信息;

提取基因的转录本信息,以及对应的基因ID;

提取基因组注释文件GFF中所有基因转录本的位置信息,以及转录本对应的基因的ID:

attachments-2018-12-1T3DdNEG5c126d45746fa.jpg

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);
上一篇下一篇

猜你喜欢

热点阅读