<commit_before>#!/usr/bin/env perl

#########################################################################################
#                                       stepDESeq2.pl
#########################################################################################
# 
#  This program removes adapter sequence. 
#
#
#########################################################################################
# AUTHORS:
#
# <NAME>, PhD 
# Jul 4, 2014
#########################################################################################


############## LIBRARIES AND PRAGMAS ################

 use List::Util qw[min max];
 use strict;
 use File::Basename;
 use Getopt::Long;
 use Pod::Usage; 

#################### VARIABLES ######################
 my $cols             = "";
 my $conds            = "";
 my $num              = "";
 my $heatmap          = "";
 my $fitType          = "";
 my $foldChange       = "";
 my $dataset          = "";
 my $padj             = "";
 my $rscriptCMD       = "";
 my $pubdir           = "";
 my $wkey             = "";
 my $outdir           = "";
 my $jobsubmit        = "";
 my $servicename      = "";
 my $help             = "";
 my $print_version    = "";
 my $version          = "1.0.0";
################### PARAMETER PARSING ####################

my $cmd=$0." ".join(" ",@ARGV); ####command line copy

GetOptions( 
    'cols=s'         => \$cols,
    'dconds=s'       => \$conds,
    'dataset=s'      => \$dataset,
    'num=s'          => \$num,
    'eheatmap=s'     => \$heatmap,
    'tfitType=s'     => \$fitType,
    'padj=s'         => \$padj,
    'foldChange=s'   => \$foldChange,
    'rscriptCMD=s'   => \$rscriptCMD,
    'outdir=s'       => \$outdir,
    'pubdir=s'       => \$pubdir,
    'wkey=s'         => \$wkey,
    'servicename=s'  => \$servicename,
    'jobsubmit=s'    => \$jobsubmit,
    'help'           => \$help, 
    'version'        => \$print_version,
) or die("Unrecognized optioins.\nFor help, run this script with -help option.\n");

if($help){
    pod2usage( {
		'-verbose' => 2, 
		'-exitval' => 1,
	} );
}

if($print_version){
  print "Version ".$version."\n";
  exit;
}

pod2usage( {'-verbose' => 0, '-exitval' => 1,} ) if ( ($cols eq "") or ($conds eq "") or ($outdir eq "") );	

################### MAIN PROGRAM ####################
#    maps the reads to the ribosome and put the files under $outdir/after_ribosome directory

my $inputdir = "$outdir/rsem";
my $input_file_suffix = "_expression_expected_count.tsv";
if (lc($dataset)!~/rsem/)
{
   $inputdir="$outdir/counts";
   $input_file_suffix = ".counts.tsv";
}

$outdir   = "$outdir/DESeq2".$dataset.$num;
`mkdir -p $outdir`;
$cols=~s/[\s\t]+//g;
$conds=~s/[\s\t]+//g;
$cols = checkCols($cols);
$cols=~s/,/\",\"/g;
$cols="c(\"$cols\")";
$cols=~s/-/\./g;
$conds=~s/,/\",\"/g;
$conds="c(\"$conds\")";

my $puboutdir   = "$pubdir/$wkey";
`mkdir -p $puboutdir`;
if (lc($dataset) =~/rsem/ )
{
   makePlot( "genes", $input_file_suffix, $inputdir, $outdir, $cols, $conds, $fitType, $heatmap, $padj, $foldChange, $puboutdir, "DESeq2".$dataset.$num, $wkey);
   makePlot( "isoforms", $input_file_suffix, $inputdir, $outdir, $cols, $conds, $fitType, $heatmap, $padj, $foldChange, $puboutdir, "DESeq2".$dataset.$num, $wkey);
}
else
{
   makePlot( $dataset, $input_file_suffix, $inputdir, $outdir, $cols, $conds, $fitType, $heatmap, $padj, $foldChange, $puboutdir, "DESeq2".$dataset.$num, $wkey);
}

`cp -R $outdir $puboutdir/.`;

sub makePlot
{
my ($type,$input_file_suffix,$inputdir, $outdir, $cols, $conds, $fitType, $heatmap, $padj_cutoff, $foldChange_cutoff, $puboutdir, $deseqdir, $wkey)=@_;
my $output = "$outdir/rscript_$type.R";
my $table = "$outdir/deseq2_$type.tsv";
my $pdfname = "$outdir/heatmap_$type.pdf";
my $alldetected = "$outdir/alldetected_$type.tsv";
my $selected_log2fc = "$outdir/selected_log2fc_$type.tsv";
my $sessioninfo = "$outdir/sessionInfo.txt";
my $inputfile=$inputdir."/".$type.$input_file_suffix;
my $col=1;
$col=2 if ($type eq "isoforms");

my $heatmapR="";
if (lc($heatmap) eq "yes")
{
$heatmapR=qq/
  f1<-res[!is.na(res\$padj) & !is.na(res\$log2FoldChange), ]
  res_selected<-f1[(f1\$padj<$padj_cutoff & abs(f1\$log2FoldChange)>log2($foldChange_cutoff)),]
  ld <- log(filtd[rownames(res_selected),]+0.1,base=2)

  cldt <- scale(t(ld), center=TRUE, scale=TRUE);
  cld <- t(cldt)

  dissimilarity <- 1 - cor(cld)
  distance <- as.dist(dissimilarity)
  pdf("$pdfname")
  plot(hclust(distance),  main="Dissimilarity = 1 - Correlation", xlab="")
  heatmap.2(cld, col=redgreen(75), scale="row",
          key=TRUE, symkey=FALSE, density.info="none",cexRow=1,cexCol=1,margins=c(12,8),trace="none",srtCol=45)
  dev.off()
/;
}

open(OUT, ">$output");
my $rscript=qq/
library("DESeq2")<commit_msg>DESeq2 fix<commit_after>#!/usr/bin/env perl

#########################################################################################
#                                       stepDESeq2.pl
#########################################################################################
# 
#  This program removes adapter sequence. 
#
#
#########################################################################################
# AUTHORS:
#
# <NAME>, PhD 
# Jul 4, 2014
#########################################################################################


############## LIBRARIES AND PRAGMAS ################

 use List::Util qw[min max];
 use strict;
 use File::Basename;
 use Getopt::Long;
 use Pod::Usage; 

#################### VARIABLES ######################
 my $cols             = "";
 my $conds            = "";
 my $num              = "";
 my $heatmap          = "";
 my $fitType          = "";
 my $foldChange       = "";
 my $dataset          = "";
 my $padj             = "";
 my $rscriptCMD       = "";
 my $pubdir           = "";
 my $wkey             = "";
 my $outdir           = "";
 my $jobsubmit        = "";
 my $servicename      = "";
 my $help             = "";
 my $print_version    = "";
 my $version          = "1.0.0";
################### PARAMETER PARSING ####################

my $cmd=$0." ".join(" ",@ARGV); ####command line copy

GetOptions( 
    'cols=s'         => \$cols,
    'dconds=s'       => \$conds,
    'dataset=s'      => \$dataset,
    'num=s'          => \$num,
    'eheatmap=s'     => \$heatmap,
    'tfitType=s'     => \$fitType,
    'padj=s'         => \$padj,
    'foldChange=s'   => \$foldChange,
    'rscriptCMD=s'   => \$rscriptCMD,
    'outdir=s'       => \$outdir,
    'pubdir=s'       => \$pubdir,
    'wkey=s'         => \$wkey,
    'servicename=s'  => \$servicename,
    'jobsubmit=s'    => \$jobsubmit,
    'help'           => \$help, 
    'version'        => \$print_version,
) or die("Unrecognized optioins.\nFor help, run this script with -help option.\n");

if($help){
    pod2usage( {
		'-verbose' => 2, 
		'-exitval' => 1,
	} );
}

if($print_version){
  print "Version ".$version."\n";
  exit;
}

pod2usage( {'-verbose' => 0, '-exitval' => 1,} ) if ( ($cols eq "") or ($conds eq "") or ($outdir eq "") );	

################### MAIN PROGRAM ####################
#    maps the reads to the ribosome and put the files under $outdir/after_ribosome directory

my $inputdir = "$outdir/rsem";
my $input_file_suffix = "_expression_expected_count.tsv";
if (lc($dataset)!~/rsem/)
{
   $inputdir="$outdir/counts";
   $input_file_suffix = ".counts.tsv";
}

$outdir   = "$outdir/DESeq2".$dataset.$num;
`mkdir -p $outdir`;
$cols=~s/[\s\t]+//g;
$conds=~s/[\s\t]+//g;
$cols = checkCols($cols);
$cols=~s/,/\",\"/g;
$cols="c(\"$cols\")";
$cols=~s/-/\./g;
$conds=~s/,/\",\"/g;
$conds="c(\"$conds\")";

my $puboutdir   = "$pubdir/$wkey";
`mkdir -p $puboutdir`;
if (lc($dataset) =~/rsem/ )
{
   makePlot( "genes", $input_file_suffix, $inputdir, $outdir, $cols, $conds, $fitType, $heatmap, $padj, $foldChange, $puboutdir, "DESeq2".$dataset.$num, $wkey);
   makePlot( "isoforms", $input_file_suffix, $inputdir, $outdir, $cols, $conds, $fitType, $heatmap, $padj, $foldChange, $puboutdir, "DESeq2".$dataset.$num, $wkey);
}
else
{
   makePlot( $dataset, $input_file_suffix, $inputdir, $outdir, $cols, $conds, $fitType, $heatmap, $padj, $foldChange, $puboutdir, "DESeq2".$dataset.$num, $wkey);
}

`cp -R $outdir $puboutdir/.`;

sub makePlot
{
my ($type,$input_file_suffix,$inputdir, $outdir, $cols, $conds, $fitType, $heatmap, $padj_cutoff, $foldChange_cutoff, $puboutdir, $deseqdir, $wkey)=@_;
my $output = "$outdir/rscript_$type.R";
my $table = "$outdir/deseq2_$type.tsv";
my $pdfname = "$outdir/heatmap_$type.pdf";
my $alldetected = "$outdir/alldetected_$type.tsv";
my $selected_log2fc = "$outdir/selected_log2fc_$type.tsv";
my $sessioninfo = "$outdir/sessionInfo.txt";
my $inputfile=$inputdir."/".$type.$input_file_suffix;
my $col=1;
$col=2 if ($type eq "isoforms");

my $heatmapR="";
if (lc($heatmap) eq "yes")
{
$heatmapR=qq/
  f1<-res[!is.na(res\$padj) & !is.na(res\$log2FoldChange), ]
  res_selected<-f1[(f1\$padj<$padj_cutoff & abs(f1\$log2FoldChange)>log2($foldChange_cutoff)),]
  if (dim(res_selected)>4 && length(conds)>2)
  {
  ld <- log(filtd[rownames(res_selected),]+0.1,base=2)

  cldt <- scale(t(ld), center=TRUE, scale=TRUE);
  cld <- t(cldt)

  dissimilarity <- 1 - cor(cld)
  distance <- as.dist(dissimilarity)
  pdf("$pdfname")
  plot(hclust(distance),  main="Dissimilarity = 1 - Correlation", xlab="")
  heatmap.2(cld, col=redgreen(75), scale="row",
          key=TRUE, symkey=FALSE, density.info="none",cexRow=1,cexCol=1,margins=c(12,8),trace="none",srtCol=45)
  dev.off()
  }
/;
}

open(OUT, ">$output");
my $rscript=qq/
library("DESeq2")
{'fraction_non_alphanumeric': 0.2981444595869588, 'fraction_numerical': 0.010378446378027047, 'mean_word_length': 3.5689655172413794, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 51, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5193808', 'n_tokens_mistral': 3719, 'n_tokens_neox': 3506, 'n_words': 645}