File: merge_rsem_n_express_for_compare.pl

package info (click to toggle)
trinityrnaseq 2.11.0%2Bdfsg-6
  • links: PTS, VCS
  • area: main
  • in suites: bullseye
  • size: 417,528 kB
  • sloc: perl: 48,420; cpp: 17,749; java: 12,695; python: 3,124; sh: 1,030; ansic: 983; makefile: 688; xml: 62
file content (87 lines) | stat: -rwxr-xr-x 1,904 bytes parent folder | download | duplicates (4)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
#!/usr/bin/env perl

use strict;
use warnings;

my $usage = "usage: $0 rsem.file eXpress.file count|FPKM\n\n";

my $rsem_file = $ARGV[0] or die $usage;
my $express_file = $ARGV[1] or die $usage;
my $method = $ARGV[2] or die $usage;

unless ($method =~ /^(count|FPKM)$/) { 
    die $usage;
}

main: {

    my %data;
    &add_RSEM_data(\%data, $rsem_file, $method);
    &add_eXpress_data(\%data, $express_file, $method);


    print join("\t", "", "RSEM", "eXpress") . "\n";
    foreach my $id (keys %data) {
        my $rsem_val = $data{$id}->{RSEM};
        unless (defined $rsem_val) {
            $rsem_val = -1;
        }
        my $express_val = $data{$id}->{eXpress};
        unless (defined $express_val) {
            $express_val = -1;
        }

        if ($rsem_val > 0 && $express_val > 0) {
            $rsem_val = "NA" if $rsem_val < 0;
            $express_val = "NA" if $express_val < 0;
            
            print join("\t", $id, $rsem_val, $express_val) . "\n";
        }
    }
    
    exit(0);
}

####
sub add_RSEM_data {
    my ($data_href, $rsem_file, $method) = @_;

    my $data_field = ($method eq 'count') ? 4 : 6;
   
    open (my $fh, $rsem_file) or die "Error, cannot open file $rsem_file";
    my $header = <$fh>;
    while (<$fh>) {
        chomp;
        my @x = split(/\t/);
        my $id = $x[0];
        my $val = $x[$data_field];

        $data_href->{$id}->{RSEM} = $val;
    }
    close $fh;

    return;
}

####
sub add_eXpress_data {
    my ($data_href, $express_file, $method) = @_;

    my $data_field = ($method eq 'count') ? 7 : 10;

    open (my $fh, $express_file) or die "Error, cannot open file $express_file";
    my $header = <$fh>;
    while (<$fh>) {
        chomp;
        my @x = split(/\t/);
        my $id = $x[1];
        my $val = $x[$data_field];

        $data_href->{$id}->{eXpress} = $val;
    }
    close $fh;

    return;
}