ratio {BGcom}R Documentation

Function for calculating the ratio T between the observed genes in common and the expected one

Description

This function take the lists of pvalues and calculate the ratio T(q) for each threshold q

Usage

ratio(n, dir, data1, data2, interval, name, pvalue)

Arguments

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

Details

This function calculate the ratio T(q) of observed number of genes in common vs expected one for each threshold q.

Value

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)

Author(s)

M. Blangiardo

References

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)

Examples

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))
        }

[Package BGcom version 1.0 Index]