• 设为首页
  • 点击收藏
  • 手机版
    手机扫一扫访问
    迪恩网络手机版
  • 关注官方公众号
    微信扫一扫关注
    公众号

【Perl示例】整合多个文件

原作者: [db:作者] 来自: [db:来源] 收藏 邀请

这个需求是在生信分析中几乎天天用到,各种语言都能实现,也都各有特点。这次以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,条条大道通罗马,高手也许几行就能解决,我这里写得有点冗余,但有些是不能省的,目的是为了更规范化,便于维护和他人阅读。


鲜花

握手

雷人

路过

鸡蛋
该文章已有0人参与评论

请发表评论

全部评论

专题导读
上一篇:
perl在win下输出中文乱码问题发布时间:2022-07-22
下一篇:
perl内置特殊变量查询发布时间:2022-07-22
热门推荐
热门话题
阅读排行榜

扫描微信二维码

查看手机版网站

随时了解更新最新资讯

139-2527-9053

在线客服(服务时间 9:00~18:00)

在线QQ客服
地址:深圳市南山区西丽大学城创智工业园
电邮:jeky_zhao#qq.com
移动电话:139-2527-9053

Powered by 互联科技 X3.4© 2001-2213 极客世界.|Sitemap