Error in plinkCoxSurv - could not find function "coxph.fit"

[zhaiting]

Hi, I am just getting started and tested run the example for plinkCoxSurv. But the output said below, and generated blank documents ("impute_example.coxph" and "impute_example.snps_removed"). I've tried uninstall and install the "survival" package but seems not working. Any suggestion and help will be appreciated! Thanks in advance!

Analyzing part 1/1... Error in checkForRemoteErrors(val) : 3 nodes produced errors; first error: could not find function "coxph.fit"

[MichelledeGroot]

Hi, I ran into the same issue, have you found a solution for this?

[aarizvi]

What version of gwasuvrvivr are you using? Are you installing it from Bioconductor or GitHub?

[MichelledeGroot]

What version of gwasuvrvivr are you using? Are you installing it from Bioconductor or GitHub?

I'm using gwasurvivr_1.1. Installed from GitHub.

[aarizvi]

Can you send your sesssionInfo()?

[aarizvi]

I just tried the function and it works on my end .... Can you revert to the Bioconductor version?

# remove.packages('gwasurvivr') - if you need to remove the package
# install.packages("BiocManager") - if you need to install BiocManager
BiocManager::install("gwasurvivr")

> packageVersion("gwasurvivr")
[1] ‘1.14.0’

library(gwasurvivr)
bed.file <- system.file(package="gwasurvivr",
                        "extdata",
                        "plink_example.bed")
covariate.file <- system.file(package="gwasurvivr", 
                              "extdata",
                              "simulated_pheno.txt")
covariate.file <- read.table(covariate.file,
                             sep=" ",
                             header=TRUE,
                             stringsAsFactors = FALSE)
covariate.file$SexFemale <- ifelse(covariate.file$sex=="female", 1L, 0L)
sample.ids <- covariate.file[covariate.file$group=="experimental",]$ID_2
plinkCoxSurv(bed.file=bed.file,
             covariate.file=covariate.file,
             id.column="ID_2",
             sample.ids=sample.ids,
             time.to.event="time",
             event="event",
             covariates=c("age", "SexFemale", "DrugTxYes"),
             inter.term=NULL,
             print.covs="only",
             out.file="impute_example",
             chunk.size=50,
             maf.filter=0.005,
             flip.dosage=TRUE,
             verbose=TRUE,
             clusterObj=NULL)  

Output:

Covariates included in the models are: age, DrugTxYes, SexFemale
52 samples are included in the analysis
Start file conversion from PLINK BED to SNP GDS ...
    BED file: '/Library/Frameworks/R.framework/Versions/4.2/Resources/library/gwasurvivr/extdata/plink_example.bed'
        SNP-major mode (Sample X SNP), 78 bytes
    FAM file: '/Library/Frameworks/R.framework/Versions/4.2/Resources/library/gwasurvivr/extdata/plink_example.fam'
    BIM file: '/Library/Frameworks/R.framework/Versions/4.2/Resources/library/gwasurvivr/extdata/plink_example.bim'
Fri Dec  9 08:25:31 2022     (store sample id, snp id, position, and chromosome)
    start writing: 100 samples, 3 SNPs ...
[==================================================] 100%, completed, 0s  
Fri Dec  9 08:25:31 2022    Done.
Optimize the access efficiency ...
Clean up the fragments of GDS file:
    open the file '/var/folders/m5/y4ws0kg56kvcd36_p_gkjxzh0000gn/T//RtmpA4qf8X/fea14aeffa79.gds' (3.3K)
    # of fragments: 42
    save to '/var/folders/m5/y4ws0kg56kvcd36_p_gkjxzh0000gn/T//RtmpA4qf8X/fea14aeffa79.gds.tmp'
    rename '/var/folders/m5/y4ws0kg56kvcd36_p_gkjxzh0000gn/T//RtmpA4qf8X/fea14aeffa79.gds.tmp' (3.0K, reduced: 264B)
    # of fragments: 20
***** Compression time ******
User:0.034
System: 0.025
Elapsed: 0.066
*****************************
Analyzing part 1/1...
0 SNPs were removed from the analysis for not meeting
the given threshold criteria or for having MAF = 0
List of removed SNPs are saved to impute_example.snps_removed
In total 3 SNPs were included in the analysis
The Cox model results output was saved to impute_example.coxph
Analysis completed on 2022-12-09 at 08:25:45

[MichelledeGroot]

Unfortunately this does not help. I have tried R-3.6.3 with gwasurvivr 1.2.0 and R-4.2.0 with gwasurvivr 1.16.0 and the example still does not work.

Could you perhaps share your session info so I can compare some more versions?

[aarizvi]

Sure

> sessionInfo()
R version 4.2.0 (2022-04-22)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Catalina 10.15.7

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] VariantAnnotation_1.42.1    Rsamtools_2.12.0           
 [3] Biostrings_2.64.1           XVector_0.36.0             
 [5] SummarizedExperiment_1.26.1 Biobase_2.56.0             
 [7] GenomicRanges_1.48.0        GenomeInfoDb_1.32.4        
 [9] IRanges_2.30.1              S4Vectors_0.34.0           
[11] MatrixGenerics_1.8.1        matrixStats_0.62.0         
[13] BiocGenerics_0.42.0        

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.8.3             lattice_0.20-45          prettyunits_1.1.1       
 [4] png_0.1-7                assertthat_0.2.1         digest_0.6.29           
 [7] utf8_1.2.2               BiocFileCache_2.4.0      R6_2.5.1                
[10] RSQLite_2.2.18           httr_1.4.3               pillar_1.7.0            
[13] zlibbioc_1.42.0          rlang_1.0.6              GenomicFeatures_1.48.4  
[16] progress_1.2.2           curl_4.3.2               rstudioapi_0.13         
[19] blob_1.2.3               Matrix_1.5-3             BiocParallel_1.30.4     
[22] stringr_1.4.0            RCurl_1.98-1.8           bit_4.0.4               
[25] biomaRt_2.52.0           DelayedArray_0.22.0      rtracklayer_1.56.1      
[28] compiler_4.2.0           pkgconfig_2.0.3          tidyselect_1.1.2        
[31] KEGGREST_1.36.3          tibble_3.1.7             GenomeInfoDbData_1.2.8  
[34] codetools_0.2-18         XML_3.99-0.11            fansi_1.0.3             
[37] crayon_1.5.1             dplyr_1.0.9              dbplyr_2.2.0            
[40] GenomicAlignments_1.32.1 bitops_1.0-7             rappdirs_0.3.3          
[43] grid_4.2.0               lifecycle_1.0.1          DBI_1.1.2               
[46] magrittr_2.0.3           cli_3.4.1                stringi_1.7.6           
[49] cachem_1.0.6             xml2_1.3.3               ellipsis_0.3.2          
[52] filelock_1.0.2           vctrs_0.4.1              generics_0.1.2          
[55] rjson_0.2.21             restfulr_0.0.15          tools_4.2.0             
[58] bit64_4.0.5              BSgenome_1.64.0          glue_1.6.2              
[61] purrr_0.3.4              hms_1.1.1                yaml_2.3.5              
[64] parallel_4.2.0           fastmap_1.1.0            AnnotationDbi_1.58.0    
[67] memoise_2.0.1            BiocIO_1.6.0       

I"m wondering why you are getting cluster errors, since this current code snippet passes a NULL cluster object

You could always try debugging using debugonce('plinkCoxSurv') and seeing where it breaks

[MichelledeGroot]

Thanks!

It seems to break on:

debug: cox.out <- t(parallel::parApply(cl = cl, 1L, X = genotypes, FUN = survFit, 
    cox.params = cox.params, print.covs = print.covs))
Browse[2]> 
Error in checkForRemoteErrors(val) : 
  2 nodes produced errors; first error: could not find function "coxph.fit"

where cl = socket cluster with 2 nodes on host ‘localhost’

[aarizvi]

Hm. Its something to do with the cluster object. For some reason the reason the package isnt being exported to your cluster objects properly.

Can you assign survFit locally?

survFit <- function(SNP, cox.params, print.covs){

    if(is.null(cox.params$pheno.file)){
        SNP <- SNP[!is.na(SNP)]
        X <- matrix(SNP, ncol = 1)
        Y <- cox.params$Y[!is.na(SNP)]
        ROWNAMES <- cox.params$ROWNAMES[!is.na(SNP)]
    }else{
        ## creating model matrix
        X <- cbind(SNP, cox.params$pheno.file)
        X <- X[!is.na(SNP),]
        Y <- cox.params$Y[!is.na(SNP)]
        ROWNAMES <- cox.params$ROWNAMES[!is.na(SNP)]
    }

        ## run fit with pre-defined parameters including INIT
        fit <- coxph.fit(X,
                         Y,
                         cox.params$STRATA,
                         cox.params$OFFSET,
                         cox.params$INIT, 
                         cox.params$CONTROL,
                         cox.params$WEIGHTS,
                         cox.params$METHOD, 
                         ROWNAMES)

        ## extract statistics
        if(print.covs=="only") {
            coef <- fit$coefficients[1]
            serr <- sqrt(diag(fit$var)[1])
            n.sample <- nrow(X)
            n.event <- sum(!grepl("[+]", as.character(Y)))
            res <- cbind(coef=coef, serr=serr, n.sample=n.sample, n.event=n.event)
            return(res)

        } else if(print.covs=="all"){

            coef <- fit$coefficients
            serr <- sqrt(diag(fit$var))
            res <- cbind(coef, serr)
            res.names <- dimnames(res)
            n.sample <- nrow(X)
            n.event <- sum(!grepl("[+]", as.character(Y)))
            res <- c(res, n.sample, n.event)
            names(res) <- c(paste(toupper(res.names[[2]][1]),
                                  res.names[[1]],
                                  sep="_"),
                            paste(toupper(res.names[[2]][2]),
                                  res.names[[1]],
                                  sep="_"),
                            "N", 
                            "N.EVENT"
                            )
            return(res)

        }
}

and change the line that assigns the fit object to from coxph.fit alone to survival::coxph.fit and see if it works ?

[MichelledeGroot]

I added survival:: to the function and added my lib paths at the beginning of the function. This seems to have resolved the coxph.fit issue.

It now fails at the same step:

Error in checkForRemoteErrors(val) : 
  2 nodes produced errors; first error: incorrect number of dimensions

[aarizvi]

Hm, a bit unclear. Not the most informative error, but usually when that comes, it means that some error occurred in the compute while running parallel.

I just added a PR to the devel branch.

devtools::install_github("suchestoncampbelllab/gwasurvivr", ref = "devel") 

Can you try with this version? Package version should 1.2

[MichelledeGroot]

Apperently my R_LIBS_USER was not set correctly. I changed it to my installation dir and now it works. Thanks for the help!

@ZzzPIPI Do you also have gwasurvivr installed in a different lib location? Maybe you have the same issue, in that case it should work after changing R_LIBS_USER to your lib dir.