美文网首页生物信息编程
Perl高级语法 多维哈希嵌套

Perl高级语法 多维哈希嵌套

作者: 214b3ff96d82 | 来源:发表于2019-08-18 10:36 被阅读124次

最近在使用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; 

始终这一行提示报错,我也不知道怎么解决
但是可以正常输出结果文件
所以,我想知道怎么才能不报错

image.png

相关文章

网友评论

    本文标题:Perl高级语法 多维哈希嵌套

    本文链接:https://www.haomeiwen.com/subject/bcfrsctx.html