Nothing
ghdist<-function(n,g=0,h=0){
#
# generate n observations from a g-and-h dist.
#
# this function is based on "Rallfun-v34.txt"
# which can be download from https://dornsife.usc.edu/labs/rwilcox/software/
x<-rnorm(n)
if (g>0){
ghdist<-(exp(g*x)-1)*exp(h*x^2/2)/g
}
if(g==0)ghdist<-x*exp(h*x^2/2)
ghdist
}
.pbos<-function(x,beta=.2,omhatx){
# Compute the percentage bend correlation between x and y.
#
# beta is the bending constant for omega sub N.
#
# this function is based on "Rallfun-v34.txt"
# which can be download from https://dornsife.usc.edu/labs/rwilcox/software/
psi<-(x-median(x))/omhatx
i1<-length(psi[psi<(-1)])
i2<-length(psi[psi>1])
sx<-ifelse(psi<(-1),0,x)
sx<-ifelse(psi>1,0,sx)
.pbos<-(sum(sx)+omhatx*(i2-i1))/(length(x)-i1-i2)
.pbos
}
fisher_transfer_test<-function(x1,z1,x0,z0,biasCorrection=TRUE){
# compute p-value with Fisher Z-transformation test
n1<-length(x1)
n0<-length(x0)
p1=cor(x1,z1,method = 'pearson')
p0=cor(x0,z0,method = 'pearson')
if(biasCorrection==TRUE){
p1=p1-p1/(2*(n1-1))
p0=p0-p0/(2*(n0-1))
}
w1=0.5*log((1+p1)/(1-p1))
w0=0.5*log((1+p0)/(1-p0))
delta=sqrt(1/(n1-3)+1/(n0-3))
T=abs(w1-w0)/delta
return(2-2*pnorm(T))
}
st1<-function(x1,z1,x0,z0){
# Compute p-value for the equal correlation test
# with Pearson correlation based on a logistic regression model
# corresponding to two independent groups
n1<-length(x1)
n0<-length(x0)
p1<-n1/(n1+n0)
p2<-1-p1
x1_mean<-mean(x1)
x0_mean<-mean(x0)
z1_mean<-mean(z1)
z0_mean<-mean(z0)
t1=(n1-1)*var(x1)
t2=(n1-1)*var(z1)
c1<-sqrt(t1*t2)
t3=(n0-1)*var(x0)
t4=(n0-1)*var(z0)
c0<-sqrt(t3*t4)
w1<-n1*(x1-x1_mean)*(z1-z1_mean)/c1
w0<-n0*(x0-x0_mean)*(z0-z0_mean)/c0
u=sum(w1)*p2-sum(w0)*p1
w=c(w1,w0)
u_var=p1*p2*var(w)*(n1+n0-1)
T=u^2/u_var
mysign=sign(u)
signedStat = mysign*T
pval=1-pchisq(T,1)
res = list(stat=T, pval=pval, signedStat = signedStat)
return(res)
}
st2<-function(x1,z1,x0,z0){
# Compute p-value for the equal correlation test
# with mad-replacing-Pearson correlation
# based on a logistic regression model
# corresponding to two independent groups
n1<-length(x1)
n0<-length(x0)
p1<-n1/(n1+n0)
p2<-1-p1
x1_median<-median(x1)
x0_median<-median(x0)
z1_median<-median(z1)
z0_median<-median(z0)
c1<-mad(x1)*mad(z1)
c0<-mad(x0)*mad(z0)
w1<-(x1-x1_median)*(z1-z1_median)/c1
w0<-(x0-x0_median)*(z0-z0_median)/c0
u=sum(w1)*p2-sum(w0)*p1
w=c(w1,w0)
u_var=p1*p2*var(w)*(n1+n0-1)
T=u^2/u_var
mysign=sign(u)
signedStat = mysign*T
pval=1-pchisq(T,1)
res = list(stat=T, pval=pval, signedStat = signedStat)
return(res)
}
st3<-function(x1,z1,x0,z0){
# Compute p-value for the equal correlation test
# with percentage bend correlation
# based on a logistic regression model
# corresponding to two independent groups
n1<-length(x1)
n0<-length(x0)
p1<-n1/(n1+n0)
p2<-1-p1
b=0.2
x1_median<-median(x1)
x0_median<-median(x0)
z1_median<-median(z1)
z0_median<-median(z0)
temp1<-sort(abs(x1- x1_median))
omhatx<-temp1[floor((1-b)*n1)]
if(omhatx==0){
omhatx=temp1[n1]
}
temp2<-sort(abs(z1- z1_median))
omhatz<-temp2[floor((1-b)*n1)]
if(omhatz==0){
omhatz=temp2[n1]
}
a1<-(x1-.pbos(x1,b,omhatx))/omhatx
b1<-(z1-.pbos(z1,b,omhatz))/omhatz
a1<-ifelse(a1<=-1,-1,a1)
a1<-ifelse(a1>=1,1,a1)
b1<-ifelse(b1<=-1,-1,b1)
b1<-ifelse(b1>=1,1,b1)
temp3<-sort(abs(x0- x0_median))
omhatx<-temp3[floor((1-b)*n0)]
if(omhatx==0){
omhatx=temp3[n0]
}
temp4<-sort(abs(z0- z0_median))
omhatz<-temp4[floor((1-b)*n0)]
if(omhatz==0){
omhatz=temp4[n0]
}
a0<-(x0-.pbos(x0,b,omhatx))/omhatx
b0<-(z0-.pbos(z0,b,omhatz))/omhatz
a0<-ifelse(a0<=-1,-1,a0)
a0<-ifelse(a0>=1,1,a0)
b0<-ifelse(b0<=-1,-1,b0)
b0<-ifelse(b0>=1,1,b0)
c1<-sqrt(sum(a1^2)*sum(b1^2))
c0<-sqrt(sum(a0^2)*sum(b0^2))
w1<-n1*a1*b1/c1
w0<-n0*a0*b0/c0
u=sum(w1)*p2-sum(w0)*p1
w=c(w1,w0)
u_var=p1*p2*var(w)*(n1+n0-1)
T=u^2/u_var
mysign=sign(u)
signedStat = mysign*T
pval=1-pchisq(T,1)
res = list(stat=T, pval=pval, signedStat = signedStat)
return(res)
}
st4<-function(x1,z1,x0,z0){
# Compute p-value for the equal correlation test
# with Spearman corrletion
# based on a logistic regression model
# corresponding to two independent groups
n1<-length(x1)
n0<-length(x0)
p1<-n1/(n1+n0)
p2<-1-p1
c1<-n1^2-1
c0<-n0^2-1
rank_x1<-rank(x1,ties.method = "average")
rank_z1<-rank(z1,ties.method = "average")
rank_x0<-rank(x0,ties.method = "average")
rank_z0<-rank(z0,ties.method = "average")
w1=(n1^2-1-6*(rank_x1-rank_z1)^2)/c1
w0=(n0^2-1-6*(rank_x0-rank_z0)^2)/c0
u=sum(w1)*p2-sum(w0)*p1
w=c(w1,w0)
u_var=p2*p1*var(w)*(n1+n0-1)
T=u^2/u_var
mysign=sign(u)
signedStat = mysign*T
pval=1-pchisq(T,1)
res = list(stat=T, pval=pval, signedStat = signedStat)
return(res)
}
st5=function (x1, z1, x0, z0) {
# created on May 21, 2020
# (1) a revised st5 function with output of test statistic,
# p-value, and signed test statistic.
# The sign of the test statistic is the same as the sign
# of U statistic. Since st5 is the average of st3 and st4,
# the sign of U statistic with the larger absolute value is returned
#
# test if the corr(x1, z1) is equal to corr(x0, z0) via a robust method
n1 <- length(x1)
n0 <- length(x0)
p1 <- n1/(n1 + n0)
p2 <- 1 - p1
b = 0.2
x1_median <- median(x1)
x0_median <- median(x0)
z1_median <- median(z1)
z0_median <- median(z0)
temp1 <- sort(abs(x1 - x1_median))
omhatx <- temp1[floor((1 - b) * n1)]
if (omhatx == 0) {
omhatx = temp1[n1]
}
temp2 <- sort(abs(z1 - z1_median))
omhatz <- temp2[floor((1 - b) * n1)]
if (omhatz == 0) {
omhatz = temp2[n1]
}
a1 <- (x1 - .pbos(x1, b, omhatx))/omhatx
b1 <- (z1 - .pbos(z1, b, omhatz))/omhatz
a1 <- ifelse(a1 <= -1, -1, a1)
a1 <- ifelse(a1 >= 1, 1, a1)
b1 <- ifelse(b1 <= -1, -1, b1)
b1 <- ifelse(b1 >= 1, 1, b1)
temp3 <- sort(abs(x0 - x0_median))
omhatx <- temp3[floor((1 - b) * n0)]
if (omhatx == 0) {
omhatx = temp3[n0]
}
temp4 <- sort(abs(z0 - z0_median))
omhatz <- temp4[floor((1 - b) * n0)]
if (omhatz == 0) {
omhatz = temp4[n0]
}
a0 <- (x0 - .pbos(x0, b, omhatx))/omhatx
b0 <- (z0 - .pbos(z0, b, omhatz))/omhatz
a0 <- ifelse(a0 <= -1, -1, a0)
a0 <- ifelse(a0 >= 1, 1, a0)
b0 <- ifelse(b0 <= -1, -1, b0)
b0 <- ifelse(b0 >= 1, 1, b0)
c1 <- sqrt(sum(a1^2) * sum(b1^2))
c0 <- sqrt(sum(a0^2) * sum(b0^2))
w11 <- n1 * a1 * b1/c1
w01 <- n0 * a0 * b0/c0
u1 = sum(w11) * p2 - sum(w01) * p1
w1 = c(w11, w01)
u_var1 = p1 * p2 * var(w1) * (n1 + n0 - 1)
T1 = u1^2/u_var1
c1 <- n1^2 - 1
c0 <- n0^2 - 1
rank_x1 <- rank(x1, ties.method = "average")
rank_z1 <- rank(z1, ties.method = "average")
rank_x0 <- rank(x0, ties.method = "average")
rank_z0 <- rank(z0, ties.method = "average")
w12 = (n1^2 - 1 - 6 * (rank_x1 - rank_z1)^2)/c1
w02 = (n0^2 - 1 - 6 * (rank_x0 - rank_z0)^2)/c0
u2 = sum(w12) * p2 - sum(w02) * p1
w2 = c(w12, w02)
u_var2 = p1 * p2 * var(w2) * (n1 + n0 - 1)
T2 = u2^2/u_var2
myT = 0.5 * (T1 + T2)
pval = 1 - stats::pchisq(myT, 1)
if(T1>T2)
{
mysign=sign(u1)
} else {
mysign = sign(u2)
}
signedStat = mysign*myT
res = list(stat=myT, pval=pval, signedStat = signedStat)
return(res)
}
st6<-function(x1,z1,x0,z0){
# Compute p-value for the equal correlation test
# with another way to combine Spearman corrletion and percentage bend correlation
# based on a multiple logistic regression model
# corresponding to two independent groups
n1<-length(x1)
n0<-length(x0)
p1<-n1/(n1+n0)
p2<-1-p1
b=0.2
x1_median<-median(x1)
x0_median<-median(x0)
z1_median<-median(z1)
z0_median<-median(z0)
temp1<-sort(abs(x1- x1_median))
omhatx<-temp1[floor((1-b)*n1)]
if(omhatx==0){
omhatx=temp1[n1]
}
temp2<-sort(abs(z1- z1_median))
omhatz<-temp2[floor((1-b)*n1)]
if(omhatz==0){
omhatz=temp2[n1]
}
a1<-(x1-.pbos(x1,b,omhatx))/omhatx
b1<-(z1-.pbos(z1,b,omhatz))/omhatz
a1<-ifelse(a1<=-1,-1,a1)
a1<-ifelse(a1>=1,1,a1)
b1<-ifelse(b1<=-1,-1,b1)
b1<-ifelse(b1>=1,1,b1)
temp3<-sort(abs(x0- x0_median))
omhatx<-temp3[floor((1-b)*n0)]
if(omhatx==0){
omhatx=temp3[n0]
}
temp4<-sort(abs(z0- z0_median))
omhatz<-temp4[floor((1-b)*n0)]
if(omhatz==0){
omhatz=temp4[n0]
}
a0<-(x0-.pbos(x0,b,omhatx))/omhatx
b0<-(z0-.pbos(z0,b,omhatz))/omhatz
a0<-ifelse(a0<=-1,-1,a0)
a0<-ifelse(a0>=1,1,a0)
b0<-ifelse(b0<=-1,-1,b0)
b0<-ifelse(b0>=1,1,b0)
c1<-sqrt(sum(a1^2)*sum(b1^2))
c0<-sqrt(sum(a0^2)*sum(b0^2))
w11<-n1*a1*b1/c1
w01<-n0*a0*b0/c0
u1=sum(w11)*p2-sum(w01)*p1
w1=c(w11,w01)
c1<-n1^2-1
c0<-n0^2-1
rank_x1<-rank(x1,ties.method = "average")
rank_z1<-rank(z1,ties.method = "average")
rank_x0<-rank(x0,ties.method = "average")
rank_z0<-rank(z0,ties.method = "average")
w12=(n1^2-1-6*(rank_x1-rank_z1)^2)/c1
w02=(n0^2-1-6*(rank_x0-rank_z0)^2)/c0
u2=sum(w12)*p2-sum(w02)*p1
w2=c(w12,w02)
w1_mean=mean(w1)
w2_mean=mean(w2)
a11=sum((w1-w1_mean)^2)
a12=sum((w1-w1_mean)*(w2-w2_mean))
a21=a12
a22=sum((w2-w2_mean)^2)
u=matrix(c(u1,u2),nrow = 1,ncol = 2,byrow = TRUE)
cov_u=p1*p2*matrix(c(a11,a12,a21,a22),nrow = 2,ncol = 2,byrow = TRUE)
st6<-tryCatch( {1-pchisq(u%*%solve(cov_u)%*%t(u),2)},error=function(e){
1-pchisq(u%*%ginv(cov_u)%*%t(u),1)
})
res = list(u=u, pval=st6, cov_u=cov_u)
return(res)
}
construct_network<-function(es,
cor_method = 'st5',
var.grp,
pseudo_adjust_cutoff = FALSE,
pAdjMethod = 'fdr',
cutoff = 0.05,
nPseudo = 25){
# input:
# es: an ExpressionSet object of microRNA dataset
# cor_method: The method for equal correlation, 'ST5' is recommand.
# var.grp: phenotype variable name indicating case-control status,0 as control, 1 as case.
# pAdjMethod:if pAdjMethod=FALSE, the function will not do mutiple testing adjustment. If pAdjMethod="fdr"/"BH"/"BY"/"holm"/"hochberg"/"hommel"/"bonferroni"/"BH"/"BY", the specific method will be used for adjusting p-value.
# cutoff: if p value is smaller than the cutoff, there will be an edge between the two nodes.
# pseudo_adjust_cutoff: if the value is true, pseudo probes will be used for adjusting the cutoff
# output:
# my_graph: obtained network in igraph object
# my_dat: obtained netork as data frame with 3 columns: edge id, node_id1,node_id2
mat1=Biobase::exprs(es)
pDat1=Biobase::pData(es)
fDat1=Biobase::fData(es)
geneid=featureNames(es)
status = pDat1[, c(var.grp)]
pos0=which(status==0) # controls
pos1=which(status==1) # cases
nProbes=dim(mat1)[1] # number of genes/probes
probes=1:nProbes
r1_pos=c()
r2_pos=c()
r1=c()
# for each probe, calculate test statistic between this probe and other probes
for(i in 1:(nProbes-1)){
temp1=rep(probes[i],time=nProbes-i)
temp2=probes[(i+1):nProbes]
r1_pos=c(r1_pos,temp1)
r2_pos=c(r2_pos,temp2)
x11=mat1[temp1,pos1]
x10=mat1[temp1,pos0]
z11=mat1[temp2,pos1]
z10=mat1[temp2,pos0]
if(i!=(nProbes-1)){
x11<-tapply(x11,rep(1:base::nrow(x11),base::ncol(x11)),function(i)i)
x10<-tapply(x10,rep(1:base::nrow(x10),base::ncol(x10)),function(i)i)
z11<-tapply(z11,rep(1:base::nrow(z11),base::ncol(z11)),function(i)i)
z10<-tapply(z10,rep(1:base::nrow(z10),base::ncol(z10)),function(i)i)
switch (cor_method,
st1 = {temp1=as.numeric(mapply(st1,x11,z11,x10,z10)[2,])},
st2 = {temp1=as.numeric(mapply(st2,x11,z11,x10,z10)[2,])},
st3 = {temp1=as.numeric(mapply(st3,x11,z11,x10,z10)[2,])},
st4 = {temp1=as.numeric(mapply(st4,x11,z11,x10,z10)[2,])},
st5 = {temp1=as.numeric(mapply(st5,x11,z11,x10,z10)[2,])},
st6 = {temp1=as.numeric(mapply(st6,x11,z11,x10,z10)[2,])},
stop("wrong cor_method!")
)
}else{
switch(cor_method,
st1={temp1=st1(x11,z11,x10,z10)$pval},
st2={temp1=st2(x11,z11,x10,z10)$pval},
st3={temp1=st3(x11,z11,x10,z10)$pval},
st4={temp1=st4(x11,z11,x10,z10)$pval},
st5={temp1=st5(x11,z11,x10,z10)$pval},
st6={temp1=st6(x11,z11,x10,z10)$pval},
stop('wrong cor_method')
)
}
r1=c(r1,temp1)
}
# p-values for pseudo genes
pvalPseudo = NULL
# decide if a pair of probes is differentially correlated
if(pseudo_adjust_cutoff==TRUE){
#nr=nrow(es) # number of probes
if(nPseudo >= nProbes/2)
{
stop("nPseudo must be smaller than half of number of genes!")
}
#nPseudo=50
# randomly pick 'nPseudo' genes as pseudo genes
pos.pick=sample(1:nProbes, size=nPseudo, replace=FALSE)
mat1.pick=mat1[pos.pick,]
fDat1.pick=fDat1[pos.pick,]
nSubj=ncol(es) # number of subjects
# for each pseudo gene gene, permute subject ids
pos=sample(1:nSubj, size=nSubj, replace=FALSE)
mat1.pick.s=mat1.pick[1,pos, drop=FALSE]
for(i in 2:nPseudo){
pos=sample(1:nSubj, size=nSubj, replace=FALSE)
mat1.pick.s=rbind(mat1.pick.s,mat1.pick[i,pos])
}
# make sure to change the row names
rownames(mat1.pick.s)=paste("perm", rownames(mat1.pick), sep=".")
rownames(fDat1.pick)=rownames(mat1.pick.s)
# also need to change the content of fDat.pick
# suppose column 'ID' is the probe id
fDat1.pick$ID=paste("perm", fDat1.pick$ID, sep=".")
# add dat.pick.s to dat
# mat2 include both original p1+p2 genes and nPseudo pseudo genes
mat2=rbind(mat1,mat1.pick.s)
fDat1$ID = NA
fDat2=rbind(fDat1,fDat1.pick)
#nrow=dim(mat1)[1] # number of genes
#nPseudo=50
#cutoff1=c()
# for each pseudo gene, we obtain its test statistic for differential correlation test
for(i in (nProbes-nPseudo+1):nProbes){
# non-pseudo genes in cases
x11=mat1[c(1:(nProbes-nPseudo)),pos1]
# non-pseudo genes in controls
x10=mat1[c(1:(nProbes-nPseudo)),pos0]
# i-the pseudo genes in cases
z11=mat1[rep(i,time=(nProbes-nPseudo)),pos1]
# i-th pseudo genes in controls
z10=mat1[rep(i,time=(nProbes-nPseudo)),pos0]
x11<-tapply(x11,rep(1:base::nrow(x11),base::ncol(x11)),function(i)i)
x10<-tapply(x10,rep(1:base::nrow(x10),base::ncol(x10)),function(i)i)
z11<-tapply(z11,rep(1:base::nrow(z11),base::ncol(z11)),function(i)i)
z10<-tapply(z10,rep(1:base::nrow(z10),base::ncol(z10)),function(i)i)
# obtain p-value for testing differential correlation
switch (cor_method,
st1 = {temp1=as.numeric(mapply(st1,x11,z11,x10,z10)[2,])},
st2 = {temp1=as.numeric(mapply(st2,x11,z11,x10,z10)[2,])},
st3 = {temp1=as.numeric(mapply(st3,x11,z11,x10,z10)[2,])},
st4 = {temp1=as.numeric(mapply(st4,x11,z11,x10,z10)[2,])},
st5 = {temp1=as.numeric(mapply(st5,x11,z11,x10,z10)[2,])},
st6 = {temp1=as.numeric(mapply(st6,x11,z11,x10,z10)[2,])},
stop("wrong cor_method!")
)
# we expect the p-values obtained based on pseudo genes are large
# because pseudo genes are non-differentially correlated.
# 'cutoff' percentile of p-value for testing differential correlation
#pvalPseudo=c(pvalPseudo, stats::quantile(temp1, prob=cutoff))
pvalPseudo=c(pvalPseudo, temp1)
}
#alpha1=stats::quantile(pvalPseudo, prob=cutoff)
alpha1 = min(pvalPseudo, na.rm=TRUE)
r1Adj = r1
}else{
alpha1=cutoff
r1Adj=stats::p.adjust(r1, pAdjMethod)
}
adjacency_matrix1=matrix(0,nrow = nProbes, ncol = nProbes)
rownames(adjacency_matrix1) = geneid
colnames(adjacency_matrix1) = geneid
pvalMat=matrix(0,nrow = nProbes, ncol = nProbes)
rownames(pvalMat) = geneid
colnames(pvalMat) = geneid
pAdjMat=matrix(0,nrow = nProbes, ncol = nProbes)
rownames(pAdjMat) = geneid
colnames(pAdjMat) = geneid
j=1
my_dat=data.frame()
# if p-value < alpha1, then the 2 genes are differentially correlated
for(i in 1:length(r1)){
if(r1Adj[i]<alpha1){
adjacency_matrix1[r1_pos[i],r2_pos[i]]=1
adjacency_matrix1[r2_pos[i],r1_pos[i]]=1
my_dat=rbind(my_dat,data.frame(edge_id=j,node_id1=r1_pos[i],node_id2=r2_pos[i]))
j=j+1
}
pvalMat[r1_pos[i],r2_pos[i]] = r1[i]
pvalMat[r2_pos[i],r1_pos[i]] = r1[i]
pAdjMat[r1_pos[i],r2_pos[i]] = r1Adj[i]
pAdjMat[r2_pos[i],r1_pos[i]] = r1Adj[i]
}
my_graph=igraph::graph_from_adjacency_matrix(as.matrix(adjacency_matrix1),
mode = "undirected",
weighted = NULL)
res = list(graph = my_graph,
network_dat = my_dat,
pvalMat = pvalMat,
pAdjMat = pAdjMat,
pvalPseudo = pvalPseudo,
alpha1 = alpha1)
invisible(res)
}
generate_data<-function(n1=50,n2=60,p1=5,p2=100){
### The function is to generate expression level matrixes of control subjects and case subjects.
### X matrix is for diseased subjects with the default sample size n1=50. Z matrix is for control subjects with the default sample size n2=60.
### X is generated from multivariate normal distribution N (0, SigmaX), where SigmaX is a block matrix ((SigmaP1, 0), (0, SigmaP2)), sigmaP1 is the p1*p1 matrix and SigmaP2 is the p2*p2 matrix. Z is generated from multivariate normal distribution N (0, SigmaZ), where SigmaZ is a block matrix ((E_P1, 0), (0, SigmaP2)) and E_P1 is p1*p1 identity matrix.
#n1=50
#n2=60
#p1=5
#p2=100
Sigma1=clusterGeneration::rcorrmatrix(d=p1)
Sigma0=clusterGeneration::rcorrmatrix(d=p2)
Ip1=diag(p1)
SigmaX=as.matrix(Matrix::bdiag(Sigma1,Sigma0))
SigmaZ=as.matrix(Matrix::bdiag(Ip1,Sigma0))
# X - n1 x (p1+p2) matrix (cases)
X=MASS::mvrnorm(n = n1, mu=numeric(p1+p2), Sigma=SigmaX, tol = 1e-6, empirical = FALSE)
# Z - n2 x (p1+p2) matrix (controls)
Z=MASS::mvrnorm(n = n2, mu=numeric(p1+p2), Sigma=SigmaZ, tol = 1e-6, empirical = FALSE)
# rows are genes
# columns are subjects
dat = t(rbind(X, Z))
geneid=paste("g", 1:nrow(dat), sep="")
sid=paste("s", 1:ncol(dat), sep="")
rownames(dat)=geneid
colnames(dat)=sid
# grp=0 indicates control; grp=1 indicates case
grp=c(rep(1, n1), rep(0, n2))
pDat = data.frame(sid=sid, grp=grp)
rownames(pDat) = sid
# memGenes=1 indicates the gene is differentially correlated with at least one another gene between cases and controls
# that is, correlation of this gene with at least one other gene in cases is different from that in controls.
# memGenes=0 indicates the gene is non-differentially correlated with any other genes between cases and controls
# that is, correlation of this gene with other genes in cases is the same as that in controls.
memGenes=c(rep(1, p1), rep(0, p2))
fDat = data.frame(geneid=geneid, memGenes=memGenes)
rownames(fDat)=geneid
es = genEset(ex=dat, pDat=pDat, fDat = fDat, annotation = "")
rownames(SigmaX)=geneid
colnames(SigmaX)=geneid
rownames(SigmaZ)=geneid
colnames(SigmaZ)=geneid
res = list(es=es, covCase = SigmaX, covCtrl = SigmaZ)
invisible(res)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.