view NEUMA-1.2.1/merge_readcount.pl @ 0:c44c43d185ef draft default tip

NEUMA-1.2.1 Uploaded
author chawhwa
date Thu, 08 Aug 2013 00:46:13 -0400
parents
children
line wrap: on
line source

#!/usr/bin/perl

if(@ARGV<2) { print "usage: $0 type[gNIR/iNIR/gReadcount] readcount_dir > output.gNIR(iNIR/gReadcount).merged\n"; exit; }
my $suffix = shift @ARGV;  #either g or i
my $readcount_dir = shift @ARGV;

my @sample_list = split/\n/,`ls -1 $readcount_dir`;
for my $i (0..$#sample_list){
 if($sample_list[$i] =~ /\.$suffix$/) { $samples{$`}=1; }
}

@sample_list = keys %samples;


for my $sample (@sample_list){
  next if($sample eq '');
  my $file =  "$readcount_dir/$sample.$suffix";

  print STDERR "$file\n";  #DEBUGGING
  open IN,$file or die "Can't open $suffix file $file\n";
  <IN>;
  while(<IN>){
    chomp;

    my $gene;
    if($suffix eq 'gNIR' || $suffix eq 'gReadcount'){
      my ($gid,$expr) = (split/\t/)[0,1];
      $gene = "$gid";
      $EXPR{$gene}{$sample} = $expr;
    }
    elsif($suffix eq 'iNIR'){
      my ($gid,$NM,$expr) = (split/\t/)[0,1,2];
      #print STDERR "$gid\t$NM\t$expr\n";
      $gene = "$gid\t$NM";
      $EXPR{$gene}{$sample} = $expr;
    }
  }
  close IN;
}



$"="\t";
if($suffix eq 'gNIR' || $suffix eq 'gReadcount') { print "gene.id\t@sample_list\n"; }
elsif($suffix eq 'iNIR') { print "gene.id\tisoform\t@sample_list\n"; }
for my $gene (keys %EXPR){
 print "$gene";
 for my $sample (@sample_list){
   if(exists ${$EXPR{$gene}}{$sample}){
     print "\t$EXPR{$gene}{$sample}";
   }
   else { print "\t0"; }
 }
 print "\n";
}