2017-08-28 63 views
1

後在安裝CATS包(2013年)和R的變化版本,並沒有什麼幾次失敗的嘗試,我決定從here調試 - 功率計算器功能和情節

的源代碼工作,我創建了一個單R從包中所有的R功能的腳本,我 跑他們,然後我就希望這個代碼將運行下面的代碼結束繪製權力的例子:

super.cats<-function(RR,MAFmax=0.5,MAFmin=0.005,by=50,rep=1536,SNPs=1E6,ncases,ncontrols,ncases2,ncontrols2,alpha=0.05/SNPs,...){ 
    powerList.O<-c() 
    powerList.J<-c() 
    powerList.R<-c() 
    powerList.F<-c() 
    power.O<-rep(0,length(RR)) 
    power.F<-rep(0,length(RR)) 
    power.J<-rep(0,length(RR)) 
    power.R<-rep(0,length(RR)) 

    MAF<-exp(seq(log(MAFmin),log(MAFmax),by=(log(MAFmax)-log(MAFmin))/by)) 
    for(nmaf in 1:length(MAF)){ 
    for(tal in 1:length(RR)){ 
     if(power.F[tal]>0.995&power.R[tal]>0.995){ 
     power.O[tal]<-1 
     power.R[tal]<-1 
     power.J[tal]<-1 
     power.F[tal]<-1 
     break 
     } 
     temp<-cats(risk=RR[tal],freq=MAF[nmaf],ncases=ncases,ncontrols=ncontrols,ncases2=ncases2,ncontrols2=ncontrols2,alpha=alpha,pimarkers=rep/SNPs,...) 
     power.O[tal]<-temp$P.one.study 
     power.J[tal]<-temp$P.joint 
     power.R[tal]<-temp$P.rep.study 
     power.F[tal]<-temp$P.first.stage 

    } 
    powerList.O<-cbind(powerList.O,power.O) 
    powerList.J<-cbind(powerList.J,power.J) 
    powerList.R<-cbind(powerList.R,power.R) 
    powerList.F<-cbind(powerList.F,power.F) 

    cat(nmaf," ") 
    } 
    cat("\n") 

    obs<-list(powerList.O=powerList.O,powerList.J=powerList.J,powerList.R=powerList.R,powerList.F=powerList.F,RR=RR,MAF=MAF,ncases=ncases,ncontrols=ncontrols,ncases2=ncases2,ncontrols2=ncontrols2,rep=rep,curve=F) 

    class(obs)<-"supercats" 
    return(obs) 
} 


############################## 

if(FALSE){ 
    #heat plot 
    rr<-seq(1,2,by=0.025) 
    c<-super.cats(rr,by=length(rr),ncases=765,ncontrols=1274,ncases2=100,ncontrols2=100,alpha=0.001,prevalence=0.01); 
    plot(c,main="power",file=NULL) 


    #curves 
    rr<-seq(1,3,by=0.05) 
    maf<-c(0.01,0.05,0.2,0.5) 
    c2<-curve.cats(rr,maf,ncases=765,ncontrols=1274,ncases2=100,ncontrols2=100,alpha=0.001,prevalence=0.01); 
    plot(c2,main="power2",ylab="Power",xlab="RR",file=NULL,col=1:4) 


} 


#### 

"cats" <- 
    function (freq=0.5,freq2=-1,ncases=500,ncontrols=500,ncases2=500,ncontrols2=500,risk=1.5,risk2=-1,pisamples=-1,prevalence=0.1,prevalence2=-1,additive=0,recessive=0,dominant=0,multiplicative=1,alpha=0.0000001,pimarkers=0.00316) 
    { model<-c(additive,recessive,dominant,multiplicative) 

    if(sum(model==1)!=1) 
    stop("chose only one model. e.i. one model must be 1 the others 0") 
    if(sum(model==0)!=3) 
    stop("chose only one model. e.i. one model must be 1 the others 0") 


    if(freq<0|freq>1) 
    stop("freq must be between 0 and 1") 
    if((freq2<0|freq2>1)&freq2!=-1) 
    stop("freq2 must be between 0 and 1 (or undefined as -1)") 
    if((pisamples<0|pisamples>1)&pisamples!=-1) 
    stop("pisamples must be between 0 and 1") 
    if((prevalence2<0|prevalence2>1)&prevalence2!=-1) 
    stop("prevalence2 must be between 0 and 1 (or undefined as -1)") 
    if(alpha<0|alpha>1) 
    stop("alpha must be between 0 and 1") 
    if(prevalence<0|prevalence>1) 
    stop("prevalence must be between 0 and 1") 
    if(pimarkers<0|pimarkers>1) 
    stop("pimarkers must be between 0 and 1") 
    if(ncases!=as.integer(ncases)|ncases<0) 
    stop("ncases must be a positive integer") 
    if(ncases2!=as.integer(ncases2)|ncases2<0) 
    stop("ncases2 must be a positive integer") 
    if(ncontrols!=as.integer(ncontrols)|ncontrols<0) 
    stop("ncontrols must be a positive integer") 
    if(ncontrols2!=as.integer(ncontrols2)|ncontrols2<0) 
    stop("ncontrols2 must be a positive integer") 
    if(risk<0) 
    stop("risk must be positive") 
    if(risk2<0&risk2!=-1) 
    stop("risk2 must be positive(or undefined as -1)") 



    res<-.Call("cats", 
      as.double(freq),as.double(freq2),as.integer(ncases),as.integer(ncontrols), 
      as.integer(ncases2),as.integer(ncontrols2),as.double(risk),as.double(risk2), 
      as.double(pisamples),as.double(prevalence),as.double(prevalence2), 
      as.integer(additive),as.integer(recessive),as.integer(dominant), 
      as.integer(multiplicative),as.double(alpha),as.double(pimarkers)) 



    options<-cbind(freq,freq2,ncases=ncases,ncontrols=ncontrols,ncases2=ncases2,ncontrols2=ncontrols2,risk,risk2,pisamples,prevalence,prevalence2,additive,recessive,dominant,multiplicative,alpha,pimarkers) 

    result<-list(P.one.study=res[1,1],P.first.stage=res[2,1],P.rep.study=res[3,1],P.joint.min=res[4,1],P.joint=res[5,1],pi=res[6,1],T.one.study=res[7,1],T.first.stage=res[8,1],T.second.stage.rep=res[9,1],T.second.stage.joint=res[10,1],E.Disease.freq.cases1=res[11,1],E.Disease.freq.controls1=res[12,1],E.Disease.freq.cases2=res[13,1],E.Disease.freq.controls2=res[14,1],options=options) 
    class(result)<-"CATS" 
    return(result) 
    } 

curve.cats<-function(RR,MAF,rep=1536,SNPs=1E6,ncases,ncontrols,ncases2,ncontrols2,alpha=0.05/SNPs,...){ 
    powerList.O<-c() 
    powerList.J<-c() 
    powerList.R<-c() 
    powerList.F<-c() 
    power.O<-rep(0,length(RR)) 
    power.F<-rep(0,length(RR)) 
    power.J<-rep(0,length(RR)) 
    power.R<-rep(0,length(RR)) 
    for(nmaf in 1:length(MAF)){ 
    for(tal in 1:length(RR)){ 
     if(power.F[tal]>0.995&power.R[tal]>0.995){ 
     power.O[tal]<-1 
     power.R[tal]<-1 
     power.J[tal]<-1 
     power.F[tal]<-1 
     break 
     } 
     tempo<-cats(risk=RR[tal],freq=MAF[nmaf],ncases=ncases,ncontrols=ncontrols,ncases2=ncases2,ncontrols2=ncontrols2,alpha=alpha,pimarkers=rep/SNPs,...) 
     power.O[tal]<-tempo$Pone.study 
     power.J[tal]<-tempo$Pjoint 
     power.R[tal]<-tempo$Prep.study 
     power.F[tal]<-tempo$Pfirst.stage 
    } 
    powerList.O<-cbind(powerList.O,power.O) 
    powerList.J<-cbind(powerList.J,power.J) 
    powerList.R<-cbind(powerList.R,power.R) 
    powerList.F<-cbind(powerList.F,power.F) 
    cat(nmaf," ") 
    } 
    cat("\n") 

    obs<-list(powerList.O=powerList.O,powerList.J=powerList.J,powerList.R=powerList.R,powerList.F=powerList.F,RR=RR,MAF=MAF,ncases=ncases,ncontrols=ncontrols,ncases2=ncases2,ncontrols2=ncontrols2,rep=rep,curve=T) 

    class(obs)<-"supercats" 
    return(obs) 
} 


lines.cats<-function(x,type="Replication",col=NULL,lty=2,...){ 
    if(type=="Joint") 
    power<-x$powerList.J 
    else if(type=="One") 
    power<-x$powerList.O 
    else if(type=="Replication") 
    power<-x$powerList.R 
    else if(type=="First") 
    power<-x$powerList.F 

    if(x$curve){ 
    if(is.null(col)) 
     col=1:length(x$MAF) 
    for(nmaf in 1:length(x$MAF)) 
     lines(x$RR,power[,nmaf],col=col[nmaf],lwd=2,lty=lty) 

    } 
    else 
    cat("only for curves \n") 
} 

rr <- seq(1,2,by=0.05) 
maf <- c(0.05,0.1,0.2,0.5) 
c2 <- curve.cats(rr,maf,ncases=600,ncontrols=600,ncases2=600,ncontrols2=600, alpha=0.000001,prevalence=0.01); 

plot(c2,type="One",main="power2",ylab="Power",xlab="RR",file=NULL,col=1:4) 
lines.cats(c2,type="Replication",lty=3) 
lines.cats(c2,type="Joint",lty=2) 
lines.cats(c2,type="First",lty=4) 
legend("left",c("One stage","Joint","Relication","First Stage"),lty=1:4,bty="n") 


### 


lines.cats<-function(x,type="Replication",col=NULL,lty=2,...){ 
    if(type=="Joint") 
    power<-x$powerList.J 
    else if(type=="One") 
    power<-x$powerList.O 
    else if(type=="Replication") 
    power<-x$powerList.R 
    else if(type=="First") 
    power<-x$powerList.F 

    if(x$curve){ 
    if(is.null(col)) 
     col=1:length(x$MAF) 
    for(nmaf in 1:length(x$MAF)) 
     lines(x$RR,power[,nmaf],col=col[nmaf],lwd=2,lty=lty) 

    } 
    else 
    cat("only for curves \n") 
} 


#### 


plot.supercats<-function(x,type="Joint",file="power.pdf",col=NULL,main=paste("POWER N=",x$ncases,":",x$ncontrols,",",x$ncases2,":",x$ncontrols2," rep=",x$rep,sep=""),...){ 
    if(type=="Joint") 
    power<-x$powerList.J 
    else if(type=="One") 
    power<-x$powerList.O 
    else if(type=="Replication") 
    power<-x$powerList.R 
    else if(type=="First") 
    power<-x$powerList.F 


    if(!is.null(file)) 
    pdf(file) 
    #curve 
    if(x$curve){ 
    if(is.null(col)) 
     col=1:length(x$MAF) 
    plot(x$RR,power[,1],ylim=c(0,1),main=main,col="transparent",...) 
    for(nmaf in 1:length(x$MAF)){ 
     lines(x$RR,power[,nmaf],col=col[nmaf],lwd=2) 
    } 
    legend(min(x$RR),1,paste("MAF=",x$MAF),col=col,lwd=2,bty="n") 
    } 
    else{ 
    #image 
    if(is.null(col)) 
     col=heat.colors(80) 

    image(x$RR,x$MAF,power,col=col,main=main,log="y",ylim=c(0.005,.5),ylab="MAF",xlab="RR",...) 
    legend("topright",paste(1:10*10,"%"),fill=col[1:10*8],bty="n") 

    } 
    if(!is.null(file)) 
    dev.off() 
} 


#### 
.onLoad=function(libname, pkgname) 
{ 
    library.dynam("CATS", pkgname, libname) 
} 
.onUnload=function(libpath) 
{ 
    library.dynam.unload("CATS", libpath) 
} 



#### 

"summary.CATS" <- 
    function(object, ...){ 
    if (!inherits(object, "CATS")) 
     stop("Not an object of class CATS!") 
    cat("Options \n") 
    ob<-t(object$options) 
    colnames(ob)<-"chosen" 
    print(ob) 

    cat("Recommended thresholds:") 
    print(structure(list("One stage Design"=object$T.one.study,"Stage 1 Threshold"=object$T.first.stage,"Replication Threshold"=object$T.second.stage.rep,"Joint Analysis Threshold"=object$T.second.stage.joint),class="power.htest")) 

    cat("Eobjectpected disesase allele frequencies") 
    print(structure(list("Cases in stage 1"=object$E.Disease.freq.cases1,"Controls in stage 1 "=object$E.Disease.freq.controls1,"Cases in stage 2"=object$E.Disease.freq.cases2,"Controls in stage 2"=object$E.Disease.freq.controls2),class="power.htest")) 

    cat("Expected Power is:") 
    print(structure(list("For a one-stage study" = signif(object$P.one.study, 
                  3), "For first stage in two-stage study" = signif(object$P.first.stage, 
                              3), "For second stage in replication analysis" = signif(object$P.rep.study, 
                                            3), "For second stage in a joint analysis" = signif(object$P.joint, 
                                                         3), pi = signif(object$pi, 3)), class = "power.htest")) 
    } 


### 


"print.CATS" <- 
    function(x, ...){ 
    if(!inherits(x,"CATS")) 
     stop("Not an object of class CATS!") 
    cat("Expected Power is;\n") 
    print(structure(list("For a one-stage study"=signif(x$P.one.study,3),"For first stage in two-stage study"=signif(x$P.first.stage,3),"For second stage in replication analysis"=signif(x$P.rep.study,3),"For second stage in a joint analysis"=signif(x$P.joint,3),"pi"=signif(x$pi,3)),class="power.htest")) 
    } 

### 

"cats" <- 
    function (freq=0.5,freq2=-1,ncases=500,ncontrols=500,ncases2=500,ncontrols2=500,risk=1.5,risk2=-1,pisamples=-1,prevalence=0.1,prevalence2=-1,additive=0,recessive=0,dominant=0,multiplicative=1,alpha=0.0000001,pimarkers=0.00316) 
    { 


    model<-c(additive,recessive,dominant,multiplicative) 

    if(sum(model==1)!=1) 
     stop("chose only one model. e.i. one model must be 1 the others 0") 
    if(sum(model==0)!=3) 
     stop("chose only one model. e.i. one model must be 1 the others 0") 


    if(freq<0|freq>1) 
     stop("freq must be between 0 and 1") 
    if((freq2<0|freq2>1)&freq2!=-1) 
     stop("freq2 must be between 0 and 1 (or undefined as -1)") 
    if((pisamples<0|pisamples>1)&pisamples!=-1) 
     stop("pisamples must be between 0 and 1") 
    if((prevalence2<0|prevalence2>1)&prevalence2!=-1) 
     stop("prevalence2 must be between 0 and 1 (or undefined as -1)") 
    if(alpha<0|alpha>1) 
     stop("alpha must be between 0 and 1") 
    if(prevalence<0|prevalence>1) 
     stop("prevalence must be between 0 and 1") 
    if(pimarkers<0|pimarkers>1) 
     stop("pimarkers must be between 0 and 1") 
    if(ncases!=as.integer(ncases)|ncases<0) 
     stop("ncases must be a positive integer") 
    if(ncases2!=as.integer(ncases2)|ncases2<0) 
     stop("ncases2 must be a positive integer") 
    if(ncontrols!=as.integer(ncontrols)|ncontrols<0) 
     stop("ncontrols must be a positive integer") 
    if(ncontrols2!=as.integer(ncontrols2)|ncontrols2<0) 
     stop("ncontrols2 must be a positive integer") 
    if(risk<0) 
     stop("risk must be positive") 
    if(risk2<0&risk2!=-1) 
     stop("risk2 must be positive(or undefined as -1)") 



    res<-.Call("cats", 
       as.double(freq),as.double(freq2),as.integer(ncases),as.integer(ncontrols), 
       as.integer(ncases2),as.integer(ncontrols2),as.double(risk),as.double(risk2), 
       as.double(pisamples),as.double(prevalence),as.double(prevalence2), 
       as.integer(additive),as.integer(recessive),as.integer(dominant), 
       as.integer(multiplicative),as.double(alpha),as.double(pimarkers),PACKAGE="CATS") 



    options<-cbind(freq,freq2,ncases=ncases,ncontrols=ncontrols,ncases2=ncases2,ncontrols2=ncontrols2,risk,risk2,pisamples,prevalence,prevalence2,additive,recessive,dominant,multiplicative,alpha,pimarkers) 

    result<-list(P.one.study=res[1,1],P.first.stage=res[2,1],P.rep.study=res[3,1],P.joint.min=res[4,1],P.joint=res[5,1],pi=res[6,1],T.one.study=res[7,1],T.first.stage=res[8,1],T.second.stage.rep=res[9,1],T.second.stage.joint=res[10,1],E.Disease.freq.cases1=res[11,1],E.Disease.freq.controls1=res[12,1],E.Disease.freq.cases2=res[13,1],E.Disease.freq.controls2=res[14,1],options=options) 
    class(result)<-"CATS" 
    return(result) 
    } 

#### EXAMPLE 
rr<-seq(1,2,by=0.05) 
maf<-c(0.05,0.1,0.2,0.5) 
c2<-curve.cats(rr,maf,ncases=600,ncontrols=600,ncases2=600, 
       ncontrols2=600,alpha=0.000001,prevalence=0.01) 

plot(c2,type="One",main="power2",ylab="Power",xlab="RR",file=NULL,col=1:4) 
lines.cats(c2,type="Replication",lty=3) 
lines.cats(c2,type="Joint",lty=2) 
lines.cats(c2,type="First",lty=4) 
legend("left",c("One stage","Joint","Relication","First Stage"),lty=1:4,bty="n") 

但我得到以下錯誤:

Error in .Call("cats", as.double(freq), as.double(freq2), as.integer(ncases), : "cats" not available for .Call() for package "CATS" Called from: cats(risk = RR[tal], freq = MAF[nmaf], ncases = ncases, ncontrols = ncontrols, 
    ncases2 = ncases2, ncontrols2 = ncontrols2, alpha = alpha, 
    pimarkers = rep/SNPs, ...) 

我試圖修改代碼,但我更改它,更多的錯誤出現。在這一點上,我將不勝感激任何形式的幫助。

更新對來自R安裝包時,我得到什麼:

adris-imac:Desktop gwallace$ R CMD INSTALL CATS_1.02.tar.gz 
* installing to library ‘/Library/Frameworks/R.framework/Versions/3.4/Resources/library’ 
* installing *source* package ‘CATS’ ... 

** libs clang -I/Library/Frameworks/R.framework/Resources/include -DNDEBUG -I/usr/local/include -fPIC -Wall -g -O2 -c CATS.c -o CATS.o In file included from CATS.c:4: ./cats.h:196:27: warning: '&&' within '||' [-Wlogical-op-parentheses] if (z > LOWER_TAIL_ONE && !upper || z > UPPER_TAIL_ZERO) 
     ~~~~~~~~~~~~~~~~~~~^~~~~~~~~ ~~ ./cats.h:196:27: note: place parentheses around the '&&' expression to silence 
     this warning if (z > LOWER_TAIL_ONE && !upper || z > UPPER_TAIL_ZERO) 
         ^
     (       ) CATS.c:86:7: error: non-void function 'cats' should return a value 
     [-Wreturn-type] 
     return ; 
    ^CATS.c:106:7: error: non-void function 'cats' should return a value 
     [-Wreturn-type] 
     return ; 
    ^CATS.c:133:7: error: non-void function 'cats' should return a value 
     [-Wreturn-type] 
     return ; 
    ^1 warning and 3 errors generated. make: *** [CATS.o] Error 1 ERROR: compilation failed for package ‘CATS’ 
* removing ‘/Library/Frameworks/R.framework/Versions/3.4/Resources/library/CATS’ adris-imac:Desktop gwallace$ 

回答

1

:,使用R CMD INSTALL CATS_1.02.tar.gz命令行安裝時出現錯誤:

install.packages("CATS_1.02.tar.gz") 
Warning in install.packages : package ‘CATS_1.02.tar.gz’ is not available (for R version 3.4.1) 
library(CATS) Error in library(CATS) : there is no package called ‘CATS’ 

更新它會出現代碼打破:

res<-.Call("cats", 
      as.double(freq),as.double(freq2),as.integer(ncases),as.integer(ncontrols), 
      as.integer(ncases2),as.integer(ncontrols2),as.double(risk),as.double(risk2), 
      as.double(pisamples),as.double(prevalence),as.double(prevalence2), 
      as.integer(additive),as.integer(recessive),as.integer(dominant), 
      as.integer(multiplicative),as.double(alpha),as.double(pimarkers)) 

.Call用於調用ex ternal C/C++代碼:

https://stat.ethz.ch/R-manual/R-devel/library/base/html/CallExternal.html

沒有這個代碼將R腳本將無法工作。

而且我測試安裝包,它似乎精細安裝:

> install.packages("CATS_1.02.tar.gz") 
> library(CATS) 
> R.version 
platform  x86_64-redhat-linux-gnu 
arch   x86_64 
os    linux-gnu 
system   x86_64, linux-gnu 
status 
major   3 
minor   4.1 
year   2017 
month   06 
day   30 
svn rev  72865 
language  R 
version.string R version 3.4.1 (2017-06-30) 
nickname  Single Candle 
> CATS::cats() 
$P.one.study 
[1] 0.961869 

$P.first.stage 
[1] 0.9806984 

$P.rep.study 
[1] 0.8297875 

$P.joint.min 
[1] 0.9999998 

$P.joint 
[1] 0.9529604 

$pi 
[1] 0.5 

$T.one.study 
[1] 5.326724 

$T.first.stage 
[1] 2.951729 

$T.second.stage.rep 
[1] 4.000192 

$T.second.stage.joint 
[1] 5.30794 

$E.Disease.freq.cases1 
[1] 0.6 

$E.Disease.freq.controls1 
[1] 0.4888889 

$E.Disease.freq.cases2 
[1] 0.6 

$E.Disease.freq.controls2 
[1] 0.4888889 

$options 
    freq freq2 ncases ncontrols ncases2 ncontrols2 risk risk2 pisamples prevalence prevalence2 additive recessive dominant multiplicative alpha pimarkers 
[1,] 0.5 -1 500  500  500  500 1.5 -1  -1  0.1   -1  0   0  0    1 1e-07 0.00316 

attr(,"class") 
[1] "CATS" 

UPDATE

根據最新鐺錯誤也許嘗試:

R CMD INSTALL --configure-args="CFLAGS=-Wno-return-type CXXFLAGS=-Wno-return-type" CATS_1.02.tar.gz 

更新2

另請嘗試將以下內容添加到~/.R/Makevars

CFLAGS=-Wno-return-type 
CXXFLAGS=-Wno-return-type 

然後重新安裝包:

R CMD INSTALL --clean --preclean CATS_1.02.tar.gz 
+0

文斯,感謝您抽出時間來看看這個。我想知道是什麼阻礙了我安裝這個軟件包。我的平臺上安裝了Single Candle版本:x86_64-apple-darwin15.6.0 – Adri

+0

安裝時是否有錯誤?如果是這樣,請用錯誤日誌更新問題。 – Vince

+0

已更新。我重新安裝了R和Rstudio,但我得到了相同的警告和錯誤。 – Adri