smagexp: MetaRNASeq.R annotate

annotate MetaRNASeq.R @ 14:e5a94bc69bd6 draft

planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160

author	sblanck
date	Tue, 16 May 2017 01:53:32 -0400
parents	93451f832736
children	ef7d98f9eb51

rev	line source
14 e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	1 #!/usr/bin/env Rscript
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	2 # setup R error handling to go to stderr
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	3 options( show.error.messages=F, error = function () { cat( geterrmessage(), file=stderr() ); q( "no", 1, F ) } )
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	4
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	5 # we need that to not crash galaxy with an UTF8 error on German LC settings.
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	6 loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	7
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	8 library("optparse")
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	9
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	10 ##### Read options
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	11 option_list=list(
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	12 make_option("--input",type="character",default="NULL",help="list of rdata objects containing eset objects"),
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	13 make_option("--result",type="character",default=NULL,help="text file containing result of the meta-analysis"),
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	14 make_option("--htmloutput",type="character",default=NULL,help="Output html report"),
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	15 make_option("--htmloutputpath",type="character",default="NULL",help="Path of output html report"),
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	16 make_option("--htmltemplate",type="character",default=NULL,help="html template)")
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	17 );
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	18
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	19 opt_parser = OptionParser(option_list=option_list);
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	20 opt = parse_args(opt_parser);
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	21
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	22 if(is.null(opt$input)){
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	23 print_help(opt_parser)
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	24 stop("input required.", call.=FALSE)
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	25 }
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	26
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	27 #loading libraries
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	28
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	29 suppressPackageStartupMessages(require(metaMA))
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	30 suppressPackageStartupMessages(require(affy))
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	31 suppressPackageStartupMessages(require(annaffy))
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	32 suppressPackageStartupMessages(require(VennDiagram))
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	33 suppressPackageStartupMessages(require(GEOquery))
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	34
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	35 listInput <- trimws( unlist( strsplit(trimws(opt$input), ",") ) )
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	36
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	37 listfiles=vector()
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	38 listfilenames=vector()
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	39
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	40 for (i in 1:length(listInput))
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	41 {
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	42 inputFileInfo <- unlist( strsplit( listInput[i], ';' ) )
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	43 listfiles=c(listfiles,inputFileInfo[1])
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	44 listfilenames=c(listfilenames,inputFileInfo[2])
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	45 }
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	46
2 93451f832736 Uploaded sblanck parents: diff changeset	47 cargs <- commandArgs()
93451f832736 Uploaded sblanck parents: diff changeset	48 cargs <- cargs[(which(cargs == "--args")+1):length(cargs)]
93451f832736 Uploaded sblanck parents: diff changeset	49 nbargs=length(cargs)
93451f832736 Uploaded sblanck parents: diff changeset	50 listfiles=vector()
93451f832736 Uploaded sblanck parents: diff changeset	51 listfilenames=vector()
93451f832736 Uploaded sblanck parents: diff changeset	52 for (i in seq(1,nbargs-6,2)) {
93451f832736 Uploaded sblanck parents: diff changeset	53 listfiles=c(listfiles,cargs[[i]])
93451f832736 Uploaded sblanck parents: diff changeset	54 listfilenames=c(listfilenames,cargs[[i+1]])
93451f832736 Uploaded sblanck parents: diff changeset	55 }
93451f832736 Uploaded sblanck parents: diff changeset	56 #mod<-cargs[[length(cargs) - 6]]
14 e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	57 outputfile <- opt$result
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	58 result.html = opt$htmloutput
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	59 html.files.path=opt$htmloutputpath
e5a94bc69bd6 planemo upload for repository https://github.com/sblanck/smagexp/tree/master/smagexp_tools commit da1436252b5db4e90c39c95140558f0f93544160 sblanck parents: 2 diff changeset	60 result.template=opt$htmltemplate
2 93451f832736 Uploaded sblanck parents: diff changeset	61
93451f832736 Uploaded sblanck parents: diff changeset	62 alpha=0.05
93451f832736 Uploaded sblanck parents: diff changeset	63
93451f832736 Uploaded sblanck parents: diff changeset	64 #print(comparison)
93451f832736 Uploaded sblanck parents: diff changeset	65
93451f832736 Uploaded sblanck parents: diff changeset	66 listData=lapply(listfiles,read.table)
93451f832736 Uploaded sblanck parents: diff changeset	67 orderData=lapply(listData, function(x) x[order(x[1]), ])
93451f832736 Uploaded sblanck parents: diff changeset	68 rawpval=lapply(orderData,function(x) x[6])
93451f832736 Uploaded sblanck parents: diff changeset	69 rawpval=lapply(rawpval, function(x) as.numeric(unlist(x)))
93451f832736 Uploaded sblanck parents: diff changeset	70
93451f832736 Uploaded sblanck parents: diff changeset	71 DE=list()
93451f832736 Uploaded sblanck parents: diff changeset	72 DE=lapply(orderData, function(x) ifelse(x[7]<=0.05,1,0))
93451f832736 Uploaded sblanck parents: diff changeset	73
93451f832736 Uploaded sblanck parents: diff changeset	74 FC=list()
93451f832736 Uploaded sblanck parents: diff changeset	75 FC=lapply(orderData, function(x) x[3])
93451f832736 Uploaded sblanck parents: diff changeset	76
93451f832736 Uploaded sblanck parents: diff changeset	77 DE=as.data.frame(DE)
93451f832736 Uploaded sblanck parents: diff changeset	78 colnames(DE)=listfilenames
93451f832736 Uploaded sblanck parents: diff changeset	79 FC=as.data.frame(FC)
93451f832736 Uploaded sblanck parents: diff changeset	80 colnames(FC)=listfilenames
93451f832736 Uploaded sblanck parents: diff changeset	81 # the comparison must only have two values and the conds must
93451f832736 Uploaded sblanck parents: diff changeset	82 # be a vector from those values, at least one of each.
93451f832736 Uploaded sblanck parents: diff changeset	83
93451f832736 Uploaded sblanck parents: diff changeset	84 #if (length(comparison) != 2) {
93451f832736 Uploaded sblanck parents: diff changeset	85 # stop("Comparison type must be a tuple: ", cargs[length(cargs) - 8])
93451f832736 Uploaded sblanck parents: diff changeset	86 #}
93451f832736 Uploaded sblanck parents: diff changeset	87
93451f832736 Uploaded sblanck parents: diff changeset	88 sink("/dev/null")
93451f832736 Uploaded sblanck parents: diff changeset	89 dir.create(html.files.path, recursive=TRUE)
93451f832736 Uploaded sblanck parents: diff changeset	90 #library(DESeq)
93451f832736 Uploaded sblanck parents: diff changeset	91 #library(HTSFilter)
93451f832736 Uploaded sblanck parents: diff changeset	92
93451f832736 Uploaded sblanck parents: diff changeset	93 #DE=list()
93451f832736 Uploaded sblanck parents: diff changeset	94 #FC=list()
93451f832736 Uploaded sblanck parents: diff changeset	95 #i=1
93451f832736 Uploaded sblanck parents: diff changeset	96
93451f832736 Uploaded sblanck parents: diff changeset	97 # Open the html output file
93451f832736 Uploaded sblanck parents: diff changeset	98 #file.conn = file(diag.html, open="w")
93451f832736 Uploaded sblanck parents: diff changeset	99
93451f832736 Uploaded sblanck parents: diff changeset	100 #writeLines( c("<html><body>"), file.conn)
93451f832736 Uploaded sblanck parents: diff changeset	101
93451f832736 Uploaded sblanck parents: diff changeset	102 # Perform deseq analysis on each study
93451f832736 Uploaded sblanck parents: diff changeset	103 #for(i in 1:length(listfiles))
93451f832736 Uploaded sblanck parents: diff changeset	104 #{
93451f832736 Uploaded sblanck parents: diff changeset	105 # f=listfiles[i]
93451f832736 Uploaded sblanck parents: diff changeset	106 # fname=listfilenames[i]
93451f832736 Uploaded sblanck parents: diff changeset	107 # study_name=unlist(strsplit(fname,"[.]"))[1]
93451f832736 Uploaded sblanck parents: diff changeset	108 # print(paste0("study.name ",study_name))
93451f832736 Uploaded sblanck parents: diff changeset	109 # d <- read.table(f, sep=" ", header=TRUE, row.names=1)
93451f832736 Uploaded sblanck parents: diff changeset	110 # conds<-sapply(strsplit(colnames(d),"[.]"),FUN=function(x) x[1])
93451f832736 Uploaded sblanck parents: diff changeset	111 # if (length(unique(conds)) != 2) {
93451f832736 Uploaded sblanck parents: diff changeset	112 # warning(as.data.frame(strsplit(colnames(d),"[.]")))
93451f832736 Uploaded sblanck parents: diff changeset	113 # stop("You can only have two columns types: ", paste(conds,collapse=" "))
93451f832736 Uploaded sblanck parents: diff changeset	114 # }
93451f832736 Uploaded sblanck parents: diff changeset	115 # if (!identical(sort(comparison), sort(unique(conds)))) {
93451f832736 Uploaded sblanck parents: diff changeset	116 # stop("Column types must use the two names from Comparison type, and vice versa. Must have at least one of each in the Column types.\nColumn types: ", cargs[2], "\n", "Comparison type: ", cargs[3])
93451f832736 Uploaded sblanck parents: diff changeset	117 # }
93451f832736 Uploaded sblanck parents: diff changeset	118 # if (length(d) != length(conds)) {
93451f832736 Uploaded sblanck parents: diff changeset	119 # stop("Number of total sample columns in counts file must correspond to the columns types field. E.g. if column types is 'kidney,kidney,liver,liver' then number of sample columns in counts file must be 4 as well.")
93451f832736 Uploaded sblanck parents: diff changeset	120 # }
93451f832736 Uploaded sblanck parents: diff changeset	121 #
93451f832736 Uploaded sblanck parents: diff changeset	122 # cds <- newCountDataSet(d, conds)
93451f832736 Uploaded sblanck parents: diff changeset	123 # cds <- estimateSizeFactors(cds)
93451f832736 Uploaded sblanck parents: diff changeset	124 #
93451f832736 Uploaded sblanck parents: diff changeset	125 # cdsBlind <- estimateDispersions( cds, method="blind" )
93451f832736 Uploaded sblanck parents: diff changeset	126 #
93451f832736 Uploaded sblanck parents: diff changeset	127 # if (length(conds) != 2) {
93451f832736 Uploaded sblanck parents: diff changeset	128 # cds <- estimateDispersions( cds )
93451f832736 Uploaded sblanck parents: diff changeset	129 # norep = FALSE
93451f832736 Uploaded sblanck parents: diff changeset	130 # }
93451f832736 Uploaded sblanck parents: diff changeset	131 #
93451f832736 Uploaded sblanck parents: diff changeset	132 # if (length(conds) == 2) {
93451f832736 Uploaded sblanck parents: diff changeset	133 # cds <- estimateDispersions( cds, method=method, sharingMode=mod, fitType="parametric" )
93451f832736 Uploaded sblanck parents: diff changeset	134 # norep = TRUE
93451f832736 Uploaded sblanck parents: diff changeset	135 # }
93451f832736 Uploaded sblanck parents: diff changeset	136 #
93451f832736 Uploaded sblanck parents: diff changeset	137 # filter<-HTSFilter(cds, plot=FALSE)
93451f832736 Uploaded sblanck parents: diff changeset	138 # cds.filter<-filter$filteredData
93451f832736 Uploaded sblanck parents: diff changeset	139 # on.index<-which(filter$on==1)
93451f832736 Uploaded sblanck parents: diff changeset	140 #
93451f832736 Uploaded sblanck parents: diff changeset	141 # res<-as.data.frame(matrix(NA,nrow=nrow(cds),ncol=ncol(cds)))
93451f832736 Uploaded sblanck parents: diff changeset	142 # nbT <- nbinomTest(cds.filter, comparison[1], comparison[2])
93451f832736 Uploaded sblanck parents: diff changeset	143 # colnames(res)<-colnames(nbT)
93451f832736 Uploaded sblanck parents: diff changeset	144 # res[on.index,]<-nbT
93451f832736 Uploaded sblanck parents: diff changeset	145 # #write.table(res[order(res$padj), ], file=outputfile, quote=FALSE, row.names=FALSE, sep="\t")
93451f832736 Uploaded sblanck parents: diff changeset	146 #
93451f832736 Uploaded sblanck parents: diff changeset	147 #
93451f832736 Uploaded sblanck parents: diff changeset	148 # temp.pval.plot = file.path( html.files.path, paste("PvalHist",i,".png",sep=""))
93451f832736 Uploaded sblanck parents: diff changeset	149 # png( temp.pval.plot, width=500, height=500 )
93451f832736 Uploaded sblanck parents: diff changeset	150 # hist(res$pval, breaks=100, col="skyblue", border="slateblue", main="")
93451f832736 Uploaded sblanck parents: diff changeset	151 # dev.off()
93451f832736 Uploaded sblanck parents: diff changeset	152 #
93451f832736 Uploaded sblanck parents: diff changeset	153 # writeLines( c("<h2>P-value histogram for ",study_name,"</h2>"), file.conn)
93451f832736 Uploaded sblanck parents: diff changeset	154 # writeLines( c("<img src='PvalHist",i,".png'><br/><br/>"), file.conn)
93451f832736 Uploaded sblanck parents: diff changeset	155 #
93451f832736 Uploaded sblanck parents: diff changeset	156 # #on enregistre la p-value
93451f832736 Uploaded sblanck parents: diff changeset	157 # rawpval[[study_name]]<-res$pval
93451f832736 Uploaded sblanck parents: diff changeset	158 # DE[[study_name]]<-ifelse(res$padj<=alpha,1,0)
93451f832736 Uploaded sblanck parents: diff changeset	159 # FC[[study_name]]<-res$log2FoldChange
93451f832736 Uploaded sblanck parents: diff changeset	160 #
93451f832736 Uploaded sblanck parents: diff changeset	161 # i=i+1
93451f832736 Uploaded sblanck parents: diff changeset	162 #}
93451f832736 Uploaded sblanck parents: diff changeset	163
93451f832736 Uploaded sblanck parents: diff changeset	164
93451f832736 Uploaded sblanck parents: diff changeset	165 # combinations
93451f832736 Uploaded sblanck parents: diff changeset	166 library(metaRNASeq)
93451f832736 Uploaded sblanck parents: diff changeset	167 fishcomb<-fishercomb(rawpval, BHth=alpha)
93451f832736 Uploaded sblanck parents: diff changeset	168 warning(length(rawpval))
93451f832736 Uploaded sblanck parents: diff changeset	169 invnormcomb<-invnorm(rawpval, nrep=c(8,8), BHth=alpha)
93451f832736 Uploaded sblanck parents: diff changeset	170 #DE[["fishercomb"]]<-ifelse(fishcomb$adjpval<=alpha,1,0)
93451f832736 Uploaded sblanck parents: diff changeset	171 #DE[["invnormcomb"]]<-ifelse(invnormcomb$adjpval<=alpha,1,0)
93451f832736 Uploaded sblanck parents: diff changeset	172
93451f832736 Uploaded sblanck parents: diff changeset	173 signsFC<-mapply(FC,FUN=function(x) sign(x))
93451f832736 Uploaded sblanck parents: diff changeset	174 sumsigns<-apply(signsFC,1,sum)
93451f832736 Uploaded sblanck parents: diff changeset	175 commonsgnFC<-ifelse(abs(sumsigns)==dim(signsFC)[2],sign(sumsigns),0)
93451f832736 Uploaded sblanck parents: diff changeset	176
93451f832736 Uploaded sblanck parents: diff changeset	177 DEresults <- data.frame(DE=DE,"DE.fishercomb"=ifelse(fishcomb$adjpval<=alpha,1,0),"DE.invnorm"=ifelse(invnormcomb$adjpval<=alpha,1,0))
93451f832736 Uploaded sblanck parents: diff changeset	178
93451f832736 Uploaded sblanck parents: diff changeset	179 unionDE <- unique(c(fishcomb$DEindices,invnormcomb$DEindices))
93451f832736 Uploaded sblanck parents: diff changeset	180 FC.selecDE <- data.frame(DEresults[unionDE,],FC[unionDE,],signFC=commonsgnFC[unionDE])
93451f832736 Uploaded sblanck parents: diff changeset	181 keepDE <- FC.selecDE[which(abs(FC.selecDE$signFC)==1),]
93451f832736 Uploaded sblanck parents: diff changeset	182
93451f832736 Uploaded sblanck parents: diff changeset	183 fishcomb_de <- rownames(keepDE)[which(keepDE[,"DE.fishercomb"]==1)]
93451f832736 Uploaded sblanck parents: diff changeset	184 invnorm_de <- rownames(keepDE)[which(keepDE[,"DE.invnorm"]==1)]
93451f832736 Uploaded sblanck parents: diff changeset	185 indstudy_de = list()
93451f832736 Uploaded sblanck parents: diff changeset	186 for (i in 1:length(listfiles)) {
93451f832736 Uploaded sblanck parents: diff changeset	187 currentIndstudy_de = rownames(keepDE)[which(keepDE[,i]==1)]
93451f832736 Uploaded sblanck parents: diff changeset	188 indstudy_de[[listfilenames[i]]]=currentIndstudy_de
93451f832736 Uploaded sblanck parents: diff changeset	189 }
93451f832736 Uploaded sblanck parents: diff changeset	190
93451f832736 Uploaded sblanck parents: diff changeset	191 IDDIRRfishcomb=IDD.IRR(fishcomb_de,indstudy_de)
93451f832736 Uploaded sblanck parents: diff changeset	192 IDDIRRinvnorm=IDD.IRR(invnorm_de,indstudy_de)
93451f832736 Uploaded sblanck parents: diff changeset	193
93451f832736 Uploaded sblanck parents: diff changeset	194 #conflits<-data.frame(ID=listData[[1]][rownames(DEresults),1],Fishercomb=DEresults[["DE.fishercomb"]],Invnormcomb=DEresults[["DE.invnorm"]],sign=commonsgnFC)
93451f832736 Uploaded sblanck parents: diff changeset	195 conflits<-data.frame(ID=listData[[1]][rownames(DEresults),1],DE=DEresults,FC=FC,signFC=commonsgnFC)
93451f832736 Uploaded sblanck parents: diff changeset	196 #write DE outputfile
93451f832736 Uploaded sblanck parents: diff changeset	197 write.table(conflits, outputfile,sep="\t",,row.names=FALSE)
93451f832736 Uploaded sblanck parents: diff changeset	198 library(VennDiagram)
93451f832736 Uploaded sblanck parents: diff changeset	199 DE_num=apply(DEresults, 2, FUN=function(x) which(x==1))
93451f832736 Uploaded sblanck parents: diff changeset	200 venn.plot<-venn.diagram(x=as.list(DE_num),filename=NULL, col="black", fill=1:length(DE_num)+1,alpha=0.6)
93451f832736 Uploaded sblanck parents: diff changeset	201 temp.venn.plot = file.path( html.files.path, paste("venn.png"))
93451f832736 Uploaded sblanck parents: diff changeset	202 png(temp.venn.plot,width=500,height=500)
93451f832736 Uploaded sblanck parents: diff changeset	203 grid.draw(venn.plot)
93451f832736 Uploaded sblanck parents: diff changeset	204 dev.off()
93451f832736 Uploaded sblanck parents: diff changeset	205
93451f832736 Uploaded sblanck parents: diff changeset	206 library(jsonlite)
93451f832736 Uploaded sblanck parents: diff changeset	207 matrixConflits=as.matrix(conflits)
93451f832736 Uploaded sblanck parents: diff changeset	208 datajson=toJSON(matrixConflits,pretty = TRUE)
93451f832736 Uploaded sblanck parents: diff changeset	209 summaryFishcombjson=toJSON(as.matrix(t(IDDIRRfishcomb)),pretty = TRUE)
93451f832736 Uploaded sblanck parents: diff changeset	210 summaryinvnormjson=toJSON(as.matrix(t(IDDIRRinvnorm)),pretty = TRUE)
93451f832736 Uploaded sblanck parents: diff changeset	211
93451f832736 Uploaded sblanck parents: diff changeset	212
93451f832736 Uploaded sblanck parents: diff changeset	213 #vennsplit=strsplit(result.venn,split="/")[[1]]
93451f832736 Uploaded sblanck parents: diff changeset	214 #venn=paste0("./",vennsplit[length(vennsplit)])
93451f832736 Uploaded sblanck parents: diff changeset	215
93451f832736 Uploaded sblanck parents: diff changeset	216
93451f832736 Uploaded sblanck parents: diff changeset	217 vennFilename="venn.png"
93451f832736 Uploaded sblanck parents: diff changeset	218 vennFile=file.path(html.files.path,vennFilename)
93451f832736 Uploaded sblanck parents: diff changeset	219 htmlfile=readChar(result.template, file.info(result.template)$size)
93451f832736 Uploaded sblanck parents: diff changeset	220 htmlfile=gsub(x=htmlfile,pattern = "###DATAJSON###",replacement = datajson, fixed = TRUE)
93451f832736 Uploaded sblanck parents: diff changeset	221 htmlfile=gsub(x=htmlfile,pattern = "###FISHSUMMARYJSON###",replacement = summaryFishcombjson, fixed = TRUE)
93451f832736 Uploaded sblanck parents: diff changeset	222 htmlfile=gsub(x=htmlfile,pattern = "###INVSUMMARYJSON###",replacement = summaryinvnormjson, fixed = TRUE)
93451f832736 Uploaded sblanck parents: diff changeset	223 htmlfile=gsub(x=htmlfile,pattern = "###VENN###",replacement = vennFilename, fixed = TRUE)
93451f832736 Uploaded sblanck parents: diff changeset	224 write(htmlfile,result.html)
93451f832736 Uploaded sblanck parents: diff changeset	225
93451f832736 Uploaded sblanck parents: diff changeset	226 #library(VennDiagram)
93451f832736 Uploaded sblanck parents: diff changeset	227 #flog.threshold(ERROR)
93451f832736 Uploaded sblanck parents: diff changeset	228 #
93451f832736 Uploaded sblanck parents: diff changeset	229 ##venn.plot<-venn.diagram(x = c(res[c(1:(length(res)-3))],meta=list(res$Meta)),filename = v, col = "black", fill = c(1:(length(res)-2)), margin=0.05, alpha = 0.6,imagetype = "png")
93451f832736 Uploaded sblanck parents: diff changeset	230 #dir.create(result.path, showWarnings = TRUE, recursive = FALSE)
93451f832736 Uploaded sblanck parents: diff changeset	231 #
93451f832736 Uploaded sblanck parents: diff changeset	232 #showVenn<-function(liste,file)
93451f832736 Uploaded sblanck parents: diff changeset	233 #{
93451f832736 Uploaded sblanck parents: diff changeset	234 # venn.plot<-venn.diagram(x = liste,
93451f832736 Uploaded sblanck parents: diff changeset	235 # filename = vennFilename, col = "black",
93451f832736 Uploaded sblanck parents: diff changeset	236 # fill = 1:length(liste)+1,
93451f832736 Uploaded sblanck parents: diff changeset	237 # margin=0.05, alpha = 0.6,imagetype = "png")
93451f832736 Uploaded sblanck parents: diff changeset	238 ## png(file);
93451f832736 Uploaded sblanck parents: diff changeset	239 ## grid.draw(venn.plot);
93451f832736 Uploaded sblanck parents: diff changeset	240 ## dev.off();
93451f832736 Uploaded sblanck parents: diff changeset	241 #
93451f832736 Uploaded sblanck parents: diff changeset	242 #}
93451f832736 Uploaded sblanck parents: diff changeset	243 #
93451f832736 Uploaded sblanck parents: diff changeset	244 #l=list()
93451f832736 Uploaded sblanck parents: diff changeset	245 #for(i in 1:length(esets))
93451f832736 Uploaded sblanck parents: diff changeset	246 #{
93451f832736 Uploaded sblanck parents: diff changeset	247 # l[[paste("study",i,sep="")]]<-res[[i]]
93451f832736 Uploaded sblanck parents: diff changeset	248 #}
93451f832736 Uploaded sblanck parents: diff changeset	249 #l[["Meta"]]=res[[length(res)-1]]
93451f832736 Uploaded sblanck parents: diff changeset	250 #showVenn(l,vennFile)
93451f832736 Uploaded sblanck parents: diff changeset	251 #file.copy(vennFilename,result.path)
93451f832736 Uploaded sblanck parents: diff changeset	252
93451f832736 Uploaded sblanck parents: diff changeset	253
93451f832736 Uploaded sblanck parents: diff changeset	254 #writeLines( c("<h2>Venn Plot</h2>"), file.conn)
93451f832736 Uploaded sblanck parents: diff changeset	255 #writeLines( c("<img src='venn.png'><br/><br/>"), file.conn)
93451f832736 Uploaded sblanck parents: diff changeset	256 #writeLines( c("</body></html>"), file.conn)
93451f832736 Uploaded sblanck parents: diff changeset	257 #close(file.conn)
93451f832736 Uploaded sblanck parents: diff changeset	258 #print("passe6")
93451f832736 Uploaded sblanck parents: diff changeset	259 #sink(NULL)

Mercurial > repos > sblanck > smagexp

annotate MetaRNASeq.R @ 14:e5a94bc69bd6 draft