提取基因组注释文件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语言学习视频课程:《perl入门》《perl高级编程》
1. 文章越来越难发?是你没发现新思路,基因家族分析发2-4分文章简单快速,学习链接:基因家族分析实操课程、基因家族文献思路解读
2. 转录组数据理解不深入?图表看不懂?点击链接学习深入解读数据结果文件,学习链接:转录组(有参)结果解读;转录组(无参)结果解读
3. 转录组数据深入挖掘技能-WGCNA,提升你的文章档次,学习链接:WGCNA-加权基因共表达网络分析
4. 转录组数据怎么挖掘?学习链接:转录组标准分析后的数据挖掘、转录组文献解读
5. 微生物16S/ITS/18S分析原理及结果解读、OTU网络图绘制、cytoscape与网络图绘制课程
6. 生物信息入门到精通必修基础课:linux系统使用、perl入门到精通、perl语言高级、R语言画图
7. 医学相关数据挖掘课程,不用做实验也能发文章:TCGA-差异基因分析、GEO芯片数据挖掘、 GEO芯片数据不同平台标准化 、GSEA富集分析课程、TCGA临床数据生存分析、TCGA-转录因子分析、TCGA-ceRNA调控网络分析
8.其他,二代测序转录组数据自主分析、NCBI数据上传、二代测序数据解读、
如果觉得我的文章对您有用,请随意打赏。你的支持将鼓励我继续创作!