RNAseqDownstream2report/RNAseq_3_cor.R

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

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 for sample group infomation.The input file containing sample name and group infomation. note colname must be like: sample	group1	group2... "),
   	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)
}

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

#check exp file is log scale
if(max(logexpr[,1])-min(logexpr[,1])>100){
  stop("Correlation anlaysis should be conducted based on expression profile on log scale.", call.=FALSE)
}
#####################
##########correlation#######
correlation<-cor(logexpr)
####finished cor

#write output
write.csv(correlation,paste(out_dir,opt$project_code,"_cor.csv",sep=""))
saveRDS(correlation,paste(out_dir,opt$project_code,"_cor.rds",sep=""))
########