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