美文网首页Perl小推车进化生物信息编程
Perl编程练习—双重哈希、哈希数组

Perl编程练习—双重哈希、哈希数组

作者: Dawn_WangTP | 来源:发表于2018-08-15 22:42 被阅读63次

    前几天晚上在练习Perl编程时突然顿悟了以前一直迷惑的哈希数组。对这一类的问题的解决瞬间开拓了新思路。这两天利用双重哈希/哈希数组解决了一些数据转换中的一些小问题。记录分享一下。

    利用哈希数组整合相同的数据名

    • 需求:在基因ID转换过程中,存在着一个GeneID对应着多个SYMBOL NAME,和description信息,需要对这多个SYMBOL做个整合,最终得到一个ID对一个SYMBOL对一个DESCRIPTION。
    • 利用哈希数组中的哈希键是唯一的特性对值做整合


      输入文件格式
      输出文件格式
    #!/usr/bin/env perl
    use strict;
    
    my %hash;
    open IN,$ARGV[0] or die "$!";
    while(<IN>){
        chomp;
        my @F = split/\t/;
        push @{$hash{$F[1]}{$F[3]}},$F[2];
    }
    foreach my $id (sort keys %hash){
        foreach my $description (sort keys %{$hash{$id}} ){
            my $name=join "//",@{$hash{$id}{$description}};
            print "$id\t$name\t$description\n";
        }
    }
    
    

    gff文件和fasta文件提取CDS序列

    • 根据gff文件中的CDS坐标提取各个转录本的核酸序列,注意存在转录本的正负链问题,每个转录本的第一个cds有frame移码问题
    • 利用双重哈希,以及程序USAGE基本用法。


      GFF文件
      FASTA
    #!/usr/bin/env perl
    use strict;
    
    my $usage=<<USAGE;
    Usage:
        perl $0 genome.fasta genome.gff3
        
    USAGE
    print $usage if(@ARGV==0);
    
    open IN ,"$ARGV[0]" or die "$!";
    my (%seq,$seq_id);
    while(<IN>){
        chomp;
        if(/^>(\S+)/){$seq_id=$1}else{
            $seq{$seq_id}.=$_;
        }
    }
    close IN;
    
    open IN,"$ARGV[1]" or die "$!";
    my(%cds_pos, %contig, %stream );
    while(<IN>){
        chomp;
        if(/\tCDS\t/){
            my @F=split/\t/;
            $F[8]=~/Parent=(.*?);/;
            $cds_pos{$1}{"$F[3]\t$F[4]\t$F[7]"}=1;
            $contig{$1}=$F[0];
            $stream{$1}=$F[6];
        }
    }
    close IN;
    
    foreach my $id(sort keys %contig){
        my $cds_seq;
        my @cds_pos_parsing=sort { $cds_pos{$id}{$a}<=>$cds_pos{$id}{$b} } keys %{$cds_pos{$id}};
        foreach (@cds_pos_parsing){
            my @F=split/\t/;
            $cds_seq.=substr($seq{$contig{$id}},$F[0]-1,$F[1]-$F[0]+1);
        }
        if($stream{$id} eq "+"){
            my $frame = $1 if ($cds_pos_parsing[0] =~/(\d+)$/);
            $cds_seq=~s/\w{$frame}//;
            print ">$id\n$cds_seq\n";
        }else{
            my $frame = $1 if ($cds_pos_parsing[-1]=~/(\d+)$/);
            $cds_seq = reverse($cds_seq);
            $cds_seq =~ tr/ATGCatgcn/TACGtacgN/;
            $cds_seq =~ s/\w{$frame}//;
            print">$id\n$cds_seq\n";
        }
    }
    
    
    

    相关文章

      网友评论

        本文标题:Perl编程练习—双重哈希、哈希数组

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