最近在使用perl操作基因组注释文件的时候,看到了某大神写的代码,里面涉及了二维哈希,以及哈希嵌套等perl高级语法。我是perl的初级水平,根本看不懂。所以,现在我把代码以及操作要求贴出来,希望有大神可以指导一下。
01. 输入文件格式
第一个文件格式
- 第一列为染色体编号
- 第二列为基因组位置
- 第三到3-6列代表每个物种在该位置的核苷酸碱基
第二个文件
- 就是一个标准的gff3文件
02. 要求输出文件
现在我就想根据基因组的注释文件GFF3,然后把第一个文件中属于CDS的序列提取出来
03. Perl脚本代码
#! /usr/bin/env perl
use strict;
use warnings;
my $gff="Bos_taurus.UMD3.1.87.gff3.check.gff";
my $lst=$ARGV[0];
my $dir="genes";
`mkdir $dir` if(!-e $dir);
my %gene;
my %pos;
open I,"< $gff";
while(<I>){
chomp;
my @a=split(/\s+/);
next unless($a[2] eq "CDS");
my ($chr,$start,$end,$strand)=($a[0],$a[3],$a[4],$a[6]);
$a[8]=~/Parent=transcript:(\w+)/;
my $id=$1;
for(my $i=$start;$i<=$end;$i++){
$pos{$chr}{$i}=$id;
$gene{$id}{pos}{$i}=0; ##这一行提示错误
}
$gene{$id}{strand}=$strand;
}
close I;
print STDERR "$gff loaded...\n";
my %seq;
open I,"cat $lst |";
my $head=<I>;
my @head=split(/\s+/,$head);
my $control=0;
while(<I>){
my @a=split(/\s+/);
my ($chr,$pos)=($a[0],$a[1]);
next unless(exists $pos{$chr}{$pos});
print STDERR "[ $control ] sites loaded...\r" if($control % 1000 == 0);
my $id=$pos{$chr}{$pos};
for(my $i=2;$i<@a;$i++){
my $species=$head[$i];
$a[$i]=uc($a[$i]);
$seq{$id}{$species}{$pos}=$a[$i];
}
$control++;
# last if($control++>10000);
}
close I;
print STDERR "$lst loaded...\n";
foreach my $id(sort keys %seq){
print STDERR "$id\r";
open O,"> $dir/$id.fa";
foreach my $species(sort keys %{$seq{$id}}){
my $nucl;
foreach my $pos(sort {$a<=>$b} keys %{$gene{$id}{pos}}){
my $base="-";
if(exists $seq{$id}{$species}{$pos}){
$base=$seq{$id}{$species}{$pos};
}
$nucl.=$base;
}
my $strand = $gene{$id}{strand};
if($strand eq "-"){
$nucl=~tr/ATCGRYMK/TAGCYRKM/;
$nucl=reverse($nucl);
}
print O ">$species\n$nucl\n";
}
close O;
}
print STDERR "\nDone\n";
04. perl 脚本报错
$gene{$id}{pos}{$i}=0;
始终这一行提示报错,我也不知道怎么解决
但是可以正常输出结果文件
所以,我想知道怎么才能不报错
网友评论