前几天晚上在练习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";
    }
}