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

NEUMA-1.2.1 Uploaded

#!/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";
}