rsem-generate-data-matrix

   1 #!/usr/bin/env perl
   2
   3 use strict;
   4
   5 if (scalar(@ARGV) == 0) {
   6     print "Usage: rsem-generate-data-matrix sampleA.[alleles/genes/isoforms].results sampleB.[alleles/genes/isoforms].results ... > output_name.matrix\n";
   7     print "All result files should have the same file type. The 'TPM' columns of every result file are extracted to form the data matrix.\n";
   8     exit(-1);
   9 }
  10
  11 my $offsite = 4; # for new file formats
  12 if ($ARGV[0] =~ /alleles.results$/) { $offsite = 5; }
  13
  14 my $line;
  15 my $n = scalar(@ARGV);
  16 my $M = -1;
  17 my @matrix = ();
  18
  19 # 0, file_name; 1, reference of expected count array; 2, reference of transcript_id/gene_id array
  20 sub loadData {
  21     open(INPUT, $_[0]);
  22     my $line = <INPUT>; # The first line contains only column names
  23     while ($line = <INPUT>) {
  24         chomp($line);
  25         my @fields = split(/\t/, $line);
  26         push(@{$_[2]}, "\"$fields[0]\"");
  27         push(@{$_[1]}, $fields[$offsite]);
  28     }
  29     close(INPUT);
  30
  31     if (scalar(@{$_[1]}) == 0) {
  32         print STDERR "Nothing is detected! $_[0] may not exist or is empty.\n";
  33         exit(-1);
  34     }
  35 }
  36
  37 #0, M; 1, reference of @ids_arr; 2, reference of @ids
  38 sub check {
  39     my $size = $_[0];
  40     for (my $i = 0; $i < $size; $i++) {
  41         if ($_[1]->[$i] ne $_[2]->[$i]) {
  42             return 0;
  43         }
  44     }
  45     return 1;
  46 }
  47
  48 my @ids_arr = ();
  49
  50 for (my $i = 0; $i < $n; $i++) {
  51     my (@ids, @ecs) = ();
  52     &loadData($ARGV[$i], \@ecs, \@ids);
  53
  54     if ($M < 0) {
  55         $M = scalar(@ids);
  56         @ids_arr = @ids;
  57     }
  58     elsif (!&check($M, \@ids_arr, \@ids)) {
  59         print STDERR "Number of lines among samples are not equal!\n";
  60         exit(-1);
  61     }
  62
  63     my $colname;
  64     if (substr($ARGV[$i], 0, 2) eq "./") { $colname = substr($ARGV[$i], 2); }
  65     else { $colname = $ARGV[$i]; }
  66     $colname = "\"$colname\"";
  67     @ecs = ($colname, @ecs);
  68     push(@matrix, \@ecs);
  69 }
  70
  71 @ids_arr = ("", @ids_arr);
  72 @matrix = (\@ids_arr, @matrix);
  73
  74 for (my $i = 0; $i <= $M; $i++) {
  75     for (my $j = 0; $j < $n; $j++) { print "$matrix[$j][$i]\t"; }
  76     print "$matrix[$n][$i]\n";
  77 }