#!/usr/bin/env Rscript
###Copyright 2019 Ying Yu from Fudan-PGx group 
# example:
# Rscript RNAseq_4_pwDEG.R  -i ballgown_geneexp_log2fpkm_floor0p01_c3r58395_2019-04-29.txt  -g group1.txt -p organoid 

# choppy report script like : @scatter-plot(dataFile='/mnt/c/Users/YY/Documents/working/choppy_report/data/zhanggroup_P1-6vsP7-13_choppy_scatterplot_degs.rds', dataType='rds', xAxis='log2FC', xTitle="log2FC",yAxis='log10p',yTitle="-log10 (p)")

suppressPackageStartupMessages(library("optparse"))

# specify our desired options in a list
# by default OptionParser will add an help option equivalent to 
# make_option(c("-h", "--help"), action="store_true", default=FALSE, 
#               help="Show this help message and exit")

option_list <- list( 
    make_option(c("-o", "--out_dir"), type="character",default="./",
        help="The output directory [default ./]"),
    make_option(c("-i", "--input"),type="character", default=NULL,
        help="The input expression files. Required!"),
    make_option(c("-g", "--sample_group"),type="character",  default=NULL,
        help="File in tab-delimited format for sample group infomation. The input file containing sample name and group infomation. note colname must be like: sample	group1	group2... Required! "),
   	make_option(c("-p", "--project_code"), type="character",default="rnaseq",
        help="Project code, which is used as prefix of output file. [default: rnaseq]"),
	make_option(c("-a", "--output_all_genes"), metavar="FALSE", default=FALSE,
        help="Output rds files for choppy contains all genes. By default, only DEGs are listed in the output rds and csv for report. NOTE choppy report may not be availble to display correctly if too many points exit. [default: FALSE]"),
	make_option(c("-b", "--low_expr_filter"), metavar="FALSE",default=FALSE,
        help="Conduct low expression filtering before DEG analysis. [default: FALSE]"),
	make_option(c("-f", "--low_expr_filter_cutoff"), type="double",default=0,metavar="number",help="Genes across all samples with expreesion lower than this value will be filtered out [default: 0]")
		)

# get command line options, if help option encountered print help and exit,
# otherwise if options not found on command line then set defaults, 
opt <- parse_args(OptionParser(option_list=option_list))

if (is.null(opt$input)){
  print_help(opt_parser)
  stop("At least one argument must be supplied (input expression file).", call.=FALSE)
}
if (is.null(opt$sample_group)){
  stop("At least one argument must be supplied (input group infomation for DEG analysis).", call.=FALSE)
}

##import files
out_dir<-paste(gsub("/$","",opt$out_dir),"/",sep="")
logexpr<-read.table(opt$input,header=T,stringsAsFactors=F,row.names=1,check.names=F)

#check exp file is log scale
if(max(logexpr[,1])-min(logexpr[,1])>100){
  stop("DEG anlaysis should be conducted based on expression profile on log scale. Please run log2 first", call.=FALSE)
}

##import sample group file and check 
sample_group<-read.table(opt$sample_group,sep="\t",header=T)

if(length(grep("group",colnames(sample_group)))==0){
stop("No group is identified in sample_group file. Make sure the head of sample_group file is like sample, group1, group2.")
}

##### low_expr_filter
if(opt$low_expr_filter==TRUE){
message("Conduct low expression filtering before DEG analysis. \n ")
message(paste("Genes across all samples with expreesion lower than", opt$low_expr_filter_cutoff, "will be filtered out. ",sep=" "))
logexpr<-logexpr[which(apply(logexpr,1,max)>=as.numeric(as.character(opt$low_expr_filter_cutoff))), ]
}

###################################
##########DEGs pairwise DEGs 
#################################
pvalue <- function(x,ncolA,ncolB) {
    obj<-try(t.test(x[1:ncolA],x[(ncolA+1):(ncolA+ncolB)],var.equal = TRUE), silent=TRUE)
   if (is(obj, "try-error")) return(1) else return(obj$p.value)
   }

message("Conducting DEG analysis. \n ")

groupn<-grep("group",colnames(sample_group))
deg<-0
for ( i in groupn){
compgroup<-combn(unique(sample_group[,i]), 2)
for ( j in 1:ncol(compgroup)){


nam<-paste(compgroup[1,j],"vs",compgroup[2,j],sep="")
groupA<-logexpr[,sample_group$sample[sample_group[,i] %in% compgroup[1,j]]]
groupB<-logexpr[,sample_group$sample[sample_group[,i] %in% compgroup[2,j]]]
group2<-cbind(groupA,groupB)
ncolA<-ncol(groupA)
ncolB<-ncol(groupB)
ttest<-apply(group2,1,function(x){pvalue(x,ncolA,ncolB)})
logfc<-rowMeans(groupA)-rowMeans(groupB)

p <- data.frame(
gene=rownames(logexpr),
versus=rep(paste(compgroup[1,j],"vs",compgroup[2,j],sep=" "),nrow(logexpr)),
ttest=ttest,
log10p=(-log10(ttest)),
logfc=logfc,
sigene= 0)

p$sigene[intersect(which(p$ttest<0.05),which(abs(p$logfc)>=1))]<-1
pdiff<-p[p$sigene==1,]

message(paste("Identify ", nrow(pdiff), " DEGs in ", nam, " comparison.",sep=""))

p$sigene<-ifelse(p$sigene==1,"DEG","nonDEG")
deg<-rbind(deg,pdiff)

#modify output
pout<-p[,c("gene","versus","logfc","log10p","sigene")]
pout$sigene<-as.factor(pout$sigene)
rownames(pout)<-c(1:nrow(pout))

if (opt$output_all_genes==FALSE){
#print only DEGs
pout<-pout[pout$sigene=="DEG",]
}
message(paste("Writing results of DEG analysis of ", nam, " comparison.",sep=""))
write.csv(pout,paste(out_dir,opt$project_code,"_",nam,"_degs.csv",sep=""),row.names=F)
saveRDS(pout,paste(out_dir,opt$project_code,"_",nam,"_degs.rds",sep=""))
}

#clear
groupA<-c()
groupB<-c()
group2<-c()
ttest<-c()
logfc<-c()
nam<-c()
}


deg<-deg[!is.na(deg[,1]),]
if(nrow(deg)==0){
message("No DEGs are identified!")
}else{

degtable<-deg[,c("gene","versus","ttest","logfc")]
colnames(degtable)<-gsub("logfc","log2FC",colnames(degtable))
colnames(degtable)<-gsub("ttest","pvalue",colnames(degtable))
degtable$pvalue<-signif(degtable$pvalue,4)
degtable$log2FC<-signif(degtable$log2FC,4)
rownames(degtable)<-c(1:nrow(degtable))

#write output
write.csv(degtable,paste(out_dir,opt$project_code,"_degs_acrossgroups.csv",sep=""),row.names=F)
### statistics number of up and down regulated genes and outpu
degtable$updown<-ifelse(degtable$log2FC>0,"up","down")
up<-degtable[degtable$updown=="up",]
down<-degtable[degtable$updown=="down",]

degstat<-rbind(
 cbind(tapply(up$updown,as.factor(as.character(up$versus)),length),"upregulated"),
 cbind(tapply(down$updown,as.factor(as.character(down$versus)),length),"downregulated")
)

degstat<-cbind(degstat,rownames(degstat))
degstat<-data.frame(degstat)
rownames(degstat)<-c(1:nrow(degstat))
colnames(degstat)<-c("number","type","versus")
degstat$number<-as.numeric(as.character(degstat$number))
########
write.csv(degstat,paste(out_dir,opt$project_code,"_degs_stats.csv",sep=""),row.names=F)
saveRDS(degstat,paste(out_dir,opt$project_code,"_degs_stats.rds",sep=""))
message("RNAseq_4_pwDEG.R finished!")
}