diff NEUMA-1.2.1/merge_LVKM_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 diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/NEUMA-1.2.1/merge_LVKM_readcount.pl	Thu Aug 08 00:46:13 2013 -0400
@@ -0,0 +1,82 @@
+#!/usr/bin/perl
+
+if(@ARGV<4) { print "usage: $0 genewise/isoformwise[g/i] EUMAcut LVKM_out_dir NR(1/0) > output.gNIR(iNIR)\n"; exit; }
+my $option = shift @ARGV;  #either g or i
+my $EUMAcut = shift @ARGV;
+my $LVKM_dir = shift @ARGV;
+my $NR_option = shift @ARGV;
+
+my $suffix = $option."LVKM";
+
+my @sample_list = split/\n/,`ls -1 $LVKM_dir`;
+for my $i (0..$#sample_list){
+ if($NR_option eq '1'){
+   if($sample_list[$i] =~ /\.\d+.-NR.[gi]LVKM$/) { $samples{$`}=1; }
+ }
+ else {
+   if($sample_list[$i] =~ /\.\d+.[gi]LVKM$/) { $samples{$`}=1; }
+ }
+}
+
+@sample_list = keys %samples;
+
+
+for my $sample (@sample_list){
+  next if($sample eq '');
+  my $file;
+  if($NR_option eq '1'){
+    $file =  "$LVKM_dir/$sample.$EUMAcut.-NR.$suffix";
+  }
+  else {
+    $file =  "$LVKM_dir/$sample.$EUMAcut.$suffix";
+  }
+
+  print STDERR "$file\n";  #DEBUGGING
+  open IN,$file or die "Can't open LVKM file $file\n";
+  <IN>;
+  while(<IN>){
+    chomp;
+
+    my $gene;
+    if($option eq 'g'){
+      my ($gid,$symbol,$origGFVK,$gEUMA) = (split/\t/)[0,1,5,7];
+      $gene = "$gid\t$symbol";
+      $temp = $origGFVK * $gEUMA;
+      if($temp == 0){
+        $EXPR{$gene}{$sample} = 0;
+      }else{
+        $EXPR{$gene}{$sample} = int(($temp + 0.5) /1000);
+      }
+    }
+    elsif($option eq 'i'){
+      my ($gid,$symbol,$NM,$finalIFVK,$iEUMA) = (split/\t/)[0,1,2,5,8];
+      $gene = "$gid\t$symbol\t$NM";
+      $temp = $finalIFVK * $iEUMA;
+      if($temp == 0){
+        $EXPR{$gene}{$sample} = 0;
+      }else{
+        $EXPR{$gene}{$sample} = int(($temp + 0.5) /1000);
+      }
+    }
+  }
+  close IN;
+}
+
+
+
+$"="\t";
+if($option eq 'g') { print "gene.id\tgene.symbol\t@sample_list\n"; }
+elsif($option eq 'i') { print "gene.id\tgene.symbol\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 "\tNA"; }
+ }
+ print "\n";
+}
+
+
+