R/gsea.R
In hxin/topOnto: ontology enrichment analysis
Defines functions
read.gct
read.cls
GSEA.GeneRanking
GSEA.EnrichmentScore.weighted
GSEA.EnrichmentScore.weighted.batch
plot.result
GSEA.HeatMapPlot
read.gmt
write.gmt
read.gct<-function(filename=NULL){
content <- readLines(filename)
content <- content[-1]
content <- content[-1]
col.names <- noquote(unlist(strsplit(content[1], "\t")))
col.names <- col.names[c(-1, -2)]
num.cols <- length(col.names)
content <- content[-1]
content <- content[which(content!="")]
num.lines <- length(content)
row.nam <- vector(length=num.lines, mode="character")
row.des <- vector(length=num.lines, mode="character")
m <- matrix(0, nrow=num.lines, ncol=num.cols)
for (i in 1:num.lines) {
line.list <- noquote(unlist(strsplit(content[i], "\t")))
row.nam[i] <- noquote(line.list[1])
row.des[i] <- noquote(line.list[2])
line.list <- line.list[c(-1, -2)]
for (j in 1:length(line.list)) {
m[i, j] <- as.numeric(line.list[j])
}
}
ds <- data.frame(m)
names(ds) <- col.names
row.names(ds) <- row.nam
return(ds)
}
read.cls<-function(file=NULL){
cls.cont <- readLines(file)
num.lines <- length(cls.cont)
class.list <- unlist(strsplit(cls.cont[[3]], " "))
s <- length(class.list)
t <- table(class.list)
l <- length(t)
phen <- vector(length=l, mode="character")
phen.label <- vector(length=l, mode="numeric")
class.v <- vector(length=s, mode="numeric")
for (i in 1:l) {
phen[i] <- noquote(names(t)[i])
phen.label[i] <- i - 1
}
for (i in 1:s) {
for (j in 1:l) {
if (class.list[i] == phen[j]) {
class.v[i] <- phen.label[j]
}
}
}
return(list(phen = phen, class.v = class.v))
}
GSEA.GeneRanking <- function(A, class.labels, nperm=10, reshuffling.type='sample.labels', permutation.type = 0, sigma.correction = "GeneCluster", fraction=1.0, replace=F, reverse.sign= F,random.seed=123456) {
# This function ranks the genes according to the signal to noise ratio for the actual phenotype and also random permutations and bootstrap
# subsamples of both the observed and random phenotypes. It uses matrix operations to implement the signal to noise calculation
# in stages and achieves fast execution speed. It supports two types of permutations: random (unbalanced) and balanced.
# It also supports subsampling and bootstrap by using masking and multiple-count variables. When "fraction" is set to 1 (default)
# the there is no subsampling or boostrapping and the matrix of observed signal to noise ratios will have the same value for
# all permutations. This is wasteful but allows to support all the multiple options with the same code. Notice that the second
# matrix for the null distribution will still have the values for the random permutations
# (null distribution). This mode (fraction = 1.0) is the defaults, the recommended one and the one used in the examples.
# It is also the one that has be tested more thoroughly. The resampling and boostrapping options are intersting to obtain
# smooth estimates of the observed distribution but its is left for the expert user who may want to perform some sanity
# checks before trusting the code.
#
# Inputs:
# A: Matrix of gene expression values (rows are genes, columns are samples)
# class.labels: Phenotype of class disticntion of interest. A vector of binary labels having first the 1's and then the 0's
# gene.labels: gene labels. Vector of probe ids or accession numbers for the rows of the expression matrix
# nperm: Number of random permutations/bootstraps to perform
# permutation.type: Permutation type: 0 = unbalanced, 1 = balanced. For experts only (default: 0)
# sigma.correction: Correction to the signal to noise ratio (Default = GeneCluster, a choice to support the way it was handled in a previous package)
# fraction: Subsampling fraction. Set to 1.0 (no resampling). For experts only (default: 1.0)
# replace: Resampling mode (replacement or not replacement). For experts only (default: F)
# reverse.sign: Reverse direction of gene list (default = F)
#
# Outputs:
# s2n.matrix: Matrix with random permuted or bootstraps signal to noise ratios (rows are genes, columns are permutations or bootstrap subsamplings
# obs.s2n.matrix: Matrix with observed signal to noise ratios (rows are genes, columns are boostraps subsamplings. If fraction is set to 1.0 then all the columns have the same values
# order.matrix: Matrix with the orderings that will sort the columns of the obs.s2n.matrix in decreasing s2n order
# obs.order.matrix: Matrix with the orderings that will sort the columns of the s2n.matrix in decreasing s2n order
#
# The Broad Institute
# SOFTWARE COPYRIGHT NOTICE AGREEMENT
# This software and its documentation are copyright 2003 by the
# Broad Institute/Massachusetts Institute of Technology.
# All rights are reserved.
#
# This software is supplied without any warranty or guaranteed support
# whatsoever. Neither the Broad Institute nor MIT can be responsible for
# its use, misuse, or functionality.
#browser()
set.seed(seed=random.seed, kind = NULL)
A <- A + 0.00000001
N <- nrow(A)
Ns <- ncol(A)
subset.mask <- matrix(0, nrow=Ns, ncol=nperm)
reshuffled.class.labels1 <- matrix(0, nrow=Ns, ncol=nperm)
reshuffled.class.labels2 <- matrix(0, nrow=Ns, ncol=nperm)
class.labels1 <- matrix(0, nrow=Ns, ncol=nperm)
class.labels2 <- matrix(0, nrow=Ns, ncol=nperm)
order.matrix <- matrix(0, nrow = N, ncol = nperm)
obs.order.matrix <- matrix(0, nrow = N, ncol = nperm)
s2n.matrix <- matrix(0, nrow = N, ncol = nperm)
obs.s2n.matrix <- matrix(0, nrow = N, ncol = nperm)
obs.gene.labels <- vector(length = N, mode="character")
obs.gene.descs <- vector(length = N, mode="character")
obs.gene.symbols <- vector(length = N, mode="character")
M1 <- matrix(0, nrow = N, ncol = nperm)
M2 <- matrix(0, nrow = N, ncol = nperm)
S1 <- matrix(0, nrow = N, ncol = nperm)
S2 <- matrix(0, nrow = N, ncol = nperm)
C <- split(class.labels, class.labels)
class1.size <- length(C[[1]])
class2.size <- length(C[[2]])
class1.index <- seq(1, class1.size, 1)
class2.index <- seq(class1.size + 1, class1.size + class2.size, 1)
for (r in 1:nperm) {
class1.subset <- sample(class1.index, size = ceiling(class1.size*fraction), replace = replace)
class2.subset <- sample(class2.index, size = ceiling(class2.size*fraction), replace = replace)
class1.subset.size <- length(class1.subset)
class2.subset.size <- length(class2.subset)
subset.class1 <- rep(0, class1.size)
for (i in 1:class1.size) {
if (is.element(class1.index[i], class1.subset)) {
subset.class1[i] <- 1
}
}
subset.class2 <- rep(0, class2.size)
for (i in 1:class2.size) {
if (is.element(class2.index[i], class2.subset)) {
subset.class2[i] <- 1
}
}
subset.mask[, r] <- as.numeric(c(subset.class1, subset.class2))
fraction.class1 <- class1.size/Ns
fraction.class2 <- class2.size/Ns
if (permutation.type == 0) { # random (unbalanced) permutation
full.subset <- c(class1.subset, class2.subset)
label1.subset <- sample(full.subset, size = Ns * fraction.class1)
reshuffled.class.labels1[, r] <- rep(0, Ns)
reshuffled.class.labels2[, r] <- rep(0, Ns)
class.labels1[, r] <- rep(0, Ns)
class.labels2[, r] <- rep(0, Ns)
for (i in 1:Ns) {
m1 <- sum(!is.na(match(label1.subset, i)))
m2 <- sum(!is.na(match(full.subset, i)))
reshuffled.class.labels1[i, r] <- m1
reshuffled.class.labels2[i, r] <- m2 - m1
if (i <= class1.size) {
class.labels1[i, r] <- m2
class.labels2[i, r] <- 0
} else {
class.labels1[i, r] <- 0
class.labels2[i, r] <- m2
}
}
} else if (permutation.type == 1) { # proportional (balanced) permutation
class1.label1.subset <- sample(class1.subset, size = ceiling(class1.subset.size*fraction.class1))
class2.label1.subset <- sample(class2.subset, size = floor(class2.subset.size*fraction.class1))
reshuffled.class.labels1[, r] <- rep(0, Ns)
reshuffled.class.labels2[, r] <- rep(0, Ns)
class.labels1[, r] <- rep(0, Ns)
class.labels2[, r] <- rep(0, Ns)
for (i in 1:Ns) {
if (i <= class1.size) {
m1 <- sum(!is.na(match(class1.label1.subset, i)))
m2 <- sum(!is.na(match(class1.subset, i)))
reshuffled.class.labels1[i, r] <- m1
reshuffled.class.labels2[i, r] <- m2 - m1
class.labels1[i, r] <- m2
class.labels2[i, r] <- 0
} else {
m1 <- sum(!is.na(match(class2.label1.subset, i)))
m2 <- sum(!is.na(match(class2.subset, i)))
reshuffled.class.labels1[i, r] <- m1
reshuffled.class.labels2[i, r] <- m2 - m1
class.labels1[i, r] <- 0
class.labels2[i, r] <- m2
}
}
}
}
# compute S2N for the random permutation matrix
##from each, randomly take X samples from the two class lables and get the mean as the null value of that gene's expression
P <- reshuffled.class.labels1 * subset.mask
n1 <- sum(P[,1])
M1 <- A %*% P
M1 <- M1/n1
gc()
A2 <- A*A
S1 <- A2 %*% P
S1 <- S1/n1 - M1*M1
S1 <- sqrt(abs((n1/(n1-1)) * S1))
gc()
P <- reshuffled.class.labels2 * subset.mask
n2 <- sum(P[,1])
M2 <- A %*% P
M2 <- M2/n2
gc()
A2 <- A*A
S2 <- A2 %*% P
S2 <- S2/n2 - M2*M2
S2 <- sqrt(abs((n2/(n2-1)) * S2))
rm(P)
rm(A2)
gc()
if (sigma.correction == "GeneCluster") { # small sigma "fix" as used in GeneCluster
S2 <- ifelse(0.2*abs(M2) < S2, S2, 0.2*abs(M2))
S2 <- ifelse(S2 == 0, 0.2, S2)
S1 <- ifelse(0.2*abs(M1) < S1, S1, 0.2*abs(M1))
S1 <- ifelse(S1 == 0, 0.2, S1)
gc()
}
M1 <- M1 - M2
rm(M2)
gc()
S1 <- S1 + S2
rm(S2)
gc()
s2n.matrix <- M1/S1
if (reverse.sign == T) {
s2n.matrix <- - s2n.matrix
}
gc()
for (r in 1:nperm) {
order.matrix[, r] <- order(s2n.matrix[, r], decreasing=T)
}
# compute S2N for the "observed" permutation matrix
P <- class.labels1 * subset.mask
n1 <- sum(P[,1])
M1 <- A %*% P
M1 <- M1/n1
gc()
A2 <- A*A
S1 <- A2 %*% P
S1 <- S1/n1 - M1*M1
S1 <- sqrt(abs((n1/(n1-1)) * S1))
gc()
P <- class.labels2 * subset.mask
n2 <- sum(P[,1])
M2 <- A %*% P
M2 <- M2/n2
gc()
A2 <- A*A
S2 <- A2 %*% P
S2 <- S2/n2 - M2*M2
S2 <- sqrt(abs((n2/(n2-1)) * S2))
rm(P)
rm(A2)
gc()
if (sigma.correction == "GeneCluster") { # small sigma "fix" as used in GeneCluster
S2 <- ifelse(0.2*abs(M2) < S2, S2, 0.2*abs(M2))
S2 <- ifelse(S2 == 0, 0.2, S2)
S1 <- ifelse(0.2*abs(M1) < S1, S1, 0.2*abs(M1))
S1 <- ifelse(S1 == 0, 0.2, S1)
gc()
}
M1 <- M1 - M2
rm(M2)
gc()
S1 <- S1 + S2
rm(S2)
gc()
obs.s2n.matrix <- M1/S1
gc()
if (reverse.sign == T) {
obs.s2n.matrix <- - obs.s2n.matrix
}
for (r in 1:nperm) {
obs.order.matrix[,r] <- order(obs.s2n.matrix[,r], decreasing=T)
}
# browser()
#rescale correl to (-1,1]
for(k in 1:ncol(s2n.matrix)){
x<-s2n.matrix[,k]
s2n.matrix[,k]<-x/max(abs(x))
}
for(k in 1:ncol(obs.s2n.matrix)){
x<-obs.s2n.matrix[,k]
obs.s2n.matrix[,k]<-x/max(abs(x))
}
if(reshuffling.type=='gene.labels'){
for (r in 1:nperm) {
order.matrix[,r]<- sample(1:nrow(obs.s2n.matrix))
s2n.matrix[,r] <- obs.s2n.matrix[,1]
}
}
o.m=cbind(obs.order.matrix[,1],order.matrix)
s2n.m = cbind(obs.s2n.matrix[,1],s2n.matrix)
s2n.m.sort.abs<-matrix(NA, nrow = nrow(s2n.m), ncol = ncol(s2n.m))
for(i in 1:ncol(s2n.m)){
s2n.m.sort.abs[,i]<-abs(s2n.m[,i][o.m[,i]])
}
return(list(#s2n.matrix = s2n.matrix,
#obs.s2n.matrix = obs.s2n.matrix,
#order.matrix = order.matrix,
#obs.order.matrix = obs.order.matrix,
o.m=o.m,
s2n.m = s2n.m,
s2n.m.sort.abs = s2n.m.sort.abs,
reshuffled.class.labels1 = reshuffled.class.labels1))
}
GSEA.EnrichmentScore.weighted <- function(gene.list, gene.set, correl.vector = NULL,gene.labels,exp.type=1,scores=NULL,lambda=1) {
#
# Computes the weighted GSEA score of gene.set in gene.list.
# The weighted score type is the exponent of the correlation
# weight: 0 (unweighted = Kolmogorov-Smirnov), 1 (weighted), and 2 (over-weighted). When the score type is 1 or 2 it is
# necessary to input the correlation vector with the values in the same order as in the gene list.
#
# Inputs:
# gene.list: The ordered gene list (e.g. integers indicating the original position in the input dataset)
# gene.set: A gene set (e.g. integers indicating the location of those genes in the input dataset)
# weighted.score.type: Type of score: weight: 0 (unweighted = Kolmogorov-Smirnov), 1 (weighted), and 2 (over-weighted)
# correl.vector: A vector with the coorelations (e.g. signal to noise scores) corresponding to the genes in the gene list
#
# Outputs:
# ES: Enrichment score (real number between -1 and +1)
# arg.ES: Location in gene.list where the peak running enrichment occurs (peak of the "mountain")
# RES: Numerical vector containing the running enrichment score for all locations in the gene list
# tag.indicator: Binary vector indicating the location of the gene sets (1's) in the gene list
#
# The Broad Institute
# SOFTWARE COPYRIGHT NOTICE AGREEMENT
# This software and its documentation are copyright 2003 by the
# Broad Institute/Massachusetts Institute of Technology.
# All rights are reserved.
#
# This software is supplied without any warranty or guaranteed support
# whatsoever. Neither the Broad Institute nor MIT can be responsible for
# its use, misuse, or functionality.
#type = 0 equ step 1 use s2n 2 use weight
browser()
gene.labels<-correl.vector
correl.vector<-sort(correl.vector,decreasing = T)
if(exp.type == 0){
weights=1/scores
t.i=c(0)
}else if(exp.type == 1){
weights=1/scores
t.i=c(1)
}else if(exp.type == 2){
weights<-lambda*(1-scores)
t.i<-match(gene.list, gene.set, nomatch=0)
t.i[t.i>0]<-weights[t.i[t.i>0]]
}
names(weights)<-names(gene.labels[gene.set])
tag.indicator <- sign(match(gene.list, gene.set, nomatch=0)) # notice that the sign is 0 (no tag) or 1 (tag)
no.tag.indicator <- 1 - tag.indicator
N <- length(gene.list)
Nh <- length(gene.set)
Nm <- N - Nh
correl.vector <- abs(correl.vector)
sum.correl.tag <- sum(tag.indicator * correl.vector^t.i)
###nomalize on gs.db
norm.tag <- 1.0/sum.correl.tag
norm.no.tag <- 1.0/Nm
RES <- cumsum(tag.indicator * correl.vector^t.i * norm.tag - no.tag.indicator * norm.no.tag)
max.ES <- max(RES)
min.ES <- min(RES)
if (max.ES > - min.ES) {
# ES <- max.ES
ES <- signif(max.ES, digits = 5)
arg.ES <- which.max(RES)
} else {
# ES <- min.ES
ES <- signif(min.ES, digits=5)
arg.ES <- which.min(RES)
}
browser()
##debug
##position_in_rank_list,gene_name,score_norm,corl,corl^score_norm
a<-names(sort(gene.labels[gene.set],decreasing = T))
names(a)<-match(a,names(sort(gene.labels,T)))
#position in rank list / gene name / correl to class label / comfident score / norm step length(corl^p*norm)
dbg<-paste(names(a),a,gene.labels[a],weights[a],correl.vector[a]^weights*norm.tag,sep='/')
#names(dbg)<-names(a)
#plot(c(0,RES,0),type='b')
##debug
gc()
return(list(ES = ES, arg.ES = arg.ES, RES = RES, indicator = tag.indicator,dbg=dbg))
}
GSEA.EnrichmentScore.weighted.batch <- function(gene.list, gene.set, correl.vector = NULL,correl.vector.sorted=NULL,gene.labels,exp.type=1,scores=NULL,lambda=1) {
# browser()
correl.vector.sorted<-correl.vector.sorted
# for(i in 1:ncol(correl.vector)){
# correl.vector.sorted[,i]<-abs(correl.vector[,i][gene.list[,i]])
# }
if(exp.type == 0){
weights=1/scores
t.i=matrix(rep(0,length(gene.list)),ncol = ncol(gene.list))
}else if(exp.type == 1){
weights=1/scores
t.i=matrix(rep(1,length(gene.list)),ncol = ncol(gene.list))
}else if(exp.type == 2){
weights<-lambda*(1-scores)
t.i<-apply(gene.list,2,match,table=gene.set,nomatch=0)
t.i<-apply(t.i,2,function(x){x[x>0]<-weights[x[x>0]];x})
}
tag.indicator<-sign(apply(gene.list,2,match,table=gene.set,nomatch=0))
no.tag.indicator <- 1 - tag.indicator
N <- nrow(gene.list)
Nh <- length(gene.set)
Nm <- N - Nh
step<-tag.indicator * correl.vector.sorted^t.i
sum.correl.tag <- colSums(step)
###nomalize on gs.db
norm.tag <- 1.0/sum.correl.tag
norm.no.tag <- 1.0/Nm
step<-sweep(step, 2, sum.correl.tag, FUN="/")
RES <- vector(mode = "list",length(norm.tag))
for(i in 1:length(norm.tag)){
RES[[i]]<-cumsum(step[,i] - no.tag.indicator[,i] * norm.no.tag)
#names(RES[[i]])<-names(correl.vector[,i][gene.list[,i]])
}
max.ES <- sapply(RES,max)
min.ES <- sapply(RES,min)
ES=vector(mode = "numeric", length = length(max.ES))
arg.ES=vector(mode = "numeric", length = length(max.ES))
for(i in 1:length(max.ES)){
if (max.ES[i] > - min.ES[i]) {
ES[i] <- signif(max.ES[i], digits = 5)
arg.ES[i] <- which.max(RES[[i]])
} else {
ES[i] <- signif(min.ES[i], digits=5)
arg.ES[i] <- which.min(RES[[i]])
}
}
#browser()
##debug
##position_in_rank_list,gene_name,score_norm,corl,corl^score_norm
# dbg<-vector(mode = "list", length = ncol(gene.list))
# gene.labels<-rownames(correl.vector)
# for(i in 1:ncol(correl.vector)){
# #gene.labels<-correl.vector[,i]
# #names(RES[[i]])<-names(gene.labels[gene.list[,i]])
# hit<-tag.indicator[,i]==1
# gname<-gene.labels[gene.list[,i][hit]]
# position<-which(hit)
# #a<-sort(gene.labels[gene.set],decreasing = T)
# #position<-match(names(a),names(gene.labels[gene.list[,i]]))
# dbg[[i]]<-paste(position,gname,correl.vector[,i][gene.list[,i][hit]],weights[gname],step[,i][position],sep='/')
# #dbg[[i]]<-paste(position,names(a),a,weights[names(a)],step[,i][position],sep='/')
# #names(dbg[[i]])=rep(c("position_in_rank_list/gene_name/correl_to_class/comfident_score/norm_step_length(corl^p*norm)"),length(dbg[[i]]))
# }
# #position in rank list / gene name / correl to class label / comfident score / norm step length(corl^p*norm)
gene.labels<-rownames(correl.vector)
hit<-tag.indicator[,1]==1
gname<-gene.labels[gene.list[,1][hit]]
position<-which(hit)
#position in rank list / gene name / correl to class label / comfident score / norm step length(corl^p*norm)
dbg<-paste(position,gname,correl.vector[,1][gene.list[,1][hit]],weights[gname],step[,1][position],sep='/')
return(list(ES = ES, arg.ES = arg.ES, RES = RES[[1]], indicator = tag.indicator[,1],dbg=dbg))
}
plot.result<-function(A=A,O=O,rl=rl,pty=pty,nom.p.val.threshold=0.05,
output.directory='~/Desktop/tmp/2/',doc.string=NULL,topgs=10,adjust.param=0.5){
if(is.null(doc.string)){
doc.string='gsea'
}
# browser()
pdf(file=paste(output.directory, doc.string, ".global.plots", sep="", collapse=""), height = 10, width = 10)
nf <- layout(matrix(c(1,1,2,3), 2, 2, byrow=T), c(1,1), c(1,1), TRUE)
Ng<-length(rl)
Ns<-ncol(A)
obs.s2n=sort(O$s2n.m[,1],T)
N=length(obs.s2n)
location <- 1:N
max.corr <- max(obs.s2n)
min.corr <- min(obs.s2n)
phen=pty$phen
x <- plot(location, obs.s2n, ylab = "Signal to Noise Ratio (S2N)", xlab = "Gene List Location", main = "Gene List Correlation (S2N) Profile", type = "l", lwd = 2, cex = 0.9, col = 1)
for (i in seq(1, N, 20)) {
lines(c(i, i), c(0, obs.s2n[i]), lwd = 3, cex = 0.9, col = colors()[12]) # shading of correlation plot
}
x <- points(location, obs.s2n, type = "l", lwd = 2, cex = 0.9, col = 1)
lines(c(1, N), c(0, 0), lwd = 2, lty = 1, cex = 0.9, col = 1) # zero correlation horizontal line
temp <- order(abs(obs.s2n), decreasing=T)
arg.correl <- temp[N]
lines(c(arg.correl, arg.correl), c(min.corr, 0.7*max.corr), lwd = 2, lty = 3, cex = 0.9, col = 1) # zero correlation vertical line
area.bias <- signif(100*(sum(obs.s2n[1:arg.correl]) + sum(obs.s2n[arg.correl:N]))/sum(abs(obs.s2n[1:N])), digits=3)
area.phen <- ifelse(area.bias >= 0, phen[1], phen[2])
delta.string <- paste("Corr. Area Bias to \"", area.phen, "\" =", abs(area.bias), "%", sep="", collapse="")
zero.crossing.string <- paste("Zero Crossing at location ", arg.correl, " (", signif(100*arg.correl/N, digits=3), " %)")
leg.txt <- c(delta.string, zero.crossing.string)
legend(x=N/10, y=max.corr, bty="n", bg = "white", legend=leg.txt, cex = 0.9)
leg.txt <- paste("\"", phen[1], "\" ", sep="", collapse="")
text(x=1, y=-0.05*max.corr, adj = c(0, 1), labels=leg.txt, cex = 0.9)
leg.txt <- paste("\"", phen[2], "\" ", sep="", collapse="")
text(x=N, y=0.05*max.corr, adj = c(1, 0), labels=leg.txt, cex = 0.9)
###########
ES.norm=t(sapply(rl,function(x){
x$ES.premut.norm
}))
obs.ES.norm<-sapply(rl,function(x){
x$ES.obs.norm
})
obs.ES<-sapply(rl,function(x){x$ES.obs})
nperm<-ncol(O$reshuffled.class.labels1)
if(Ng>0){
phi.densities.pos <- matrix(0, nrow=512, ncol=nperm)
phi.densities.neg <- matrix(0, nrow=512, ncol=nperm)
obs.phi.densities.pos <- matrix(0, nrow=512, ncol=nperm)
obs.phi.densities.neg <- matrix(0, nrow=512, ncol=nperm)
phi.density.mean.pos <- vector(length=512, mode = "numeric")
phi.density.mean.neg <- vector(length=512, mode = "numeric")
obs.phi.density.mean.pos <- vector(length=512, mode = "numeric")
obs.phi.density.mean.neg <- vector(length=512, mode = "numeric")
phi.density.median.pos <- vector(length=512, mode = "numeric")
phi.density.median.neg <- vector(length=512, mode = "numeric")
obs.phi.density.median.pos <- vector(length=512, mode = "numeric")
obs.phi.density.median.neg <- vector(length=512, mode = "numeric")
x.coor.pos <- vector(length=512, mode = "numeric")
x.coor.neg <- vector(length=512, mode = "numeric")
#browser()
for (i in 1:nperm) {
pos.phi <- ES.norm[ES.norm[, i] >= 0, i]
if (length(pos.phi) > 2) {
temp <- density(pos.phi, adjust=adjust.param, n = 512, from=0, to=3.5)
} else {
temp <- list(x = 3.5*(seq(1, 512) - 1)/512, y = rep(0.001, 512))
}
phi.densities.pos[, i] <- temp$y
norm.factor <- sum(phi.densities.pos[, i])
phi.densities.pos[, i] <- phi.densities.pos[, i]/norm.factor
if (i == 1) {
x.coor.pos <- temp$x
}
neg.phi <- ES.norm[ES.norm[, i] < 0, i]
if (length(neg.phi) > 2) {
temp <- density(neg.phi, adjust=adjust.param, n = 512, from=-3.5, to=0)
} else {
temp <- list(x = 3.5*(seq(1, 512) - 1)/512, y = rep(0.001, 512))
}
phi.densities.neg[, i] <- temp$y
norm.factor <- sum(phi.densities.neg[, i])
phi.densities.neg[, i] <- phi.densities.neg[, i]/norm.factor
if (i == 1) {
x.coor.neg <- temp$x
}
pos.phi <- obs.ES.norm[obs.ES.norm >= 0]
if (length(pos.phi) > 2) {
temp <- density(pos.phi, adjust=adjust.param, n = 512, from=0, to=3.5)
} else {
temp <- list(x = 3.5*(seq(1, 512) - 1)/512, y = rep(0.001, 512))
}
obs.phi.densities.pos[, i] <- temp$y
norm.factor <- sum(obs.phi.densities.pos[, i])
obs.phi.densities.pos[, i] <- obs.phi.densities.pos[, i]/norm.factor
neg.phi <- obs.ES.norm[obs.ES.norm < 0]
if (length(neg.phi)> 2) {
temp <- density(neg.phi, adjust=adjust.param, n = 512, from=-3.5, to=0)
} else {
temp <- list(x = 3.5*(seq(1, 512) - 1)/512, y = rep(0.001, 512))
}
obs.phi.densities.neg[, i] <- temp$y
norm.factor <- sum(obs.phi.densities.neg[, i])
obs.phi.densities.neg[, i] <- obs.phi.densities.neg[, i]/norm.factor
}
phi.density.mean.pos <- apply(phi.densities.pos, 1, mean)
phi.density.mean.neg <- apply(phi.densities.neg, 1, mean)
obs.phi.density.mean.pos <- apply(obs.phi.densities.pos, 1, mean)
obs.phi.density.mean.neg <- apply(obs.phi.densities.neg, 1, mean)
phi.density.median.pos <- apply(phi.densities.pos, 1, median)
phi.density.median.neg <- apply(phi.densities.neg, 1, median)
obs.phi.density.median.pos <- apply(obs.phi.densities.pos, 1, median)
obs.phi.density.median.neg <- apply(obs.phi.densities.neg, 1, median)
x <- c(x.coor.neg, x.coor.pos)
x.plot.range <- range(x)
y1 <- c(phi.density.mean.neg, phi.density.mean.pos)
y2 <- c(obs.phi.density.mean.neg, obs.phi.density.mean.pos)
y.plot.range <- c(-0.3*max(c(y1, y2)), max(c(y1, y2)))
#print(c(y.plot.range, max(c(y1, y2)), max(y1), max(y2)))
plot(x, y1, xlim = x.plot.range, ylim = 1.5*y.plot.range, type = "l", lwd = 2, col = 2, xlab = "NES", ylab = "P(NES)", main = "Global Observed and Null Densities (Area Normalized)")
y1.point <- y1[seq(1, length(x), 2)]
y2.point <- y2[seq(2, length(x), 2)]
x1.point <- x[seq(1, length(x), 2)]
x2.point <- x[seq(2, length(x), 2)]
points(x, y1, type = "l", lwd = 2, col = colors()[555])
points(x, y2, type = "l", lwd = 2, col = colors()[29])
for (i in 1:Ng) {
col <- ifelse(obs.ES.norm[i] > 0, 2, 3)
lines(c(obs.ES.norm[i], obs.ES.norm[i]), c(-0.2*max(c(y1, y2)), 0), lwd = 1, lty = 1, col = 1)
}
leg.txt <- paste("Neg. ES: \"", phen[2], " \" ", sep="", collapse="")
text(x=x.plot.range[1], y=-0.25*max(c(y1, y2)), adj = c(0, 1), labels=leg.txt, cex = 0.9)
leg.txt <- paste(" Pos. ES: \"", phen[1], "\" ", sep="", collapse="")
text(x=x.plot.range[2], y=-0.25*max(c(y1, y2)), adj = c(1, 1), labels=leg.txt, cex = 0.9)
leg.txt <- c("Null Density", "Observed Density", "Observed NES values")
c.vec <- c(colors()[555], colors()[29], 1)
lty.vec <- c(1, 1, 1)
lwd.vec <- c(2, 2, 2)
legend(x=0, y=1.5*y.plot.range[2], bty="n", bg = "white", legend=leg.txt, lty = lty.vec, lwd = lwd.vec, col = c.vec, cex = 0.9)
B <- A[order(obs.s2n,decreasing = T),]
if (N > 300) {
C <- rbind(B[1:100,], rep(0, Ns), rep(0, Ns), B[(floor(N/2) - 50 + 1):(floor(N/2) + 50),], rep(0, Ns), rep(0, Ns), B[(N - 100 + 1):N,])
}
rm(B)
GSEA.HeatMapPlot(V = C, col.labels = pty$class.v, col.classes = pty$phen, main = "Heat Map for Genes in Dataset")
dev.off()
}
p.vals=sapply(rl,function(x){x$p})
tag.frac<-sapply(rl,function(x){
if(x$ES.obs>0)
sum(x$GSEA.results$indicator[1:x$GSEA.results$arg.ES[1]])/sum(x$GSEA.results$indicator)
else
sum(x$GSEA.results$indicator[x$GSEA.results$arg.ES[1]:length(x$GSEA.results$indicator)])/sum(x$GSEA.results$indicator)
})
gene.frac<-sapply(rl,function(x){
if(x$ES.obs>0)
x$GSEA.results$arg.ES[1]/length(x$GSEA.results$indicator)
else
1-x$GSEA.results$arg.ES[1]/length(x$GSEA.results$indicator)
})
#browser()
report1 <- data.frame(cbind(names(rl), Term(ONTTERM)[names(rl)],sapply(rl,function(x){sum(x$GSEA.results$indicator)}), round(obs.ES,3),
round(obs.ES.norm,3), round(p.vals,3), round(tag.frac,3), round(gene.frac,3)),row.names = NULL,stringsAsFactors = F)
names(report1) <- c("GS","DEF","SIZE", "ES", "NES", "p", "Tag %", "Gene %")
# print(report)
report2 <- report1
report.index2 <- order(obs.ES.norm, decreasing=T)
for (i in 1:Ng) {
report2[i,] <- report1[report.index2[i],]
}
report3 <- report1
report.index3 <- order(obs.ES.norm, decreasing=F)
for (i in 1:Ng) {
report3[i,] <- report1[report.index3[i],]
}
phen1.rows <- length(obs.ES.norm[obs.ES.norm >= 0])
phen2.rows <- length(obs.ES.norm[obs.ES.norm < 0])
report.phen1 <- report2[1:phen1.rows,]
report.phen2 <- report3[1:phen2.rows,]
if (output.directory != "") {
if (phen1.rows > 0) {
filename <- paste(output.directory, doc.string, ".SUMMARY.RESULTS.REPORT.", pty$phen[1],".txt", sep="", collapse="")
write.table(report.phen1, file = filename, quote=F, row.names=F, sep = "\t")
}
if (phen2.rows > 0) {
filename <- paste(output.directory, doc.string, ".SUMMARY.RESULTS.REPORT.", pty$phen[2],".txt", sep="", collapse="")
write.table(report.phen2, file = filename, quote=F, row.names=F, sep = "\t")
}
}
####plot for each gene set
if (topgs > floor(Ng/2)) {
topgs <- floor(Ng/2)
}
#browser()
result<-list()
result[['report']][[pty$phen[1]]]<-report.phen1
result[['report']][[pty$phen[2]]]<-report.phen2
#browser()
for (i in 1:Ng) {
# result[[names(rl)[i]]][['p']]<-p.vals[i]
# result[[names(rl)[i]]][['ES']]<-obs.ES[i]
# result[[names(rl)[i]]][['ES.norm']]<-obs.ES.norm[i]
# if(names(rl)[i]=='DOID:0014667'){
# browser()
# 1
# }
if (p.vals[i] <= nom.p.val.threshold ||
(is.element(i, c(order(obs.ES.norm,decreasing = T)[1:topgs], order(obs.ES.norm,decreasing = T)[(Ng - topgs + 1): Ng]))))
{
# produce report per gene set
obs.index<-which(rl[[i]]$GSEA.results$indicator==1)
gene.s2n<-obs.s2n[obs.index]
gene.names<-names(gene.s2n)
gene.RES<-rl[[i]]$GSEA.results$RES[obs.index]
loc <- match(i,order(obs.ES.norm,decreasing = T))
if (rl[[i]]$'ES.obs' >= 0) {
core.enrichment<-obs.index <= rl[[i]]$GSEA.results$arg.ES[1]
loc <- match(i,order(obs.ES.norm,decreasing = T))
phen.tag <- phen[1]
} else {
core.enrichment<-obs.index >= rl[[i]]$GSEA.results$arg.ES[1]
loc <- Ng - match(i,order(obs.ES.norm,decreasing = T)) + 1
phen.tag <- phen[2]
}
# browser()
gene.report <- data.frame(cbind(c(1:length(gene.names)), gene.names, obs.index, gene.s2n,
round(1/as.numeric(sapply(rl[[i]]$GSEA.results$dbg,function(x){strsplit(x,split = '/',fixed = T)[[1]][4]}))),
sapply(rl[[i]]$GSEA.results$dbg,function(x){strsplit(x,split = '/',fixed = T)[[1]][4]}),
sapply(rl[[i]]$GSEA.results$dbg,function(x){strsplit(x,split = '/',fixed = T)[[1]][5]}),
gene.RES, core.enrichment),row.names = NULL,stringsAsFactors = F)
names(gene.report) <- c("gene.number", "GENE","LIST LOC", "S2N","CS","expont p","STEP","RES","CORE_ENRICHMENT")
result[[phen.tag]][[names(rl)[i]]][['report']]<-gene.report
result[[phen.tag]][[names(rl)[i]]][['p']]<-p.vals[[i]]
result[[phen.tag]][[names(rl)[i]]][['ES']]<-obs.ES[[i]]
result[[phen.tag]][[names(rl)[i]]][['ES.norm']]<-obs.ES.norm[[i]]
result[[phen.tag]][[names(rl)[i]]][['ES.norm']]<-obs.ES.norm[[i]]
result[[phen.tag]][[names(rl)[i]]][['rl']]<-rl
result[[phen.tag]][[names(rl)[i]]][['obs.s2n']]<- obs.s2n
result[[phen.tag]][[names(rl)[i]]][['obs.ES']]<- obs.ES
result[[phen.tag]][[names(rl)[i]]][['N']]<- N
result[[phen.tag]][[names(rl)[i]]][['obs.index']]<- obs.index
# print(gene.report)
if (output.directory != "") {
filename <- paste(output.directory, doc.string, ".", phen.tag, ".", loc,".", names(rl)[i], ".report.", ".txt", sep="", collapse="")
write.table(gene.report, file = filename, quote=F, row.names=F, sep = "\t")
gs.filename <- paste(output.directory, doc.string, ".", phen.tag, ".", loc,".", names(rl)[i], ".plot", ".pdf", sep="", collapse="")
pdf(file=gs.filename, height = 6, width = 14)
result[[phen.tag]][[names(rl)[i]]][['plot']]<-gs.filename
}
nf <- layout(matrix(c(1,2,3), 1, 3, byrow=T), 1, c(1, 1, 1), TRUE)
ind <- 1:N
min.RES <- min(rl[[i]]$GSEA.results$RES)
max.RES <- max(rl[[i]]$GSEA.results$RES)
if (max.RES < 0.3) max.RES <- 0.3
if (min.RES > -0.3) min.RES <- -0.3
delta <- (max.RES - min.RES)*0.50
min.plot <- min.RES - 2*delta
max.plot <- max.RES
max.corr <- max(obs.s2n)
min.corr <- min(obs.s2n)
Obs.correl.vector.norm <- (obs.s2n - min.corr)/(max.corr - min.corr)*1.25*delta + min.plot
zero.corr.line <- (- min.corr/(max.corr - min.corr))*1.25*delta + min.plot
col <- ifelse(obs.ES[i] > 0, 2, 4)
# Running enrichment plot
sub.string <- paste("Number of genes: ", N, " (in list), ", length(gene.names), " (in gene set)", sep = "", collapse="")
main.string <- paste("Gene Set ", i, ":", names(rl)[i])
plot(ind, rl[[i]]$GSEA.results$RES, main = main.string, sub = sub.string, xlab = "Gene List Index", ylab = "Running Enrichment Score (RES)", xlim=c(1, N), ylim=c(min.plot, max.plot), type = "l", lwd = 2, cex = 1, col = col)
for (j in seq(1, N, 20)) {
lines(c(j, j), c(zero.corr.line, Obs.correl.vector.norm[j]), lwd = 1, cex = 1, col = colors()[12]) # shading of correlation plot
}
lines(c(1, N), c(0, 0), lwd = 1, lty = 2, cex = 1, col = 1) # zero RES line
lines(c(rl[[i]]$GSEA.results$arg.ES[1], rl[[i]]$GSEA.results$arg.ES[1]), c(min.plot, max.plot), lwd = 1, lty = 3, cex = 1, col = col) # max enrichment vertical line
for(j in obs.index){
lines(c(j, j), c(min.plot + 1.25*delta, min.plot + 1.75*delta), lwd = 1, lty = 1, cex = 1, col = 1) # enrichment tags
}
lines(ind, Obs.correl.vector.norm, type = "l", lwd = 1, cex = 1, col = 1)
lines(c(1, N), c(zero.corr.line, zero.corr.line), lwd = 1, lty = 1, cex = 1, col = 1) # zero correlation horizontal line
temp <- order(abs(obs.s2n), decreasing=T)
arg.correl <- temp[N]
lines(c(arg.correl, arg.correl), c(min.plot, max.plot), lwd = 1, lty = 3, cex = 1, col = 3) # zero crossing correlation vertical line
leg.txt <- paste("\"", phen[1], "\" ", sep="", collapse="")
text(x=1, y=min.plot, adj = c(0, 0), labels=leg.txt, cex = 1.0)
leg.txt <- paste("\"", phen[2], "\" ", sep="", collapse="")
text(x=N, y=min.plot, adj = c(1, 0), labels=leg.txt, cex = 1.0)
adjx <- ifelse(obs.ES[i] > 0, 0, 1)
leg.txt <- paste("Peak at ", rl[[i]]$GSEA.results$arg.ES[1], sep="", collapse="")
text(x=rl[[i]]$GSEA.results$arg.ES[1], y=min.plot + 1.8*delta, adj = c(adjx, 0), labels=leg.txt, cex = 1.0)
leg.txt <- paste("Zero crossing at ", arg.correl, sep="", collapse="")
text(x=arg.correl, y=min.plot + 1.95*delta, adj = c(adjx, 0), labels=leg.txt, cex = 1.0)
# nominal p-val histogram
sub.string <- paste("ES =", signif(obs.ES[i], digits = 3), " NES =", signif(obs.ES.norm[i], digits=3), " Nom. p-val=", signif(rl[[i]]$p, digits = 3))
temp <- density(rl[[i]]$ES.premut, adjust=adjust.param)
x.plot.range <- range(temp$x)
y.plot.range <- c(-0.125*max(temp$y), 1.5*max(temp$y))
plot(temp$x, temp$y, type = "l", sub = sub.string, xlim = x.plot.range, ylim = y.plot.range, lwd = 2, col = 2, main = "Gene Set Null Distribution", xlab = "ES", ylab="P(ES)")
x.loc <- which.min(abs(temp$x - obs.ES[i]))
lines(c(obs.ES[i], obs.ES[i]), c(0, temp$y[x.loc]), lwd = 2, lty = 1, cex = 1, col = 1)
lines(x.plot.range, c(0, 0), lwd = 1, lty = 1, cex = 1, col = 1)
leg.txt <- c("Gene Set Null Density", "Observed Gene Set ES value")
c.vec <- c(2, 1)
lty.vec <- c(1, 1)
lwd.vec <- c(2, 2)
legend(x=-0.2, y=y.plot.range[2], bty="n", bg = "white", legend=leg.txt, lty = lty.vec, lwd = lwd.vec, col = c.vec, cex = 1.0)
leg.txt <- paste("Neg. ES \"", phen[2], "\" ", sep="", collapse="")
text(x=x.plot.range[1], y=-0.1*max(temp$y), adj = c(0, 0), labels=leg.txt, cex = 1.0)
leg.txt <- paste(" Pos. ES: \"", phen[1], "\" ", sep="", collapse="")
text(x=x.plot.range[2], y=-0.1*max(temp$y), adj = c(1, 0), labels=leg.txt, cex = 1.0)
# create pinkogram for each gene set
pinko<-A[names(obs.s2n[obs.index]),]
##in case only one gene in the set!
if(is.null(rownames(pinko))){
pinko<-t(as.matrix(pinko))
rownames(pinko)<-names(obs.s2n[obs.index])
}
##
#browser()
GSEA.HeatMapPlot(V = pinko, row.names = rownames(pinko), col.labels = pty$class.v, col.classes = pty$phen, col.names = colnames(pinko), main =" Heat Map for Genes in Gene Set", xlab=" ", ylab=" ")
dev.off()
} # if p.vals thres
} # loop over gene sets
return(result)
}
GSEA.HeatMapPlot <- function(V, row.names = F, col.labels, col.classes, col.names = F, main = " ", xlab=" ", ylab=" ") {
#
# Plots a heatmap "pinkogram" of a gene expression matrix including phenotype vector and gene, sample and phenotype labels
#
# The Broad Institute
# SOFTWARE COPYRIGHT NOTICE AGREEMENT
# This software and its documentation are copyright 2003 by the
# Broad Institute/Massachusetts Institute of Technology.
# All rights are reserved.
#
# This software is supplied without any warranty or guaranteed support
# whatsoever. Neither the Broad Institute nor MIT can be responsible for
# its use, misuse, or functionality.
n.rows <- length(V[,1])
n.cols <- length(V[1,])
row.mean <- apply(V, MARGIN=1, FUN=mean)
row.sd <- apply(V, MARGIN=1, FUN=sd)
row.n <- length(V[,1])
for (i in 1:n.rows) {
if (row.sd[i] == 0) {
V[i,] <- 0
} else {
V[i,] <- (V[i,] - row.mean[i])/(0.5 * row.sd[i])
}
V[i,] <- ifelse(V[i,] < -6, -6, V[i,])
V[i,] <- ifelse(V[i,] > 6, 6, V[i,])
}
mycol <- c("#0000FF", "#0000FF", "#4040FF", "#7070FF", "#8888FF", "#A9A9FF", "#D5D5FF", "#EEE5EE", "#FFAADA", "#FF9DB0", "#FF7080", "#FF5A5A", "#FF4040", "#FF0D1D", "#FF0000") # blue-pinkogram colors. The first and last are the colors to indicate the class vector (phenotype). This is the 1998-vintage, pre-gene cluster, original pinkogram color map
mid.range.V <- mean(range(V)) - 0.1
heatm <- matrix(0, nrow = n.rows + 1, ncol = n.cols)
heatm[1:n.rows,] <- V[seq(n.rows, 1, -1),]
heatm[n.rows + 1,] <- ifelse(col.labels == 0, 7, -7)
image(1:n.cols, 1:(n.rows + 1), t(heatm), col=mycol, axes=FALSE, main=main, xlab= xlab, ylab=ylab)
if (length(row.names) >= 1) {
numC <- nchar(row.names)
size.row.char <- 35/(n.rows + 5)
size.col.char <- 25/(n.cols + 5)
maxl <- floor(n.rows/1.6)
for (i in 1:n.rows) {
row.names[i] <- substr(row.names[i], 1, maxl)
}
row.names <- c(row.names[seq(n.rows, 1, -1)], "Class")
axis(2, at=1:(n.rows + 1), labels=row.names, adj= 0.5, tick=FALSE, las = 1, cex.axis=size.row.char, font.axis=2, line=-1)
}
if (length(col.names) > 1) {
axis(1, at=1:n.cols, labels=col.names, tick=FALSE, las = 3, cex.axis=size.col.char, font.axis=2, line=-1)
}
C <- split(col.labels, col.labels)
class1.size <- length(C[[1]])
class2.size <- length(C[[2]])
axis(3, at=c(floor(class1.size/2),class1.size + floor(class2.size/2)), labels=col.classes, tick=FALSE, las = 1, cex.axis=1.25, font.axis=2, line=-1)
return()
}
read.gmt<-function(file='/home/xin/Downloads/GSEA-P-R/GeneSetDatabases/HDO_o.gmt'){
l<-list()
gs.db<-readLines(file)
for (i in 1:length(gs.db)) {
temp<-strsplit(gs.db[[i]], "\t")[[1]]
l[[temp[1]]]<-temp[-c(1,2)]
}
l
}
write.gmt<-function(term2geneID,file){
def<-Term(ONTTERM)
s=''
for(i in 1:length(term2geneID)){
n<-names(term2geneID)[i]
s<-paste(s,paste(n,def[n],paste(term2geneID[[i]],collapse = '\t'),sep = '\t'),'\n',sep='')
}
sink(file)
cat(s)
sink()
}
hxin/topOnto documentation built on May 17, 2019, 9:15 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions read.gct read.cls GSEA.GeneRanking GSEA.EnrichmentScore.weighted GSEA.EnrichmentScore.weighted.batch plot.result GSEA.HeatMapPlot read.gmt write.gmt
read.gct<-function(filename=NULL){
content <- readLines(filename)
content <- content[-1]
content <- content[-1]
col.names <- noquote(unlist(strsplit(content[1], "\t")))
col.names <- col.names[c(-1, -2)]
num.cols <- length(col.names)
content <- content[-1]
content <- content[which(content!="")]
num.lines <- length(content)
row.nam <- vector(length=num.lines, mode="character")
row.des <- vector(length=num.lines, mode="character")
m <- matrix(0, nrow=num.lines, ncol=num.cols)
for (i in 1:num.lines) {
line.list <- noquote(unlist(strsplit(content[i], "\t")))
row.nam[i] <- noquote(line.list[1])
row.des[i] <- noquote(line.list[2])
line.list <- line.list[c(-1, -2)]
for (j in 1:length(line.list)) {
m[i, j] <- as.numeric(line.list[j])
}
}
ds <- data.frame(m)
names(ds) <- col.names
row.names(ds) <- row.nam
return(ds)
}
read.cls<-function(file=NULL){
cls.cont <- readLines(file)
num.lines <- length(cls.cont)
class.list <- unlist(strsplit(cls.cont[[3]], " "))
s <- length(class.list)
t <- table(class.list)
l <- length(t)
phen <- vector(length=l, mode="character")
phen.label <- vector(length=l, mode="numeric")
class.v <- vector(length=s, mode="numeric")
for (i in 1:l) {
phen[i] <- noquote(names(t)[i])
phen.label[i] <- i - 1
}
for (i in 1:s) {
for (j in 1:l) {
if (class.list[i] == phen[j]) {
class.v[i] <- phen.label[j]
}
}
}
return(list(phen = phen, class.v = class.v))
}
GSEA.GeneRanking <- function(A, class.labels, nperm=10, reshuffling.type='sample.labels', permutation.type = 0, sigma.correction = "GeneCluster", fraction=1.0, replace=F, reverse.sign= F,random.seed=123456) {
# This function ranks the genes according to the signal to noise ratio for the actual phenotype and also random permutations and bootstrap
# subsamples of both the observed and random phenotypes. It uses matrix operations to implement the signal to noise calculation
# in stages and achieves fast execution speed. It supports two types of permutations: random (unbalanced) and balanced.
# It also supports subsampling and bootstrap by using masking and multiple-count variables. When "fraction" is set to 1 (default)
# the there is no subsampling or boostrapping and the matrix of observed signal to noise ratios will have the same value for
# all permutations. This is wasteful but allows to support all the multiple options with the same code. Notice that the second
# matrix for the null distribution will still have the values for the random permutations
# (null distribution). This mode (fraction = 1.0) is the defaults, the recommended one and the one used in the examples.
# It is also the one that has be tested more thoroughly. The resampling and boostrapping options are intersting to obtain
# smooth estimates of the observed distribution but its is left for the expert user who may want to perform some sanity
# checks before trusting the code.
#
# Inputs:
# A: Matrix of gene expression values (rows are genes, columns are samples)
# class.labels: Phenotype of class disticntion of interest. A vector of binary labels having first the 1's and then the 0's
# gene.labels: gene labels. Vector of probe ids or accession numbers for the rows of the expression matrix
# nperm: Number of random permutations/bootstraps to perform
# permutation.type: Permutation type: 0 = unbalanced, 1 = balanced. For experts only (default: 0)
# sigma.correction: Correction to the signal to noise ratio (Default = GeneCluster, a choice to support the way it was handled in a previous package)
# fraction: Subsampling fraction. Set to 1.0 (no resampling). For experts only (default: 1.0)
# replace: Resampling mode (replacement or not replacement). For experts only (default: F)
# reverse.sign: Reverse direction of gene list (default = F)
#
# Outputs:
# s2n.matrix: Matrix with random permuted or bootstraps signal to noise ratios (rows are genes, columns are permutations or bootstrap subsamplings
# obs.s2n.matrix: Matrix with observed signal to noise ratios (rows are genes, columns are boostraps subsamplings. If fraction is set to 1.0 then all the columns have the same values
# order.matrix: Matrix with the orderings that will sort the columns of the obs.s2n.matrix in decreasing s2n order
# obs.order.matrix: Matrix with the orderings that will sort the columns of the s2n.matrix in decreasing s2n order
#
# The Broad Institute
# SOFTWARE COPYRIGHT NOTICE AGREEMENT
# This software and its documentation are copyright 2003 by the
# Broad Institute/Massachusetts Institute of Technology.
# All rights are reserved.
#
# This software is supplied without any warranty or guaranteed support
# whatsoever. Neither the Broad Institute nor MIT can be responsible for
# its use, misuse, or functionality.
#browser()
set.seed(seed=random.seed, kind = NULL)
A <- A + 0.00000001
N <- nrow(A)
Ns <- ncol(A)
subset.mask <- matrix(0, nrow=Ns, ncol=nperm)
reshuffled.class.labels1 <- matrix(0, nrow=Ns, ncol=nperm)
reshuffled.class.labels2 <- matrix(0, nrow=Ns, ncol=nperm)
class.labels1 <- matrix(0, nrow=Ns, ncol=nperm)
class.labels2 <- matrix(0, nrow=Ns, ncol=nperm)
order.matrix <- matrix(0, nrow = N, ncol = nperm)
obs.order.matrix <- matrix(0, nrow = N, ncol = nperm)
s2n.matrix <- matrix(0, nrow = N, ncol = nperm)
obs.s2n.matrix <- matrix(0, nrow = N, ncol = nperm)
obs.gene.labels <- vector(length = N, mode="character")
obs.gene.descs <- vector(length = N, mode="character")
obs.gene.symbols <- vector(length = N, mode="character")
M1 <- matrix(0, nrow = N, ncol = nperm)
M2 <- matrix(0, nrow = N, ncol = nperm)
S1 <- matrix(0, nrow = N, ncol = nperm)
S2 <- matrix(0, nrow = N, ncol = nperm)
C <- split(class.labels, class.labels)
class1.size <- length(C[[1]])
class2.size <- length(C[[2]])
class1.index <- seq(1, class1.size, 1)
class2.index <- seq(class1.size + 1, class1.size + class2.size, 1)
for (r in 1:nperm) {
class1.subset <- sample(class1.index, size = ceiling(class1.size*fraction), replace = replace)
class2.subset <- sample(class2.index, size = ceiling(class2.size*fraction), replace = replace)
class1.subset.size <- length(class1.subset)
class2.subset.size <- length(class2.subset)
subset.class1 <- rep(0, class1.size)
for (i in 1:class1.size) {
if (is.element(class1.index[i], class1.subset)) {
subset.class1[i] <- 1
}
}
subset.class2 <- rep(0, class2.size)
for (i in 1:class2.size) {
if (is.element(class2.index[i], class2.subset)) {
subset.class2[i] <- 1
}
}
subset.mask[, r] <- as.numeric(c(subset.class1, subset.class2))
fraction.class1 <- class1.size/Ns
fraction.class2 <- class2.size/Ns
if (permutation.type == 0) { # random (unbalanced) permutation
full.subset <- c(class1.subset, class2.subset)
label1.subset <- sample(full.subset, size = Ns * fraction.class1)
reshuffled.class.labels1[, r] <- rep(0, Ns)
reshuffled.class.labels2[, r] <- rep(0, Ns)
class.labels1[, r] <- rep(0, Ns)
class.labels2[, r] <- rep(0, Ns)
for (i in 1:Ns) {
m1 <- sum(!is.na(match(label1.subset, i)))
m2 <- sum(!is.na(match(full.subset, i)))
reshuffled.class.labels1[i, r] <- m1
reshuffled.class.labels2[i, r] <- m2 - m1
if (i <= class1.size) {
class.labels1[i, r] <- m2
class.labels2[i, r] <- 0
} else {
class.labels1[i, r] <- 0
class.labels2[i, r] <- m2
}
}
} else if (permutation.type == 1) { # proportional (balanced) permutation
class1.label1.subset <- sample(class1.subset, size = ceiling(class1.subset.size*fraction.class1))
class2.label1.subset <- sample(class2.subset, size = floor(class2.subset.size*fraction.class1))
reshuffled.class.labels1[, r] <- rep(0, Ns)
reshuffled.class.labels2[, r] <- rep(0, Ns)
class.labels1[, r] <- rep(0, Ns)
class.labels2[, r] <- rep(0, Ns)
for (i in 1:Ns) {
if (i <= class1.size) {
m1 <- sum(!is.na(match(class1.label1.subset, i)))
m2 <- sum(!is.na(match(class1.subset, i)))
reshuffled.class.labels1[i, r] <- m1
reshuffled.class.labels2[i, r] <- m2 - m1
class.labels1[i, r] <- m2
class.labels2[i, r] <- 0
} else {
m1 <- sum(!is.na(match(class2.label1.subset, i)))
m2 <- sum(!is.na(match(class2.subset, i)))
reshuffled.class.labels1[i, r] <- m1
reshuffled.class.labels2[i, r] <- m2 - m1
class.labels1[i, r] <- 0
class.labels2[i, r] <- m2
}
}
}
}
# compute S2N for the random permutation matrix
##from each, randomly take X samples from the two class lables and get the mean as the null value of that gene's expression
P <- reshuffled.class.labels1 * subset.mask
n1 <- sum(P[,1])
M1 <- A %*% P
M1 <- M1/n1
gc()
A2 <- A*A
S1 <- A2 %*% P
S1 <- S1/n1 - M1*M1
S1 <- sqrt(abs((n1/(n1-1)) * S1))
gc()
P <- reshuffled.class.labels2 * subset.mask
n2 <- sum(P[,1])
M2 <- A %*% P
M2 <- M2/n2
gc()
A2 <- A*A
S2 <- A2 %*% P
S2 <- S2/n2 - M2*M2
S2 <- sqrt(abs((n2/(n2-1)) * S2))
rm(P)
rm(A2)
gc()
if (sigma.correction == "GeneCluster") { # small sigma "fix" as used in GeneCluster
S2 <- ifelse(0.2*abs(M2) < S2, S2, 0.2*abs(M2))
S2 <- ifelse(S2 == 0, 0.2, S2)
S1 <- ifelse(0.2*abs(M1) < S1, S1, 0.2*abs(M1))
S1 <- ifelse(S1 == 0, 0.2, S1)
gc()
}
M1 <- M1 - M2
rm(M2)
gc()
S1 <- S1 + S2
rm(S2)
gc()
s2n.matrix <- M1/S1
if (reverse.sign == T) {
s2n.matrix <- - s2n.matrix
}
gc()
for (r in 1:nperm) {
order.matrix[, r] <- order(s2n.matrix[, r], decreasing=T)
}
# compute S2N for the "observed" permutation matrix
P <- class.labels1 * subset.mask
n1 <- sum(P[,1])
M1 <- A %*% P
M1 <- M1/n1
gc()
A2 <- A*A
S1 <- A2 %*% P
S1 <- S1/n1 - M1*M1
S1 <- sqrt(abs((n1/(n1-1)) * S1))
gc()
P <- class.labels2 * subset.mask
n2 <- sum(P[,1])
M2 <- A %*% P
M2 <- M2/n2
gc()
A2 <- A*A
S2 <- A2 %*% P
S2 <- S2/n2 - M2*M2
S2 <- sqrt(abs((n2/(n2-1)) * S2))
rm(P)
rm(A2)
gc()
if (sigma.correction == "GeneCluster") { # small sigma "fix" as used in GeneCluster
S2 <- ifelse(0.2*abs(M2) < S2, S2, 0.2*abs(M2))
S2 <- ifelse(S2 == 0, 0.2, S2)
S1 <- ifelse(0.2*abs(M1) < S1, S1, 0.2*abs(M1))
S1 <- ifelse(S1 == 0, 0.2, S1)
gc()
}
M1 <- M1 - M2
rm(M2)
gc()
S1 <- S1 + S2
rm(S2)
gc()
obs.s2n.matrix <- M1/S1
gc()
if (reverse.sign == T) {
obs.s2n.matrix <- - obs.s2n.matrix
}
for (r in 1:nperm) {
obs.order.matrix[,r] <- order(obs.s2n.matrix[,r], decreasing=T)
}
# browser()
#rescale correl to (-1,1]
for(k in 1:ncol(s2n.matrix)){
x<-s2n.matrix[,k]
s2n.matrix[,k]<-x/max(abs(x))
}
for(k in 1:ncol(obs.s2n.matrix)){
x<-obs.s2n.matrix[,k]
obs.s2n.matrix[,k]<-x/max(abs(x))
}
if(reshuffling.type=='gene.labels'){
for (r in 1:nperm) {
order.matrix[,r]<- sample(1:nrow(obs.s2n.matrix))
s2n.matrix[,r] <- obs.s2n.matrix[,1]
}
}
o.m=cbind(obs.order.matrix[,1],order.matrix)
s2n.m = cbind(obs.s2n.matrix[,1],s2n.matrix)
s2n.m.sort.abs<-matrix(NA, nrow = nrow(s2n.m), ncol = ncol(s2n.m))
for(i in 1:ncol(s2n.m)){
s2n.m.sort.abs[,i]<-abs(s2n.m[,i][o.m[,i]])
}
return(list(#s2n.matrix = s2n.matrix,
#obs.s2n.matrix = obs.s2n.matrix,
#order.matrix = order.matrix,
#obs.order.matrix = obs.order.matrix,
o.m=o.m,
s2n.m = s2n.m,
s2n.m.sort.abs = s2n.m.sort.abs,
reshuffled.class.labels1 = reshuffled.class.labels1))
}
GSEA.EnrichmentScore.weighted <- function(gene.list, gene.set, correl.vector = NULL,gene.labels,exp.type=1,scores=NULL,lambda=1) {
#
# Computes the weighted GSEA score of gene.set in gene.list.
# The weighted score type is the exponent of the correlation
# weight: 0 (unweighted = Kolmogorov-Smirnov), 1 (weighted), and 2 (over-weighted). When the score type is 1 or 2 it is
# necessary to input the correlation vector with the values in the same order as in the gene list.
#
# Inputs:
# gene.list: The ordered gene list (e.g. integers indicating the original position in the input dataset)
# gene.set: A gene set (e.g. integers indicating the location of those genes in the input dataset)
# weighted.score.type: Type of score: weight: 0 (unweighted = Kolmogorov-Smirnov), 1 (weighted), and 2 (over-weighted)
# correl.vector: A vector with the coorelations (e.g. signal to noise scores) corresponding to the genes in the gene list
#
# Outputs:
# ES: Enrichment score (real number between -1 and +1)
# arg.ES: Location in gene.list where the peak running enrichment occurs (peak of the "mountain")
# RES: Numerical vector containing the running enrichment score for all locations in the gene list
# tag.indicator: Binary vector indicating the location of the gene sets (1's) in the gene list
#
# The Broad Institute
# SOFTWARE COPYRIGHT NOTICE AGREEMENT
# This software and its documentation are copyright 2003 by the
# Broad Institute/Massachusetts Institute of Technology.
# All rights are reserved.
#
# This software is supplied without any warranty or guaranteed support
# whatsoever. Neither the Broad Institute nor MIT can be responsible for
# its use, misuse, or functionality.
#type = 0 equ step 1 use s2n 2 use weight
browser()
gene.labels<-correl.vector
correl.vector<-sort(correl.vector,decreasing = T)
if(exp.type == 0){
weights=1/scores
t.i=c(0)
}else if(exp.type == 1){
weights=1/scores
t.i=c(1)
}else if(exp.type == 2){
weights<-lambda*(1-scores)
t.i<-match(gene.list, gene.set, nomatch=0)
t.i[t.i>0]<-weights[t.i[t.i>0]]
}
names(weights)<-names(gene.labels[gene.set])
tag.indicator <- sign(match(gene.list, gene.set, nomatch=0)) # notice that the sign is 0 (no tag) or 1 (tag)
no.tag.indicator <- 1 - tag.indicator
N <- length(gene.list)
Nh <- length(gene.set)
Nm <- N - Nh
correl.vector <- abs(correl.vector)
sum.correl.tag <- sum(tag.indicator * correl.vector^t.i)
###nomalize on gs.db
norm.tag <- 1.0/sum.correl.tag
norm.no.tag <- 1.0/Nm
RES <- cumsum(tag.indicator * correl.vector^t.i * norm.tag - no.tag.indicator * norm.no.tag)
max.ES <- max(RES)
min.ES <- min(RES)
if (max.ES > - min.ES) {
# ES <- max.ES
ES <- signif(max.ES, digits = 5)
arg.ES <- which.max(RES)
} else {
# ES <- min.ES
ES <- signif(min.ES, digits=5)
arg.ES <- which.min(RES)
}
browser()
##debug
##position_in_rank_list,gene_name,score_norm,corl,corl^score_norm
a<-names(sort(gene.labels[gene.set],decreasing = T))
names(a)<-match(a,names(sort(gene.labels,T)))
#position in rank list / gene name / correl to class label / comfident score / norm step length(corl^p*norm)
dbg<-paste(names(a),a,gene.labels[a],weights[a],correl.vector[a]^weights*norm.tag,sep='/')
#names(dbg)<-names(a)
#plot(c(0,RES,0),type='b')
##debug
gc()
return(list(ES = ES, arg.ES = arg.ES, RES = RES, indicator = tag.indicator,dbg=dbg))
}
GSEA.EnrichmentScore.weighted.batch <- function(gene.list, gene.set, correl.vector = NULL,correl.vector.sorted=NULL,gene.labels,exp.type=1,scores=NULL,lambda=1) {
# browser()
correl.vector.sorted<-correl.vector.sorted
# for(i in 1:ncol(correl.vector)){
# correl.vector.sorted[,i]<-abs(correl.vector[,i][gene.list[,i]])
# }
if(exp.type == 0){
weights=1/scores
t.i=matrix(rep(0,length(gene.list)),ncol = ncol(gene.list))
}else if(exp.type == 1){
weights=1/scores
t.i=matrix(rep(1,length(gene.list)),ncol = ncol(gene.list))
}else if(exp.type == 2){
weights<-lambda*(1-scores)
t.i<-apply(gene.list,2,match,table=gene.set,nomatch=0)
t.i<-apply(t.i,2,function(x){x[x>0]<-weights[x[x>0]];x})
}
tag.indicator<-sign(apply(gene.list,2,match,table=gene.set,nomatch=0))
no.tag.indicator <- 1 - tag.indicator
N <- nrow(gene.list)
Nh <- length(gene.set)
Nm <- N - Nh
step<-tag.indicator * correl.vector.sorted^t.i
sum.correl.tag <- colSums(step)
###nomalize on gs.db
norm.tag <- 1.0/sum.correl.tag
norm.no.tag <- 1.0/Nm
step<-sweep(step, 2, sum.correl.tag, FUN="/")
RES <- vector(mode = "list",length(norm.tag))
for(i in 1:length(norm.tag)){
RES[[i]]<-cumsum(step[,i] - no.tag.indicator[,i] * norm.no.tag)
#names(RES[[i]])<-names(correl.vector[,i][gene.list[,i]])
}
max.ES <- sapply(RES,max)
min.ES <- sapply(RES,min)
ES=vector(mode = "numeric", length = length(max.ES))
arg.ES=vector(mode = "numeric", length = length(max.ES))
for(i in 1:length(max.ES)){
if (max.ES[i] > - min.ES[i]) {
ES[i] <- signif(max.ES[i], digits = 5)
arg.ES[i] <- which.max(RES[[i]])
} else {
ES[i] <- signif(min.ES[i], digits=5)
arg.ES[i] <- which.min(RES[[i]])
}
}
#browser()
##debug
##position_in_rank_list,gene_name,score_norm,corl,corl^score_norm
# dbg<-vector(mode = "list", length = ncol(gene.list))
# gene.labels<-rownames(correl.vector)
# for(i in 1:ncol(correl.vector)){
# #gene.labels<-correl.vector[,i]
# #names(RES[[i]])<-names(gene.labels[gene.list[,i]])
# hit<-tag.indicator[,i]==1
# gname<-gene.labels[gene.list[,i][hit]]
# position<-which(hit)
# #a<-sort(gene.labels[gene.set],decreasing = T)
# #position<-match(names(a),names(gene.labels[gene.list[,i]]))
# dbg[[i]]<-paste(position,gname,correl.vector[,i][gene.list[,i][hit]],weights[gname],step[,i][position],sep='/')
# #dbg[[i]]<-paste(position,names(a),a,weights[names(a)],step[,i][position],sep='/')
# #names(dbg[[i]])=rep(c("position_in_rank_list/gene_name/correl_to_class/comfident_score/norm_step_length(corl^p*norm)"),length(dbg[[i]]))
# }
# #position in rank list / gene name / correl to class label / comfident score / norm step length(corl^p*norm)
gene.labels<-rownames(correl.vector)
hit<-tag.indicator[,1]==1
gname<-gene.labels[gene.list[,1][hit]]
position<-which(hit)
#position in rank list / gene name / correl to class label / comfident score / norm step length(corl^p*norm)
dbg<-paste(position,gname,correl.vector[,1][gene.list[,1][hit]],weights[gname],step[,1][position],sep='/')
return(list(ES = ES, arg.ES = arg.ES, RES = RES[[1]], indicator = tag.indicator[,1],dbg=dbg))
}
plot.result<-function(A=A,O=O,rl=rl,pty=pty,nom.p.val.threshold=0.05,
output.directory='~/Desktop/tmp/2/',doc.string=NULL,topgs=10,adjust.param=0.5){
if(is.null(doc.string)){
doc.string='gsea'
}
# browser()
pdf(file=paste(output.directory, doc.string, ".global.plots", sep="", collapse=""), height = 10, width = 10)
nf <- layout(matrix(c(1,1,2,3), 2, 2, byrow=T), c(1,1), c(1,1), TRUE)
Ng<-length(rl)
Ns<-ncol(A)
obs.s2n=sort(O$s2n.m[,1],T)
N=length(obs.s2n)
location <- 1:N
max.corr <- max(obs.s2n)
min.corr <- min(obs.s2n)
phen=pty$phen
x <- plot(location, obs.s2n, ylab = "Signal to Noise Ratio (S2N)", xlab = "Gene List Location", main = "Gene List Correlation (S2N) Profile", type = "l", lwd = 2, cex = 0.9, col = 1)
for (i in seq(1, N, 20)) {
lines(c(i, i), c(0, obs.s2n[i]), lwd = 3, cex = 0.9, col = colors()[12]) # shading of correlation plot
}
x <- points(location, obs.s2n, type = "l", lwd = 2, cex = 0.9, col = 1)
lines(c(1, N), c(0, 0), lwd = 2, lty = 1, cex = 0.9, col = 1) # zero correlation horizontal line
temp <- order(abs(obs.s2n), decreasing=T)
arg.correl <- temp[N]
lines(c(arg.correl, arg.correl), c(min.corr, 0.7*max.corr), lwd = 2, lty = 3, cex = 0.9, col = 1) # zero correlation vertical line
area.bias <- signif(100*(sum(obs.s2n[1:arg.correl]) + sum(obs.s2n[arg.correl:N]))/sum(abs(obs.s2n[1:N])), digits=3)
area.phen <- ifelse(area.bias >= 0, phen[1], phen[2])
delta.string <- paste("Corr. Area Bias to \"", area.phen, "\" =", abs(area.bias), "%", sep="", collapse="")
zero.crossing.string <- paste("Zero Crossing at location ", arg.correl, " (", signif(100*arg.correl/N, digits=3), " %)")
leg.txt <- c(delta.string, zero.crossing.string)
legend(x=N/10, y=max.corr, bty="n", bg = "white", legend=leg.txt, cex = 0.9)
leg.txt <- paste("\"", phen[1], "\" ", sep="", collapse="")
text(x=1, y=-0.05*max.corr, adj = c(0, 1), labels=leg.txt, cex = 0.9)
leg.txt <- paste("\"", phen[2], "\" ", sep="", collapse="")
text(x=N, y=0.05*max.corr, adj = c(1, 0), labels=leg.txt, cex = 0.9)
###########
ES.norm=t(sapply(rl,function(x){
x$ES.premut.norm
}))
obs.ES.norm<-sapply(rl,function(x){
x$ES.obs.norm
})
obs.ES<-sapply(rl,function(x){x$ES.obs})
nperm<-ncol(O$reshuffled.class.labels1)
if(Ng>0){
phi.densities.pos <- matrix(0, nrow=512, ncol=nperm)
phi.densities.neg <- matrix(0, nrow=512, ncol=nperm)
obs.phi.densities.pos <- matrix(0, nrow=512, ncol=nperm)
obs.phi.densities.neg <- matrix(0, nrow=512, ncol=nperm)
phi.density.mean.pos <- vector(length=512, mode = "numeric")
phi.density.mean.neg <- vector(length=512, mode = "numeric")
obs.phi.density.mean.pos <- vector(length=512, mode = "numeric")
obs.phi.density.mean.neg <- vector(length=512, mode = "numeric")
phi.density.median.pos <- vector(length=512, mode = "numeric")
phi.density.median.neg <- vector(length=512, mode = "numeric")
obs.phi.density.median.pos <- vector(length=512, mode = "numeric")
obs.phi.density.median.neg <- vector(length=512, mode = "numeric")
x.coor.pos <- vector(length=512, mode = "numeric")
x.coor.neg <- vector(length=512, mode = "numeric")
#browser()
for (i in 1:nperm) {
pos.phi <- ES.norm[ES.norm[, i] >= 0, i]
if (length(pos.phi) > 2) {
temp <- density(pos.phi, adjust=adjust.param, n = 512, from=0, to=3.5)
} else {
temp <- list(x = 3.5*(seq(1, 512) - 1)/512, y = rep(0.001, 512))
}
phi.densities.pos[, i] <- temp$y
norm.factor <- sum(phi.densities.pos[, i])
phi.densities.pos[, i] <- phi.densities.pos[, i]/norm.factor
if (i == 1) {
x.coor.pos <- temp$x
}
neg.phi <- ES.norm[ES.norm[, i] < 0, i]
if (length(neg.phi) > 2) {
temp <- density(neg.phi, adjust=adjust.param, n = 512, from=-3.5, to=0)
} else {
temp <- list(x = 3.5*(seq(1, 512) - 1)/512, y = rep(0.001, 512))
}
phi.densities.neg[, i] <- temp$y
norm.factor <- sum(phi.densities.neg[, i])
phi.densities.neg[, i] <- phi.densities.neg[, i]/norm.factor
if (i == 1) {
x.coor.neg <- temp$x
}
pos.phi <- obs.ES.norm[obs.ES.norm >= 0]
if (length(pos.phi) > 2) {
temp <- density(pos.phi, adjust=adjust.param, n = 512, from=0, to=3.5)
} else {
temp <- list(x = 3.5*(seq(1, 512) - 1)/512, y = rep(0.001, 512))
}
obs.phi.densities.pos[, i] <- temp$y
norm.factor <- sum(obs.phi.densities.pos[, i])
obs.phi.densities.pos[, i] <- obs.phi.densities.pos[, i]/norm.factor
neg.phi <- obs.ES.norm[obs.ES.norm < 0]
if (length(neg.phi)> 2) {
temp <- density(neg.phi, adjust=adjust.param, n = 512, from=-3.5, to=0)
} else {
temp <- list(x = 3.5*(seq(1, 512) - 1)/512, y = rep(0.001, 512))
}
obs.phi.densities.neg[, i] <- temp$y
norm.factor <- sum(obs.phi.densities.neg[, i])
obs.phi.densities.neg[, i] <- obs.phi.densities.neg[, i]/norm.factor
}
phi.density.mean.pos <- apply(phi.densities.pos, 1, mean)
phi.density.mean.neg <- apply(phi.densities.neg, 1, mean)
obs.phi.density.mean.pos <- apply(obs.phi.densities.pos, 1, mean)
obs.phi.density.mean.neg <- apply(obs.phi.densities.neg, 1, mean)
phi.density.median.pos <- apply(phi.densities.pos, 1, median)
phi.density.median.neg <- apply(phi.densities.neg, 1, median)
obs.phi.density.median.pos <- apply(obs.phi.densities.pos, 1, median)
obs.phi.density.median.neg <- apply(obs.phi.densities.neg, 1, median)
x <- c(x.coor.neg, x.coor.pos)
x.plot.range <- range(x)
y1 <- c(phi.density.mean.neg, phi.density.mean.pos)
y2 <- c(obs.phi.density.mean.neg, obs.phi.density.mean.pos)
y.plot.range <- c(-0.3*max(c(y1, y2)), max(c(y1, y2)))
#print(c(y.plot.range, max(c(y1, y2)), max(y1), max(y2)))
plot(x, y1, xlim = x.plot.range, ylim = 1.5*y.plot.range, type = "l", lwd = 2, col = 2, xlab = "NES", ylab = "P(NES)", main = "Global Observed and Null Densities (Area Normalized)")
y1.point <- y1[seq(1, length(x), 2)]
y2.point <- y2[seq(2, length(x), 2)]
x1.point <- x[seq(1, length(x), 2)]
x2.point <- x[seq(2, length(x), 2)]
points(x, y1, type = "l", lwd = 2, col = colors()[555])
points(x, y2, type = "l", lwd = 2, col = colors()[29])
for (i in 1:Ng) {
col <- ifelse(obs.ES.norm[i] > 0, 2, 3)
lines(c(obs.ES.norm[i], obs.ES.norm[i]), c(-0.2*max(c(y1, y2)), 0), lwd = 1, lty = 1, col = 1)
}
leg.txt <- paste("Neg. ES: \"", phen[2], " \" ", sep="", collapse="")
text(x=x.plot.range[1], y=-0.25*max(c(y1, y2)), adj = c(0, 1), labels=leg.txt, cex = 0.9)
leg.txt <- paste(" Pos. ES: \"", phen[1], "\" ", sep="", collapse="")
text(x=x.plot.range[2], y=-0.25*max(c(y1, y2)), adj = c(1, 1), labels=leg.txt, cex = 0.9)
leg.txt <- c("Null Density", "Observed Density", "Observed NES values")
c.vec <- c(colors()[555], colors()[29], 1)
lty.vec <- c(1, 1, 1)
lwd.vec <- c(2, 2, 2)
legend(x=0, y=1.5*y.plot.range[2], bty="n", bg = "white", legend=leg.txt, lty = lty.vec, lwd = lwd.vec, col = c.vec, cex = 0.9)
B <- A[order(obs.s2n,decreasing = T),]
if (N > 300) {
C <- rbind(B[1:100,], rep(0, Ns), rep(0, Ns), B[(floor(N/2) - 50 + 1):(floor(N/2) + 50),], rep(0, Ns), rep(0, Ns), B[(N - 100 + 1):N,])
}
rm(B)
GSEA.HeatMapPlot(V = C, col.labels = pty$class.v, col.classes = pty$phen, main = "Heat Map for Genes in Dataset")
dev.off()
}
p.vals=sapply(rl,function(x){x$p})
tag.frac<-sapply(rl,function(x){
if(x$ES.obs>0)
sum(x$GSEA.results$indicator[1:x$GSEA.results$arg.ES[1]])/sum(x$GSEA.results$indicator)
else
sum(x$GSEA.results$indicator[x$GSEA.results$arg.ES[1]:length(x$GSEA.results$indicator)])/sum(x$GSEA.results$indicator)
})
gene.frac<-sapply(rl,function(x){
if(x$ES.obs>0)
x$GSEA.results$arg.ES[1]/length(x$GSEA.results$indicator)
else
1-x$GSEA.results$arg.ES[1]/length(x$GSEA.results$indicator)
})
#browser()
report1 <- data.frame(cbind(names(rl), Term(ONTTERM)[names(rl)],sapply(rl,function(x){sum(x$GSEA.results$indicator)}), round(obs.ES,3),
round(obs.ES.norm,3), round(p.vals,3), round(tag.frac,3), round(gene.frac,3)),row.names = NULL,stringsAsFactors = F)
names(report1) <- c("GS","DEF","SIZE", "ES", "NES", "p", "Tag %", "Gene %")
# print(report)
report2 <- report1
report.index2 <- order(obs.ES.norm, decreasing=T)
for (i in 1:Ng) {
report2[i,] <- report1[report.index2[i],]
}
report3 <- report1
report.index3 <- order(obs.ES.norm, decreasing=F)
for (i in 1:Ng) {
report3[i,] <- report1[report.index3[i],]
}
phen1.rows <- length(obs.ES.norm[obs.ES.norm >= 0])
phen2.rows <- length(obs.ES.norm[obs.ES.norm < 0])
report.phen1 <- report2[1:phen1.rows,]
report.phen2 <- report3[1:phen2.rows,]
if (output.directory != "") {
if (phen1.rows > 0) {
filename <- paste(output.directory, doc.string, ".SUMMARY.RESULTS.REPORT.", pty$phen[1],".txt", sep="", collapse="")
write.table(report.phen1, file = filename, quote=F, row.names=F, sep = "\t")
}
if (phen2.rows > 0) {
filename <- paste(output.directory, doc.string, ".SUMMARY.RESULTS.REPORT.", pty$phen[2],".txt", sep="", collapse="")
write.table(report.phen2, file = filename, quote=F, row.names=F, sep = "\t")
}
}
####plot for each gene set
if (topgs > floor(Ng/2)) {
topgs <- floor(Ng/2)
}
#browser()
result<-list()
result[['report']][[pty$phen[1]]]<-report.phen1
result[['report']][[pty$phen[2]]]<-report.phen2
#browser()
for (i in 1:Ng) {
# result[[names(rl)[i]]][['p']]<-p.vals[i]
# result[[names(rl)[i]]][['ES']]<-obs.ES[i]
# result[[names(rl)[i]]][['ES.norm']]<-obs.ES.norm[i]
# if(names(rl)[i]=='DOID:0014667'){
# browser()
# 1
# }
if (p.vals[i] <= nom.p.val.threshold ||
(is.element(i, c(order(obs.ES.norm,decreasing = T)[1:topgs], order(obs.ES.norm,decreasing = T)[(Ng - topgs + 1): Ng]))))
{
# produce report per gene set
obs.index<-which(rl[[i]]$GSEA.results$indicator==1)
gene.s2n<-obs.s2n[obs.index]
gene.names<-names(gene.s2n)
gene.RES<-rl[[i]]$GSEA.results$RES[obs.index]
loc <- match(i,order(obs.ES.norm,decreasing = T))
if (rl[[i]]$'ES.obs' >= 0) {
core.enrichment<-obs.index <= rl[[i]]$GSEA.results$arg.ES[1]
loc <- match(i,order(obs.ES.norm,decreasing = T))
phen.tag <- phen[1]
} else {
core.enrichment<-obs.index >= rl[[i]]$GSEA.results$arg.ES[1]
loc <- Ng - match(i,order(obs.ES.norm,decreasing = T)) + 1
phen.tag <- phen[2]
}
# browser()
gene.report <- data.frame(cbind(c(1:length(gene.names)), gene.names, obs.index, gene.s2n,
round(1/as.numeric(sapply(rl[[i]]$GSEA.results$dbg,function(x){strsplit(x,split = '/',fixed = T)[[1]][4]}))),
sapply(rl[[i]]$GSEA.results$dbg,function(x){strsplit(x,split = '/',fixed = T)[[1]][4]}),
sapply(rl[[i]]$GSEA.results$dbg,function(x){strsplit(x,split = '/',fixed = T)[[1]][5]}),
gene.RES, core.enrichment),row.names = NULL,stringsAsFactors = F)
names(gene.report) <- c("gene.number", "GENE","LIST LOC", "S2N","CS","expont p","STEP","RES","CORE_ENRICHMENT")
result[[phen.tag]][[names(rl)[i]]][['report']]<-gene.report
result[[phen.tag]][[names(rl)[i]]][['p']]<-p.vals[[i]]
result[[phen.tag]][[names(rl)[i]]][['ES']]<-obs.ES[[i]]
result[[phen.tag]][[names(rl)[i]]][['ES.norm']]<-obs.ES.norm[[i]]
result[[phen.tag]][[names(rl)[i]]][['ES.norm']]<-obs.ES.norm[[i]]
result[[phen.tag]][[names(rl)[i]]][['rl']]<-rl
result[[phen.tag]][[names(rl)[i]]][['obs.s2n']]<- obs.s2n
result[[phen.tag]][[names(rl)[i]]][['obs.ES']]<- obs.ES
result[[phen.tag]][[names(rl)[i]]][['N']]<- N
result[[phen.tag]][[names(rl)[i]]][['obs.index']]<- obs.index
# print(gene.report)
if (output.directory != "") {
filename <- paste(output.directory, doc.string, ".", phen.tag, ".", loc,".", names(rl)[i], ".report.", ".txt", sep="", collapse="")
write.table(gene.report, file = filename, quote=F, row.names=F, sep = "\t")
gs.filename <- paste(output.directory, doc.string, ".", phen.tag, ".", loc,".", names(rl)[i], ".plot", ".pdf", sep="", collapse="")
pdf(file=gs.filename, height = 6, width = 14)
result[[phen.tag]][[names(rl)[i]]][['plot']]<-gs.filename
}
nf <- layout(matrix(c(1,2,3), 1, 3, byrow=T), 1, c(1, 1, 1), TRUE)
ind <- 1:N
min.RES <- min(rl[[i]]$GSEA.results$RES)
max.RES <- max(rl[[i]]$GSEA.results$RES)
if (max.RES < 0.3) max.RES <- 0.3
if (min.RES > -0.3) min.RES <- -0.3
delta <- (max.RES - min.RES)*0.50
min.plot <- min.RES - 2*delta
max.plot <- max.RES
max.corr <- max(obs.s2n)
min.corr <- min(obs.s2n)
Obs.correl.vector.norm <- (obs.s2n - min.corr)/(max.corr - min.corr)*1.25*delta + min.plot
zero.corr.line <- (- min.corr/(max.corr - min.corr))*1.25*delta + min.plot
col <- ifelse(obs.ES[i] > 0, 2, 4)
# Running enrichment plot
sub.string <- paste("Number of genes: ", N, " (in list), ", length(gene.names), " (in gene set)", sep = "", collapse="")
main.string <- paste("Gene Set ", i, ":", names(rl)[i])
plot(ind, rl[[i]]$GSEA.results$RES, main = main.string, sub = sub.string, xlab = "Gene List Index", ylab = "Running Enrichment Score (RES)", xlim=c(1, N), ylim=c(min.plot, max.plot), type = "l", lwd = 2, cex = 1, col = col)
for (j in seq(1, N, 20)) {
lines(c(j, j), c(zero.corr.line, Obs.correl.vector.norm[j]), lwd = 1, cex = 1, col = colors()[12]) # shading of correlation plot
}
lines(c(1, N), c(0, 0), lwd = 1, lty = 2, cex = 1, col = 1) # zero RES line
lines(c(rl[[i]]$GSEA.results$arg.ES[1], rl[[i]]$GSEA.results$arg.ES[1]), c(min.plot, max.plot), lwd = 1, lty = 3, cex = 1, col = col) # max enrichment vertical line
for(j in obs.index){
lines(c(j, j), c(min.plot + 1.25*delta, min.plot + 1.75*delta), lwd = 1, lty = 1, cex = 1, col = 1) # enrichment tags
}
lines(ind, Obs.correl.vector.norm, type = "l", lwd = 1, cex = 1, col = 1)
lines(c(1, N), c(zero.corr.line, zero.corr.line), lwd = 1, lty = 1, cex = 1, col = 1) # zero correlation horizontal line
temp <- order(abs(obs.s2n), decreasing=T)
arg.correl <- temp[N]
lines(c(arg.correl, arg.correl), c(min.plot, max.plot), lwd = 1, lty = 3, cex = 1, col = 3) # zero crossing correlation vertical line
leg.txt <- paste("\"", phen[1], "\" ", sep="", collapse="")
text(x=1, y=min.plot, adj = c(0, 0), labels=leg.txt, cex = 1.0)
leg.txt <- paste("\"", phen[2], "\" ", sep="", collapse="")
text(x=N, y=min.plot, adj = c(1, 0), labels=leg.txt, cex = 1.0)
adjx <- ifelse(obs.ES[i] > 0, 0, 1)
leg.txt <- paste("Peak at ", rl[[i]]$GSEA.results$arg.ES[1], sep="", collapse="")
text(x=rl[[i]]$GSEA.results$arg.ES[1], y=min.plot + 1.8*delta, adj = c(adjx, 0), labels=leg.txt, cex = 1.0)
leg.txt <- paste("Zero crossing at ", arg.correl, sep="", collapse="")
text(x=arg.correl, y=min.plot + 1.95*delta, adj = c(adjx, 0), labels=leg.txt, cex = 1.0)
# nominal p-val histogram
sub.string <- paste("ES =", signif(obs.ES[i], digits = 3), " NES =", signif(obs.ES.norm[i], digits=3), " Nom. p-val=", signif(rl[[i]]$p, digits = 3))
temp <- density(rl[[i]]$ES.premut, adjust=adjust.param)
x.plot.range <- range(temp$x)
y.plot.range <- c(-0.125*max(temp$y), 1.5*max(temp$y))
plot(temp$x, temp$y, type = "l", sub = sub.string, xlim = x.plot.range, ylim = y.plot.range, lwd = 2, col = 2, main = "Gene Set Null Distribution", xlab = "ES", ylab="P(ES)")
x.loc <- which.min(abs(temp$x - obs.ES[i]))
lines(c(obs.ES[i], obs.ES[i]), c(0, temp$y[x.loc]), lwd = 2, lty = 1, cex = 1, col = 1)
lines(x.plot.range, c(0, 0), lwd = 1, lty = 1, cex = 1, col = 1)
leg.txt <- c("Gene Set Null Density", "Observed Gene Set ES value")
c.vec <- c(2, 1)
lty.vec <- c(1, 1)
lwd.vec <- c(2, 2)
legend(x=-0.2, y=y.plot.range[2], bty="n", bg = "white", legend=leg.txt, lty = lty.vec, lwd = lwd.vec, col = c.vec, cex = 1.0)
leg.txt <- paste("Neg. ES \"", phen[2], "\" ", sep="", collapse="")
text(x=x.plot.range[1], y=-0.1*max(temp$y), adj = c(0, 0), labels=leg.txt, cex = 1.0)
leg.txt <- paste(" Pos. ES: \"", phen[1], "\" ", sep="", collapse="")
text(x=x.plot.range[2], y=-0.1*max(temp$y), adj = c(1, 0), labels=leg.txt, cex = 1.0)
# create pinkogram for each gene set
pinko<-A[names(obs.s2n[obs.index]),]
##in case only one gene in the set!
if(is.null(rownames(pinko))){
pinko<-t(as.matrix(pinko))
rownames(pinko)<-names(obs.s2n[obs.index])
}
##
#browser()
GSEA.HeatMapPlot(V = pinko, row.names = rownames(pinko), col.labels = pty$class.v, col.classes = pty$phen, col.names = colnames(pinko), main =" Heat Map for Genes in Gene Set", xlab=" ", ylab=" ")
dev.off()
} # if p.vals thres
} # loop over gene sets
return(result)
}
GSEA.HeatMapPlot <- function(V, row.names = F, col.labels, col.classes, col.names = F, main = " ", xlab=" ", ylab=" ") {
#
# Plots a heatmap "pinkogram" of a gene expression matrix including phenotype vector and gene, sample and phenotype labels
#
# The Broad Institute
# SOFTWARE COPYRIGHT NOTICE AGREEMENT
# This software and its documentation are copyright 2003 by the
# Broad Institute/Massachusetts Institute of Technology.
# All rights are reserved.
#
# This software is supplied without any warranty or guaranteed support
# whatsoever. Neither the Broad Institute nor MIT can be responsible for
# its use, misuse, or functionality.
n.rows <- length(V[,1])
n.cols <- length(V[1,])
row.mean <- apply(V, MARGIN=1, FUN=mean)
row.sd <- apply(V, MARGIN=1, FUN=sd)
row.n <- length(V[,1])
for (i in 1:n.rows) {
if (row.sd[i] == 0) {
V[i,] <- 0
} else {
V[i,] <- (V[i,] - row.mean[i])/(0.5 * row.sd[i])
}
V[i,] <- ifelse(V[i,] < -6, -6, V[i,])
V[i,] <- ifelse(V[i,] > 6, 6, V[i,])
}
mycol <- c("#0000FF", "#0000FF", "#4040FF", "#7070FF", "#8888FF", "#A9A9FF", "#D5D5FF", "#EEE5EE", "#FFAADA", "#FF9DB0", "#FF7080", "#FF5A5A", "#FF4040", "#FF0D1D", "#FF0000") # blue-pinkogram colors. The first and last are the colors to indicate the class vector (phenotype). This is the 1998-vintage, pre-gene cluster, original pinkogram color map
mid.range.V <- mean(range(V)) - 0.1
heatm <- matrix(0, nrow = n.rows + 1, ncol = n.cols)
heatm[1:n.rows,] <- V[seq(n.rows, 1, -1),]
heatm[n.rows + 1,] <- ifelse(col.labels == 0, 7, -7)
image(1:n.cols, 1:(n.rows + 1), t(heatm), col=mycol, axes=FALSE, main=main, xlab= xlab, ylab=ylab)
if (length(row.names) >= 1) {
numC <- nchar(row.names)
size.row.char <- 35/(n.rows + 5)
size.col.char <- 25/(n.cols + 5)
maxl <- floor(n.rows/1.6)
for (i in 1:n.rows) {
row.names[i] <- substr(row.names[i], 1, maxl)
}
row.names <- c(row.names[seq(n.rows, 1, -1)], "Class")
axis(2, at=1:(n.rows + 1), labels=row.names, adj= 0.5, tick=FALSE, las = 1, cex.axis=size.row.char, font.axis=2, line=-1)
}
if (length(col.names) > 1) {
axis(1, at=1:n.cols, labels=col.names, tick=FALSE, las = 3, cex.axis=size.col.char, font.axis=2, line=-1)
}
C <- split(col.labels, col.labels)
class1.size <- length(C[[1]])
class2.size <- length(C[[2]])
axis(3, at=c(floor(class1.size/2),class1.size + floor(class2.size/2)), labels=col.classes, tick=FALSE, las = 1, cex.axis=1.25, font.axis=2, line=-1)
return()
}
read.gmt<-function(file='/home/xin/Downloads/GSEA-P-R/GeneSetDatabases/HDO_o.gmt'){
l<-list()
gs.db<-readLines(file)
for (i in 1:length(gs.db)) {
temp<-strsplit(gs.db[[i]], "\t")[[1]]
l[[temp[1]]]<-temp[-c(1,2)]
}
l
}
write.gmt<-function(term2geneID,file){
def<-Term(ONTTERM)
s=''
for(i in 1:length(term2geneID)){
n<-names(term2geneID)[i]
s<-paste(s,paste(n,def[n],paste(term2geneID[[i]],collapse = '\t'),sep = '\t'),'\n',sep='')
}
sink(file)
cat(s)
sink()
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.