千家信息网

perl如何输出基因的位置信息

发表于:2025-01-19 作者:千家信息网编辑
千家信息网最后更新 2025年01月19日,这篇文章主要介绍了perl如何输出基因的位置信息的相关知识,内容详细易懂,操作简单快捷,具有一定借鉴价值,相信大家阅读完这篇perl如何输出基因的位置信息文章都会有所收获,下面我们一起来看看吧。per
千家信息网最后更新 2025年01月19日perl如何输出基因的位置信息

这篇文章主要介绍了perl如何输出基因的位置信息的相关知识,内容详细易懂,操作简单快捷,具有一定借鉴价值,相信大家阅读完这篇perl如何输出基因的位置信息文章都会有所收获,下面我们一起来看看吧。

perl输出基因的位置信息按照基因所在染色体,和位置信息排序

我们在整理基因组的gff文件,想输出基因的位置信息,以及基因所对应的多个转录本信息,需要对基因按照染色体排序,这里使用到了perl里面hash按照值来排序,而且还用了两个值基因型排序。示例代码如下,以下代码可以从gff文件当中提取所有基因的位置信息以及对应的多个转录本信息:

perl代码如下:

#!/usr/bin/perl -wuse strict;use Cwd qw(abs_path getcwd);use Getopt::Long;use Data::Dumper;die "perl $0  " unless(@ARGV==2);my$gff=$ARGV[0];my%gene=();my%gene_region=();my%mRNA2Gene=();my%Gene2mRNA=();open IN,"$gff" or die "$!";open OUT ,">$ARGV[1]" or die "$!";print OUT "#gene_ID\tchr\tstart\tend\tstrand\ttranscript_id\n";while(){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}){push @{$Gene2mRNA{$pid}},$id;}#print "mRNA:$id\n";}}close(IN);#多层排序,先按染色体排序,再按基因位置排序for my$id(sort {$gene_region{$a}->[0] cmp $gene_region{$b}->[0] or $gene_region{$a}->[1] <=> $gene_region{$b}->[1] } keys %gene_region){print OUT "$id\t".join("\t",(@{$gene_region{$id}},sort @{$Gene2mRNA{$id}})    )."\n";}close(OUT);

关于"perl如何输出基因的位置信息"这篇文章的内容就介绍到这里,感谢各位的阅读!相信大家对"perl如何输出基因的位置信息"知识都有一定的了解,大家如果还想学习更多知识,欢迎关注行业资讯频道。

0