regulon.feats=function(fm, genelist, cnv_annot=NULL, filter.val=3, filtertypes=NULL){
fmname=rownames(fm)
if(!is.null(filtertypes)){
fmname=fmname[!substr(fmname, 1,6)%in%filtertypes]
}
feature.gene=suppressWarnings(do.call(rbind, strsplit(fmname, ":"))[,3])
if(any(grepl(":", genelist)))genelist[grepl(":", genelist)]=do.call(rbind, strsplit(genelist[grepl(":", genelist)], ":"))[,3]
if(!is.null(cnv_annot)){
l.genes=strsplit(cnv_annot[,2], ",")
names(l.genes)=cnv_annot[,1]
st=stack(l.genes)
st=st[st[,1]%in%genelist,]
st[,2]=as.character(st[,2])
}
nested.features=function(gene, fmname, feature.gene){
direct=fmname[feature.gene%in%gene]
if(!is.null(cnv_annot))direct=c(direct, st[st[,1]%in%gene,2])
# order this to numeric, then to gexp-cnv-meth-gnab
n=unique(direct)
types=substr(n, 1, 6)
ord=c("N:GEXP","B:GEXP","B:GNAB", "N:CNVR", "B:CNVR","N:METH", "N:MIRN", "N:SAMP", "B:SAMP", "N:CLIN", "B:CLIN", "N:DRUG")
ord=ord[ord%in%types]
r=unlist(lapply(ord, function(o) n[types%in%o]))
return(r)
}
feats=lapply(genelist, nested.features, fmname, feature.gene)
# Filter if not enough data:
m=fm[rownames(fm)%in%unlist(feats),]
s=is.na(m)
m=m[!rowSums(s)==dim(m)[1]&rowSums(!s)>=filter.val, !colSums(s)==dim(m)[2]]
feats=lapply(feats, function(n)n[n%in%rownames(m)])
names(feats)=genelist
feats=feats[!unlist(lapply(feats, length))<1]
return(feats)
}
get.feat.within.region=function(INPUT, FEAT){
temp=tempfile()
temp2=tempfile()
# write out bed format temp file:
write.table(INPUT, temp, quote = F, row.names = F, col.names = F, sep="\t")
# intersect with features
system(paste0("bedtools intersect -a ",FEAT, " -b ",temp, " -wa -wb > ", temp2)) # genes in region
# read in R
d = read.delim(temp2, header=F, stringsAsFactors = F)
d[,4]=gsub(":.*.", "", d[,4])
d2=d[,c(10, 4)]
unlink(temp)
unlink(temp2)
return(d2)
}
fun.correlation.regulon=function(i, l.regulon.gene, fm, filter.val=3, method.cor="spearman"){
feats=l.regulon.gene[[i]]
genename=names(l.regulon.gene)[i]
x=data.matrix(t(as.matrix(fm[match(feats,rownames(fm)),,drop=F])))
class(x) <- "numeric"
if(dim(x)[2]<2)return(NULL) # this one did not have pairs to test
# filter if binary has less than n samples
if(any(grepl("^B:", colnames(x)))){
check=colSums(x[,grepl("^B:", colnames(x)),drop=F], na.rm = T)
x=x[,!colnames(x)%in%names(check[check<filter.val]),drop=F]
}
# filter if numeric has less than n non-zero samples
if(any(grepl("^N:", colnames(x)))){
check=colSums(x[,grepl("^N:", colnames(x)),drop=F]!=0, na.rm = T)
x=x[,!colnames(x)%in%names(check[check<filter.val]),drop=F]
}
if(dim(x)[2]<2)return(NULL) # this one did not have pairs to test
cor_rho=cor(x,use = "pairwise.complete.obs", method = method.cor)
cor_pval=cor.test.p.m(x, method = method.cor, use = "pairwise.complete.obs")
results=flat_cor_mat(cor_rho, cor_pval)
# respect original order:
results=do.call(rbind, lapply(colnames(x), function(n)results[results[,1]%in%n,]))
rownames(results)=NULL
# filter pairwise result:
results=results[!duplicated(t(apply(results[,1:2], 1, sort))),]
# add gene to cnvr:
# results[grepl(":CNVR:", results[,1]),1]=gsub(":CNVR:", paste0(":CNVR:",genename, "@"), results[grepl(":CNVR:", results[,1]),1])
# results[grepl(":CNVR:", results[,2]),2]=gsub(":CNVR:", paste0(":CNVR:",genename, "@"), results[grepl(":CNVR:", results[,2]),2])
return(results)
}
fun.correlation.extrafeatures=function(i, l.regulon.gene, fm, extrafeatures, filter.val=3, method.cor="spearman"){
feats=l.regulon.gene[[i]]
genename=names(l.regulon.gene)[i]
x=data.matrix(t(as.matrix(fm[match(feats,rownames(fm)),,drop=F])))
class(x) <- "numeric"
# filter if binary has less than n samples
if(any(grepl("^B:", colnames(x)))){
check=colSums(x[,grepl("^B:", colnames(x)),drop=F], na.rm = T)
x=x[,!colnames(x)%in%names(check[check<filter.val]),drop=F]
}
# filter if numeric has less than n non-zero samples
if(any(grepl("^N:", colnames(x)))){
check=colSums(x[,grepl("^N:", colnames(x)),drop=F]!=0, na.rm = T)
x=x[,!colnames(x)%in%names(check[check<filter.val]),drop=F]
}
# remove if in both lists
extrafeatures=extrafeatures[!extrafeatures%in%feats]
# Filter if not enough data (less than 3 not NA):
m2=fm[rownames(fm)%in%extrafeatures,]
s=is.na(m2)
m2=m2[!rowSums(s)==dim(m2)[2]&rowSums(!s)>=filter.val, ]
y=data.matrix(t(as.matrix(m2)))
class(y) <- "numeric"
y=cbind(x,y)
# filter if binary has less than n 3 samples
if(any(grepl("^B:", colnames(y)))){
check=colSums(y[,grepl("^B:", colnames(y)),drop=F], na.rm = T)
y=y[,!colnames(y)%in%names(check[check<filter.val])]
}
# filter if numeric has less than n non-zero samples
if(any(grepl("^N:", colnames(y)))){
check=colSums(y[,grepl("^N:", colnames(y)),drop=F]!=0, na.rm = T)
y=y[,!colnames(y)%in%names(check[check<filter.val]),drop=F]
}
if(dim(x)[2]==0|dim(y)[2]==0)return(NULL) # no pairs to test
cor_rho=cor(x,y,use = "pairwise.complete.obs", method = method.cor)
# compute cor pval
cor_pval=t(apply(x, 2, function(v1){
apply(y, 2, function(v2){
if(sum(complete.cases(cbind(v1, v2)))<filter.val)return(NA)
cor.test.p(v1, v2, method = method.cor, use = "pairwise.complete.obs")
})
}))
results=flat_cor_mat(data.matrix(cor_rho), data.matrix(cor_pval))
results=do.call(rbind, lapply(colnames(x), function(n)results[results[,1]%in%n,]))
rownames(results)=NULL
# remove regulonfeatures, they are coming from another test:
results=results[!(results[,1]%in%colnames(x)&results[,2]%in%colnames(x)),]
# filter pairwise result:
results=results[!duplicated(t(apply(results[,1:2], 1, sort))),]
# add gene to cnvr: BUG for wilcox
# results[grepl(":CNVR:", results[,1]),1]=gsub(":CNVR:", paste0(":CNVR:",genename, "@"), results[grepl(":CNVR:", results[,1]),1])
return(results)
}
p.adjust.datatype=function(results, adjust.method="BH", log10=F, prettyNumbers=T, orderdata=T){
# number of features tested:
datatypes=suppressWarnings(do.call(rbind, strsplit(results[,2], ":"))[,2]) # sometimes more than 3 columns
datatypes2=suppressWarnings(do.call(rbind, strsplit(results[,1], ":"))[,2])
datatypes3=apply(cbind(datatypes2,datatypes),1,paste, collapse=":")
datatypes=datatypes3
results=results[order(datatypes),]
datatypes=datatypes[order(datatypes)]
vals=seq(datatypes)
vals=unlist(lapply(unique(datatypes), function(d)p.adjust(results$p[datatypes%in%d], method = adjust.method)))
if(log10){
vals=abs(log10(vals))
}
if(prettyNumbers){
vals=prettyNum(signif(vals,2))
results[,3]=prettyNum(signif(results[,3],2))
results[,4]=prettyNum(signif(results[,4],2))
}else{
vals=as.numeric(vals)
results[,3]=as.numeric(results[,3])
results[,4]=as.numeric(results[,4])
}
results$adj.p=vals
results$signifCode=""
results$signifCode[as.numeric(vals)<0.1]="."
results$signifCode[as.numeric(vals)<0.05]="*"
results$signifCode[as.numeric(vals)<0.01]="**"
results$signifCode[as.numeric(vals)<0.001]="***"
results$signifCode[as.numeric(vals)<0.0001]="****"
results$datapairs=datatypes
if(orderdata)results=results[order(results[,1], datatypes, -as.numeric(results$cor)>0, as.numeric(results$adj.p), decreasing = F),]
return(results)
}
ft.alt <- function(a, b, c, d, alternative="greater") {
as.numeric(fisher.test(matrix(c(a, b, c, d), 2), alternative = alternative)$p.value)
}
fisher.test.pwpw=function(results, find, fm, fisher.alternative, cores){
results.o=results
# remove gene anno if CNV type, to match with FM
results.o[,1]=gsub("CNVR:.*.@", "CNVR:", results.o[,1])
results.o[,2]=gsub("CNVR:.*.@", "CNVR:", results.o[,2])
featsA=as.list(as.data.frame(t(fm[results.o[find,1],])))
featsB=as.list(as.data.frame(t(fm[results.o[find,2],])))
grid=data.frame(do.call(rbind, mclapply(seq(find), function(i)c(table(featsA[[i]], featsB[[i]])), mc.cores=cores)))
colnames(grid)=c("a", "b", "c", "d")
if(fisher.alternative=="two.sided"){
p.ft <- with(grid, HighSpeedStats::ultrafastfet(a, b, c, d))
results$p[find]=p.ft
}else{
p.ft <- with(grid, mapply(ft.alt, a, b, c, d, fisher.alternative))
results$p[find]=p.ft
}
results$test.method[find]=paste0("Fisher´s exact (", fisher.alternative, ")")
return(results)
}
# pairwise.correlation=function(l.regulon.gene, fm, extrafeatures=NULL, filter.val=3, cores=5, adjust.method="BH", method.cor="spearman", prettyNumbers=T, use.fisher=T, fisher.alternative="two.sided"){
#
# #****************************** Test within regulon ******************************
# r=parallel::mclapply(seq(l.regulon.gene),fun.correlation.regulon,l.regulon.gene, fm, filter.val=filter.val,method.cor=method.cor, mc.cores=cores)
# results=do.call(rbind, r)
# results=results[!is.na(results$p),,drop=F]
#
# if(!is.null(results)){
# results$test.method=paste(method.cor, "correlation")
# results$test.group="gene.features"
#
# # if any pairs are binary-binary, add fisher test p-value
# fa=substr(results[,1], 0, 1)=="B"
# fb=substr(results[,2], 0, 1)=="B"
# find=which(fa&fb)
# findl=fa&fb
#
# if(use.fisher&any(findl)){
# results=fisher.test.pwpw(results, find, fm, fisher.alternative, cores)
# results1=p.adjust.datatype(results[findl,], orderdata=T, adjust.method=adjust.method, prettyNumbers=prettyNumbers)
# results2=p.adjust.datatype(results[!findl,], orderdata=T, adjust.method=adjust.method, prettyNumbers=prettyNumbers)
# results=rbind(results1, results2)
# }else{
# results=p.adjust.datatype(results, orderdata=T, adjust.method=adjust.method, prettyNumbers=prettyNumbers)
# }
# }
# #*********************************************************************************
#
# # Test against extrafeatures if they exist:
# if(!is.null(extrafeatures)){
# r=parallel::mclapply(seq(l.regulon.gene),fun.correlation.extrafeatures, l.regulon.gene, fm, extrafeatures=extrafeatures, filter.val=filter.val,method.cor=method.cor, mc.cores=cores)
# results.e=do.call(rbind, r)
# results.e=results.e[!is.na(results.e$p),]
# results.e$test.method=paste(method.cor, "correlation")
# results.e$test.group="extra.features"
#
# # if any pairs are binary-binary, add fisher test p-value
# fa=substr(results.e[,1], 0, 1)=="B"
# fb=substr(results.e[,2], 0, 1)=="B"
# find=which(fa&fb)
# findl=fa&fb
#
# if(use.fisher&any(findl)){
# results.e=fisher.test.pwpw(results.e, find,fm, fisher.alternative, cores)
# results1=p.adjust.datatype(results.e[findl,], orderdata=T, adjust.method=adjust.method, prettyNumbers=prettyNumbers)
# if(sum(!findl)>0){
# results2=p.adjust.datatype(results.e[!findl,], orderdata=T, adjust.method=adjust.method, prettyNumbers=prettyNumbers)
# results.e=rbind(results1, results2)
# }else{
# results.e=results1
# }
# }else{
# results.e=p.adjust.datatype(results.e, orderdata=T, adjust.method=adjust.method, prettyNumbers=prettyNumbers)
# }
# results=rbind(results, results.e)
# }
#
# if(adjust.method=="BH")adj="FDR"
# if(adjust.method=="bonferroni")adj="FWER"
# colnames(results)[7]=adj
# results[grepl("Fisher", results$test.method),3]=""
#
# # order columnts:
# cols=c("featureA", "featureB", "cor", "p", adj, "signifCode", "test.method", "test.group", "datapairs")
# results=results[,cols]
# return(results)
# }
filter.pairwise.res=function(results, filter.table=NULL, FDR=0.1){
if(is.null(filter.table)){
filter.table=data.frame("pair"=unique(results$datapairs), "FDR"=rep(FDR, length(unique(results$datapairs))), stringsAsFactors = F)
type1=gsub(":.*.", "", filter.table[,1])
type2=gsub(".*.:", "", filter.table[,1])
# usually fewer features
filter.table$FDR[type1=="GNAB"|type2=="GNAB"]=0.25 # mutations can be interesting, high P-value kept
filter.table$FDR[type1=="SAMP"&type2=="SAMP"]=0.1 # these can be very different
filter.table$FDR[type1=="CLIN"&type2=="CLIN"]=0.1 # these can be very different
# exclude, not usually meaningful
filter.table$FDR[type1=="CNVR"&type2=="CNVR"]=0 # remove, highly correlated
filter.table$FDR[type1=="METH"&type2=="METH"]=0 # remove, highly correlated
# SAMP special rules:
filter.table$FDR[type1=="SAMP"&type2=="GEXP"]=1e-30
filter.table$FDR[type1=="SAMP"&type2=="MIRN"]=1e-30
filter.table$FDR[type1=="SAMP"&type2=="METH"]=1e-8
print("default table used: ")
print(filter.table)
}
pairs=results$datapairs
p=unique(pairs)
# make sure all pairs are found in filter.table
if(!all(p%in%filter.table[,1]))stop(paste("Not all pairs were found: ", p[!p%in%filter.table[,1],]))
filter.table=filter.table[match(p, filter.table[,1]),]
# p-value per pairs:
vals=do.call(rbind, lapply(seq(p), function(i){
r=results[pairs%in%p[i],]
filt=as.numeric(r[,5])<as.numeric(filter.table[i,2])
r=r[filt,]
}))
return(vals)
}
pairwise.correlation=function(l.regulon.gene, fm, extrafeatures=NULL, filter.val=3, cores=5, adjust.method="BH", method.cor="spearman", prettyNumbers=T, use.fisher=T, use.wilcox=T, fisher.alternative="two.sided", wilcox.alternative="two.sided"){
#****************************** Test within regulon ******************************
r=parallel::mclapply(seq(l.regulon.gene),fun.correlation.regulon,l.regulon.gene, fm, filter.val=filter.val,method.cor=method.cor, mc.cores=cores)
results=do.call(rbind, r)
results=results[!is.na(results$p),,drop=F]
if(!is.null(results)){
results$test.method=paste(method.cor, "correlation")
results$test.group="gene.features"
# if any pairs are binary-binary, add fisher test p-value
fa=substr(results[,1], 0, 1)=="B"
fb=substr(results[,2], 0, 1)=="B"
find=which(fa&fb)
findl=fa&fb
# if any pairs are binary-numeric, add wilcoxon test p-value
fc=substr(results[,1], 0, 1)=="N"
fd=substr(results[,2], 0, 1)=="B"
fe=substr(results[,1], 0, 1)=="B"
ff=substr(results[,2], 0, 1)=="N"
findw=fc&fd|fe&ff
# compute p-value depending on test type:
if(use.fisher&any(findl)){
# update B vs B tests:
results=fisher.test.pwpw(results, find, fm, fisher.alternative, cores)
}
# compute p-value depending on test type:
if(use.wilcox&any(findw)) {
# update B vs N tests:
results=test.feats.wilcox(features = results, rows = findw, data = fm, wilcox.alternative = wilcox.alternative, cores = cores)
}
# go through each test and adjust p-value per data type:
results=do.call(rbind, lapply(unique(results$test.method), function(type){
p.adjust.datatype(results[results$test.method%in%type,], orderdata=T, adjust.method=adjust.method, prettyNumbers=prettyNumbers)
}))
}
#*********************************************************************************
# Test against extrafeatures if they exist:
if(!is.null(extrafeatures)){
r=parallel::mclapply(seq(l.regulon.gene),fun.correlation.extrafeatures, l.regulon.gene, fm, extrafeatures=extrafeatures, filter.val=filter.val,method.cor=method.cor, mc.cores=cores)
results.e=do.call(rbind, r)
results.e=results.e[!is.na(results.e$p),]
results.e$test.method=paste(method.cor, "correlation")
results.e$test.group="extra.features"
# if any pairs are binary-binary, add fisher test p-value
fa=substr(results.e[,1], 0, 1)=="B"
fb=substr(results.e[,2], 0, 1)=="B"
find=which(fa&fb)
findl=fa&fb
# if any pairs are binary-numeric, add wilcoxon test p-value
fc=substr(results.e[,1], 0, 1)=="N"
fd=substr(results.e[,2], 0, 1)=="B"
fe=substr(results.e[,1], 0, 1)=="B"
ff=substr(results.e[,2], 0, 1)=="N"
findw=fc&fd|fe&ff
# compute p-value depending on test type:
if(use.fisher&any(findl)){
# update B vs B tests:
results.e=fisher.test.pwpw(results.e, find, fm, fisher.alternative, cores)
}
# compute p-value depending on test type:
if(use.wilcox&any(findw)) {
# update B vs N tests:
results.e=test.feats.wilcox(features = results.e, rows = findw, data = fm, wilcox.alternative = wilcox.alternative, cores = cores)
}
# go through each test and adjust p-value per data type
results.e=do.call(rbind, lapply(unique(results.e$test.method), function(type){
p.adjust.datatype(results.e[results.e$test.method%in%type,], orderdata=T, adjust.method=adjust.method, prettyNumbers=prettyNumbers)
}))
results=rbind(results, results.e)
}
if(adjust.method=="BH")adj="FDR"
if(adjust.method=="bonferroni")adj="FWER"
colnames(results)[7]=adj
results[grepl("Fisher", results$test.method),3]=""
# order columnts:
cols=c("featureA", "featureB", "cor", "p", adj, "signifCode", "test.method", "test.group", "datapairs")
results=results[,cols]
return(results)
}
filter.pairwise.res=function(results, filter.table=NULL, FDR=0.1){
if(is.null(filter.table)){
filter.table=data.frame("pair"=unique(results$datapairs), "FDR"=rep(FDR, length(unique(results$datapairs))), stringsAsFactors = F)
type1=gsub(":.*.", "", filter.table[,1])
type2=gsub(".*.:", "", filter.table[,1])
# usually fewer features
filter.table$FDR[type1=="GNAB"|type2=="GNAB"]=0.25 # mutations can be interesting, high P-value kept
filter.table$FDR[type1=="SAMP"&type2=="SAMP"]=0.1 # these can be very different
filter.table$FDR[type1=="CLIN"&type2=="CLIN"]=0.1 # these can be very different
# exclude, not usually meaningful
filter.table$FDR[type1=="CNVR"&type2=="CNVR"]=0 # remove, highly correlated
filter.table$FDR[type1=="METH"&type2=="METH"]=0 # remove, highly correlated
# SAMP special rules:
filter.table$FDR[type1=="SAMP"&type2=="GEXP"]=1e-30
filter.table$FDR[type1=="SAMP"&type2=="MIRN"]=1e-30
filter.table$FDR[type1=="SAMP"&type2=="METH"]=1e-8
print("default table used: ")
print(filter.table)
}
pairs=results$datapairs
p=unique(pairs)
# make sure all pairs are found in filter.table
if(!all(p%in%filter.table[,1]))stop(paste("Not all pairs were found: ", p[!p%in%filter.table[,1],]))
filter.table=filter.table[match(p, filter.table[,1]),]
# p-value per pairs:
vals=do.call(rbind, lapply(seq(p), function(i){
r=results[pairs%in%p[i],]
filt=as.numeric(r[,5])<as.numeric(filter.table[i,2])&!is.na(r[,5])
r=r[filt,]
}))
return(vals)
}
test.feats.wilcox <- function(features, rows=NA, data, wilcox.alternative="two.sided", cores=1) {
if(length(rows)>=1) {
pairs <- features[rows,]
} else pairs <- features
p=unlist(parallel::mclapply(1:nrow(pairs), function(i){
binary=substr(pairs[i,2], 0, 1)=="B"
if(binary) {
x <- as.numeric(data[pairs[i,1],])
y <- as.logical(as.numeric(data[pairs[i,2],]))
}else {
x <- as.numeric(data[pairs[i,2],])
y <- as.logical(as.numeric(data[pairs[i,1],]))
}
a <- vector()
b <- vector()
for(j in 1:length(y)) {
if(is.na(y[j])) next
if(y[j]) a <- c(a, x[j])
if(!y[j]) b <- c(b, x[j])
}
if((is.numeric(a) && is.numeric(b))&&(length(a)>0&&length(b>0))&&sum(is.na(a))!=length(a)&&sum(is.na(b))!=length(b)) {
p <- wilcox.test(a,b, alternative = wilcox.alternative)$p.value
} else {
p=NA # not enough data
}
return(p)
}, mc.cores=cores))
pairs$test.method=paste0("Wilcoxon test (", wilcox.alternative, ")")
pairs$p<- as.numeric(as.character(p))
# return original table but with updated p and test.method
features[rows,]=pairs
return(features)
}
Datasets
Models
