ratio {BGcom} | R Documentation |
This function take the lists of pvalues and calculate the ratio T(q) for each threshold q
ratio(n, dir, data1, data2, interval, name, pvalue)
n |
length of the vector of simulated pvalue |
dir |
directory for storing the plots |
data1,data2 |
Lists of pvalues to be compared |
interval |
The interval between two threshold |
name |
The name to be used in the plots |
pvalue |
Indicate if the data are pvalues (TRUE) or posterior probability (FALSE). If they are posterior probability they are transformed in pvalues automatically |
This function calculate the ratio T(q) of observed number of genes in common vs expected one for each threshold q.
Returns an object of class list with the ratio and the threshold. In particular:
m |
maximum T |
r |
threshold associated to the max T |
common |
number of features corresponding to the threshold r |
l1 |
differentially expressed features in the first experiments corresponding to the threshold r |
l2 |
differentially expressed features in the second experiments corresponding to the threshold r |
ratios |
vector or T values for each threshold |
thresh.ratios |
threshold corresponding to T values |
L1 |
differentially expressed features in the first experiment corresponding to the T values |
L2 |
differentially expressed features in the second experiment corresponding to the T values |
int |
features in common corresponding to the T values |
interval |
interval on the p-value scale defined by the user (default is 0.01) |
name |
names to be used in the plots (defined by the user) |
M. Blangiardo
Stone et al.(1988), Investigations of excess environmental risks around putative sources: statistical problems and a proposed test,Statistics in Medicine, 7, 649-660. M.Blangiardo and S.Richardson Statistical tools for synthesizing lists of differentially expressed features in related experiments (under revision)
data = simulation(n=2000,GammaA=1,GammaB=1,r1=0.5,r2=0.8,DEfirst=1000,DEsecond=800,DEcommon=700) T<- ratio(n=length(data$Pval1),data1=data$Pval1,data2=data$Pval2,interval=0.01,dir="D:/",name="CompData1Data2",pvalue=FALSE) ## The function is currently defined as function(n,dir,data1,data2,interval,name,pvalue){ if(pvalue==TRUE){ data1=data1; data2=data2 } if(pvalue==FALSE){ data1=1-data1; data2=1-data2 } ID=seq(1,n) l=length(seq(0,1,interval)) threshold = cbind(rep(NA,l),seq(0,1,interval)) int=c() L1=c() L2=c() for(i in 1:l){ threshold[i,1] <- (length(ID[data1<=threshold[i,2] & data2<=threshold[i,2]])) int[i]<- threshold[i,1] } threshold.expected = cbind(rep(NA,l),rep(NA,l),rep(NA,l),seq(0,1,interval)) for(i in 1:l){ threshold.expected[i,1] <- length(ID[data1<=threshold.expected[i,4]])/n; threshold.expected[i,2]<- length(ID[data2<=threshold.expected[i,4]])/n; threshold.expected[i,3] <- threshold.expected[i,1]*threshold.expected[i,2]*n L1[i] <- threshold.expected[i,1]*n L2[i] <- threshold.expected[i,2]*n } #Calculate the ratio for the number of observed genes/number of expected ones expected = threshold.expected[,3] observed = threshold[,1] ratios = matrix(0,l,1) for(i in 1:l){ ratios[i,1] <- observed[i]/expected[i] } ratios=ratios[int>0,1] thresh.ratios=threshold[int>0,2] L1=L1[int>0] L2=L2[int>0] int=int[int>0] ps.options(paper="a4",horizontal=TRUE) setwd(dir) ps.options(horizontal=FALSE) postscript(paste("Ratio","_",name,".ps")) plot(thresh.ratios,ratios,type="l", ylab= "T",xlab="Pvalue",main="",yaxt="n",cex.main=0.7,cex.axis=1.2,ylim=c(0,(max(ratios,na.rm=TRUE)+sd(ratios,na.rm=TRUE)))) axis(2, at = c(seq(0,max(ratios,na.rm=TRUE)+sd(ratios,na.rm=TRUE),0.5),max(ratios,na.rm=TRUE)), labels = c(seq(0,max(ratios,na.rm=TRUE)+sd(ratios,na.rm=TRUE),0.5),round(max(ratios,na.rm=TRUE),1)), tick = TRUE,cex.axis=1.2) m=max(ratios,na.rm=TRUE) r=thresh.ratios[ratios==m] if(length(r)>1) r=round(mean(r),2) common = threshold[threshold[,2]==r,1] l1 = threshold.expected[threshold.expected[,4]==r,1]*n l2 = threshold.expected[threshold.expected[,4]==r,2]*n legend(x=0.3,y=m/6,legend= expression(paste(O["1+"],"(q*)","=")),bty="n") legend(x=0.3,y=m/6*2,legend= expression(paste(O["+1"],"(q*)","=")),bty="n") legend(x=0.3,y=m/6*4,legend= expression(paste(O[11],"(q*)","=")),bty="n") legend(x=0.3,y=m/6*5,legend= expression(paste("q*","=")),bty="n") legend(x=0.3,y=m/6*6,legend= expression(paste(T("q*"),"=")),bty="n") legend(x=0.45,y=m/6,legend= l1,bty="n") legend(x=0.45,y=m/6*2,legend= l2,bty="n") legend(x=0.45,y=m/6*4,legend= common,bty="n") legend(x=0.4,y=m/6*5,legend= r,bty="n") legend(x=0.4,y=m/6*6,legend= round(m,2),bty="n") dev.off() return(list(m=m,r=r,common=common,l1=l1,l2=l2,ratios=ratios,thresh.ratios=thresh.ratios,L1=L1,L2=L2,int=int,interval=interval,name=name)) }