这个需求是在生信分析中几乎天天用到,各种语言都能实现,也都各有特点。这次以perl 为例。
已知
文件CT-VS-CON.All.xls 为全部蛋白表达矩阵及其差异分析结果。
文件Homo_sapiens.ko 为蛋白KEGG注释结果。
文件Homo_sapiens.fa 为蛋白鉴定数据库(有的序列以多行展示)。
需求
将以上三表整合为一个表,输出需要的信息。
实现
#! /usr/bin/perl -w
use strict;
=pod
Description: combine table
Author:
Created:
Version:
=cut
use Getopt::Long;
#use Bio::SeqIO;
my ($exp,$ko,$fa,$help,$outdir);
GetOptions(
"exp:s" => \$exp,
"ko:s" => \$ko,
"fa:s" => \$fa,
"outdir:s" => \$outdir,
"help|?" => \$help
);
$outdir ||= ".";
if (!defined $exp || !defined $ko || !defined $fa || defined $help){
die << "USAGE";
description: combine table
usage: perl $0 [options]
options:
-exp <file> * all proteins expression matrix
-ko <file> * all proteins ko annotation table (protein=>ko)
-fa <file> all proteins sequence
-outdir <path> output directory, default is current directory "."
-help|? help information
eg:
perl $0 -exp a-VS-b.All.xls -ko All.ko -fa All.fa -outdir .
USAGE
}
my %ko;
open KO, "<$ko" or die $!;
while(<KO>){
chomp;
if($_ =~ /^#/){next;}
my @K=split/\t/;
if(defined $K[1]){
$ko{$K[0]}=$K[1];
}else{
$ko{$K[0]}="-";
}
}
close KO;
my %hash;
my $keys;
open DB, "< $fa" or die $!;
while(<DB>){
chomp;
if($_ =~ /^>(.*?)\s/){ #非贪婪匹配
$keys = $1;
}else{
$hash{$keys} .= $_;
}
}
close DB;
#foreach my $id(keys %hash){
# print "$id\n$hash{$id}\n";
#}
my %exp;
my ($ratio,$class,$des,$pvalue);
open OUT, ">out.txt" or die $!;
open EXP, "< $exp" or die $!;
while(<EXP>){
chomp;
my @tmp=split/\t/,$_;
if($.==1){
for(my $i=0;$i<=$#tmp;$i++){
if($tmp[$i] =~ /ratio/i){$ratio=$i;}
if($tmp[$i] =~ /pvalue/i){$pvalue=$i;}
if($tmp[$i] =~ /class/i){$class=$i;}
if($tmp[$i] =~ /Description/i){$des=$i;}
}
print OUT "Protein\tratio\tpvalue\tClass\tDescription\tlog2FC\tKEGG\tSequence\n";
next;
}
$tmp[$class] =~ s/Non/None/;
$exp{$tmp[0]}="$tmp[0]\t$tmp[$ratio]\t$tmp[$pvalue]\t$tmp[$class]\t$tmp[$des]\t";
$exp{$tmp[0]} .= log2($tmp[$ratio]);
}
close EXP;
my ($kegg,$seq);
foreach my $id (keys %exp){
if(exists $ko{$id}){
$kegg=$ko{$id};
}else{
$kegg="-";
}
if(exists $hash{$id}){
$seq=$hash{$id};
}else{$seq="-";}
print OUT "$exp{$id}\t$kegg\t$seq\n";
}
close OUT;
############################
# subroutine
############################
sub log2 {
my $n = shift;
return log($n)/log(2);
}
Perl 编写的特点是only write ,条条大道通罗马,高手也许几行就能解决,我这里写得有点冗余,但有些是不能省的,目的是为了更规范化,便于维护和他人阅读。
|
请发表评论