Hello,
I am attempting to do an association analysis with GWASTools. I have followed the example script for making the GenotypeData object, and when I call variables with get__() I see exactly what I would expect. I am using imputed genotypes. The genotype file that I am working with is only a small subset of chromosome 1 (I wanted to troubleshoot with a small dataset before scaling up), and the phenotype includes only males, although there are females in the file. Females just have NA for all phenotypes except sex. When I call getSex(), I see a vector of Ms and Fs as I would expect.
My R commands are listed below:
ncfile <- tempfile()
snpfile <- tempfile()
scanfile <- tempfile()
probfile <- file.path("/path1/dcct_1kg_chr1_1_0")
sampfile <- file.path("/path2/EDcph.sample")
ncdfImputedDosage(input.files=c(probfile,sampfile), ncdf.filename=ncfile,chromosome=1, input.type="IMPUTE2",input.dosage=FALSE,
snp.annot.filename=snpfile, scan.annot.filename=scanfile, verbose=TRUE)
nc <- NcdfGenotypeReader(ncfile)
scanAnnot <- getobj(scanfile)
snpAnnot <- getobj(snpfile)
genoData <- GenotypeData(nc,scanAnnot=scanAnnot, snpAnnot=snpAnnot)
# Two different association tests give me different errors, which makes me think something is wrong with the GenotypeData object, but I can't figure out what:
> assocTestCPH(genoData, event="pheno", time.to.event="time_to_event", covars=c("Age"),
+ chromosome.set = 1, block.size=500, verbose=TRUE, outfile="try_cph.out")
Error in `[<-`(`*tmp*`, sample.dat$sex == "F", YchromCode(genoData), value = FALSE) :
subscript out of bounds
> mydat <- assocTestRegression(genoData, outcome="pheno", model.type="logistic", dosage=TRUE)
Determining chromosome and SNP information...
Gene action check...
Covariate check...
Interaction Variable check...
Models checks...
Determining number of SNP blocks for Chromosome 1 ...
Determining which Scans to keep for Chromosome 1 ...
Beginning calculations for Chromosome 1 ...
Error in R_nc_inq_varndims: NetCDF: Not a valid ID
Error in varndims.ncdf(nc, varid) : error returned from C call
Here is the sessionInfo:
R version 3.1.1 (2014-07-10)
Platform: x86_64-apple-darwin13.1.0 (64-bit)
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] parallel stats graphics grDevices utils datasets methods base
other attached packages:
[1] GWASTools_1.10.1 sandwich_2.3-2 gdsfmt_1.1.0.1 ncdf_1.6.8 Biobase_2.24.0 BiocGenerics_0.10.0
loaded via a namespace (and not attached):
[1] DBI_0.3.1 DNAcopy_1.38.1 grid_3.1.1 GWASExactHW_1.01 lattice_0.20-29 lmtest_0.9-33 quantsmooth_1.30.0 RSQLite_0.11.4 splines_3.1.1 survival_2.37-7
[11] tools_3.1.1 zoo_1.7-11
I would greatly appreciate any suggestions.
-Melody Palmer, University of Washington Dept of Urology, mryner@uw.edu
