RNAseqDownstream2report/RNAseq_1_ballgown.R

#!/usr/bin/env Rscript
# example:
# Rscript RNAseq_1_ballgown.R -o /home/yuying/rnaseqreport_test -i ./ballgown/ -l FALSE -p test 
# Rscript RNAseq_1_ballgown.R -o /home/yuying/rnaseqreport_test -i ./ballgown/ 


suppressPackageStartupMessages(library("optparse"))
suppressPackageStartupMessages(library("ballgown"))

# 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 directory input of expression files. It is output from ballgown software."),
    make_option(c("-f", "--floor_value"),metavar="number",default=0.01,
        help="A number to add to each value before log2 transformation to avoid infinite value.[default: 0.01]"),
    make_option(c("-l", "--log2_norm"), metavar="TRUE", default=TRUE,
        help="Perform log2 transformation on FPKM value. [default: TRUE]"),
	make_option(c("-p", "--project_code"), type="character",default="rnaseq",
        help="Project code, which is used as prefix of output file. [default: rnaseq]")
		)

# 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 file).", call.=FALSE)
}

#generate FPKM expression profile from ballgown outputs
geballgown_expr <- ballgown(dataDir =  opt$input ,samplePattern = ".*",meas = "all")
expr <- gexpr(geballgown_expr)     
 message("finish ballgown\n")

#remove _1P and FPKM from colnames, _1P is from alicloud app, FPKM is added due to default output of stringtie/ballgown. 
nam<-colnames(expr)
nam<-gsub("_1P$","",nam)
nam<-gsub("^FPKM.","",nam)
colnames(expr) <- nam

out_dir<-paste(gsub("/$","",opt$out_dir),"/",sep="")

if(opt$log2_norm==TRUE){
message("start log2 transformation\n")

logexpr<-apply(expr,2,function(x){log2(x+as.numeric(opt$floor_value))})
logexpr_out<-cbind(rownames(logexpr),round(logexpr,3))
colnames(logexpr_out)[1]<-"Gene"

message("output log2 expression file\n")

write.table(logexpr_out,file = paste(out_dir,opt$project_code,"_geneexp_log2fpkm_floor0p01_c",ncol(logexpr),"r",nrow(logexpr),"_",Sys.Date(),".txt",sep=""),sep="\t",row.names=F,quote=F)
}else{ 

#output expression file with fpkm 
expr<-cbind(rownames(expr),round(expr,3))
colnames(expr)[1]<-"Gene"

message("output fpkm expression file\n")

write.table(expr,file = paste(out_dir,opt$project_code,"_geneexp_fpkm_c",ncol(expr),"r",nrow(expr),"_",Sys.Date(),".txt",sep=""),sep="\t",row.names=F,quote=F)
}