##################################################################################
#Andy Rampersaud
#05.26.2015
#06.30.2015: Tisha edited to include edgeR results
#This R script is used to parse the DESeq and EdgeR output files
#Specifically the *_forSEGEXUpload.txt generated by DESeq_1* jobs
#The *_forSEGEXUpload*.txt is parsed to determine the number of differential genes
#Standard fold change and p-value cutoffs are used but they can be user-defined values
#08.04.15: Andy updated to include counting method in output file name
#----------------------------------------------------------------------------------
#Notes:
#<DiffExp_Output> 			: the *_forSEGEXUpload.txt file
#<fold_change_cutoff> 			: the |fold change| minimum for differential genes
#<padj_cutoff> 				: the padj maximum for differential genes
#<differential_expression_method>	: the R package used (DESeq or EdgeR)
#<count_method>				: the counting method (GeneBody, Exonic_Only, Intronic_Only) (Use the COL_SUFFIX variable from the DiffExp_* job)
#----------------------------------------------------------------------------------
#Usage: 
#Rscript Diff_Genes.R <DiffExp_Output> <fold_change_cutoff> <padj_cutoff> <differential_expression_method> <count_method>	
#Example command to run script:
#Rscript Diff_Genes.R Male_G83_M1M2_vs_Female_G83_M3M4_DiffExp_v2_GeneBody_forSEGEXUpload_DESeq.txt 2 0.05 DESeq GeneBody
#OR
#Rscript Diff_Genes.R Male_G83_M1M2_vs_Female_G83_M3M4_DiffExp_v2_GeneBody_forSEGEXUpload_EdgeR.txt 2 0.05 EdgeR GeneBody
##################################################################################
#---------------------------------------------------------------------------
#Pass arguments from sh into R
args <- commandArgs(trailingOnly = TRUE)
DiffExp_Output <- args[1]
fold_change_cutoff <- as.numeric(args[2])
padj_cutoff <- as.numeric(args[3])
postfix <- args[4]
count_method <- args[5]
#---------------------------------------------------------------------------
if(length(args) != 5){
print("Need 5 arguments:")
print("Usage: Rscript Diff_Genes.R <DiffExp_Output> <fold_change_cutoff> <padj_cutoff> <differential_expression_method> <count_method>	")
}
#---------------------------------------------------------------------------
print("Print arguments:")
print("-----------------")
print("Differential expression_Output:")
paste(DiffExp_Output,sep="")
print("fold_change_cutoff:")
paste(fold_change_cutoff,sep="")
print("padj_cutoff:")
paste(padj_cutoff,sep="")
print("postfix:")
paste(postfix,sep="")
print("count_method:")
paste(count_method,sep="")
print("-----------------")
#---------------------------------------------------------------------------
#Use a pattern match in an "if" statement 
#to remind user to use the  *.annotated file
#Spaces/brackets matter for "if" statement syntax
#Need invisible() to avoid printing NULL
invisible(if ( grepl("_forSEGEXUpload_DESeq.txt", DiffExp_Output)  || grepl("_forSEGEXUpload_EdgeR.txt", DiffExp_Output))  {
#Do nothing
} else {
print("WARNING: DiffExp_Output is not the *_forSEGEXUpload_DESeq.txt or *_forSEGEXUpload_EdgeR.txt file!")
}#End of else statement
)#End of invisible
#---------------------------------------------------------------------------
#Need to set the dir so that this script can work in any folder:
dir <- getwd()
setwd(dir)
#---------------------------------------------------------------------------
#Check to make sure *_forSEGEXUpload_DESeq.txt file exists
#Need invisible() to avoid printing NULL
invisible(if ( length(list.files(pattern = "_forSEGEXUpload_DESeq.txt")) >= 1 || length(list.files(pattern = "_forSEGEXUpload_EdgeR.txt")) >= 1)  {
#Do nothing
} else {
print("DESeq or EdgeR job most likely failed.")
print("Quitting R now.")
quit()
}#End of else statement
)#End of invisible
#---------------------------------------------------------------------------
DiffExp_data <- read.table(file=DiffExp_Output, sep="\t", header=TRUE, as.is=TRUE, fill = TRUE)
#colnames(DESeq_data)
#[1] "id"               "ratio"            "foldChange"      
#[4] "rpkm_mean_Female" "rpkm_mean_Male"   "padj"            
#[7] "pval"      
#dim(DESeq_data)
#[1] 24197     7
#View(DESeq_data)
############################################################
print("Differential Gene Counts")
print("Number of significant differential genes with positive fold change (Up Genes):")
#Make sure cutoff values are numeric:
fold_change_cutoff <- as.numeric(fold_change_cutoff)
padj_cutoff <- as.numeric(padj_cutoff)
Up_Genes <- DiffExp_data[DiffExp_data$"foldChange" > fold_change_cutoff & DiffExp_data$"padj" < padj_cutoff,]
dim(Up_Genes)
#[1] 121   7
#G83 Test data:
#[1] 33  7
#View(Up_Genes_DESeq)
print("Number of significant differential genes with negative fold change (Down Genes):")
#Need the negative fold change:
fold_change_cutoff_neg <- as.numeric(paste("-",fold_change_cutoff, sep=""))
Down_Genes <- DiffExp_data[DiffExp_data$"foldChange" < fold_change_cutoff_neg & DiffExp_data$"padj" < padj_cutoff,]
dim(Down_Genes)
#[1] 169   7
#G83 Test data:
#[1] 33  7
#View(Down_Genes_DESeq)
############################################################
#---------------------------------------------------------------------------
#Sort data_frames
Up_Genes <- Up_Genes[order(-Up_Genes$"foldChange",Up_Genes$"padj"),]
Down_Genes <- Down_Genes[order(Down_Genes$"foldChange",Down_Genes$"padj"),]
#---------------------------------------------------------------------------
#Reformat the data frames for readability:
#---------------------------------------------------------------------------
#Need p-values as numeric:
Up_Genes$pval <- as.numeric(as.character(Up_Genes$pval))
Up_Genes$padj <- as.numeric(as.character(Up_Genes$padj))
#Format the p-values separately:
Up_Genes$pval <- format(Up_Genes$pval, digits=2, nsmall=2)
Up_Genes$padj <- format(Up_Genes$padj, digits=2, nsmall=2)
#---------------------------------------------------------------------------
#Need p-values as numeric:
Down_Genes$pval <- as.numeric(as.character(Down_Genes$pval))
Down_Genes$padj <- as.numeric(as.character(Down_Genes$padj))
#Format the p-values separately:
Down_Genes$pval <- format(Down_Genes$pval, digits=2, nsmall=2)
Down_Genes$padj <- format(Down_Genes$padj, digits=2, nsmall=2)
#---------------------------------------------------------------------------
Up_Genes <- format(Up_Genes, digits=2, scientific=FALSE)
Down_Genes <- format(Down_Genes, digits=2, scientific=FALSE)
#---------------------------------------------------------------------------
#Make text files of differential genes
write.table(Up_Genes, file = paste("Up_Genes_", postfix, "_", count_method,".txt", sep =""), quote=FALSE, sep = "\t",row.names = FALSE,col.names = TRUE)
write.table(Down_Genes, file = paste("Down_Genes_", postfix, "_", count_method,".txt", sep =""), quote=FALSE, sep = "\t",row.names = FALSE,col.names = TRUE)
#---------------------------------------------------------------------------
#Save output files to a variable:
Output_File1 <- paste("Up_Genes_", postfix, "_", count_method,".txt", sep ="")
Output_File2 <- paste("Down_Genes_", postfix, "_", count_method,".txt", sep ="")
#---------------------------------------------------------------------------
print(paste("Check out ",Output_File1, "!", sep=""))
print(paste("Check out ",Output_File2, "!", sep=""))
##################################################################################