Read from vcz

.Rbuildignore

@@ -1 +1,3 @@
 thirdParty/cget
+thirdParty/tensorstore-0.1.65
+src/build

.gitignore

@@ -1,11 +1,14 @@
+**/.DS_Store
 /archived
 thirdParty/cget
 thirdParty/zlib-1.2.8
 thirdParty/bgen/.lock-waf_linux2_build
 thirdParty/bgen/build
 thirdParty/bgen/.hg
+src/*.dylib
 src/*.o
 src/*.so
 src/backup/*
+src/build
 extdata/input/nfam_100_nindep_0_step1_includeMoreRareVariants_poly.raw
 *.swp

R/Geno.R

@@ -24,6 +24,7 @@ checkGenoInput = function(bgenFile = "",
                  vcfFile = "",
 		 vcfFileIndex = "",
                  vcfField = "DS",
+                 vczFile = "",
                  savFile = "",
 		 savFileIndex = "",
                  sampleFile = "",
@@ -60,6 +61,9 @@ checkGenoInput = function(bgenFile = "",
         }
         dosageFileType = "vcf"
 
+    }else if (vczFile != "") {
+      Check_File_Exist(vczFile)
+      dosageFileType <- "vcz"
     }else if (savFile != "") {
         Check_File_Exist(savFile, "savFile")
 	#Check_File_Exist(savFileIndex, "savFileIndex")
@@ -124,6 +128,7 @@ setGenoInput = function(bgenFile = "",
                  vcfFile = "",
                  vcfFileIndex = "",
                  vcfField = "DS",
+                 vczFile = "",
                  savFile = "",
                  savFileIndex = "",
                  sampleFile = "",
@@ -142,6 +147,7 @@ setGenoInput = function(bgenFile = "",
                  vcfFile = vcfFile,
                  vcfFileIndex = vcfFileIndex,
                  vcfField = vcfField,
+                 vczFile = vczFile,
                  savFile = savFile,
                  savFileIndex = savFileIndex,
                  bedFile = bedFile,
@@ -437,6 +443,11 @@ if(FALSE){
  #   markerInfo[,POS:=NULL]
  # }
 
+  if (dosageFileType == "vcz") {
+    setVCZobjInCPP(vczFile, t_SampleInModel = sampleInModel)
+    markerInfo = NULL
+  }
+
   genoList = list(markerInfo = markerInfo, genoType = dosageFileType, anyInclude = anyInclude)
   return(genoList)
 }

R/RcppExports.R

@@ -81,6 +81,10 @@ setVCFobjInCPP <- function(t_vcfFileName, t_vcfFileIndex, t_vcfField, t_SampleIn
     invisible(.Call('_SAIGE_setVCFobjInCPP', PACKAGE = 'SAIGE', t_vcfFileName, t_vcfFileIndex, t_vcfField, t_SampleInModel))
 }
 
+setVCZobjInCPP <- function(t_vczFileName, t_SampleInModel) {
+    invisible(.Call('_SAIGE_setVCZobjInCPP', PACKAGE = 'SAIGE', t_vczFileName, t_SampleInModel))
+}
+
 setSAIGEobjInCPP <- function(t_XVX, t_XXVX_inv, t_XV, t_XVX_inv_XV, t_Sigma_iXXSigma_iX, t_X, t_S_a, t_res, t_mu2, t_mu, t_varRatio_sparse, t_varRatio_null, t_cateVarRatioMinMACVecExclude, t_cateVarRatioMaxMACVecInclude, t_SPA_Cutoff, t_tauvec, t_traitType, t_y, t_impute_method, t_flagSparseGRM, t_isFastTest, t_pval_cutoff_for_fastTest, t_locationMat, t_valueVec, t_dimNum, t_isCondition, t_condition_genoIndex, t_is_Firth_beta, t_pCutoffforFirth, t_offset, t_resout) {
     invisible(.Call('_SAIGE_setSAIGEobjInCPP', PACKAGE = 'SAIGE', t_XVX, t_XXVX_inv, t_XV, t_XVX_inv_XV, t_Sigma_iXXSigma_iX, t_X, t_S_a, t_res, t_mu2, t_mu, t_varRatio_sparse, t_varRatio_null, t_cateVarRatioMinMACVecExclude, t_cateVarRatioMaxMACVecInclude, t_SPA_Cutoff, t_tauvec, t_traitType, t_y, t_impute_method, t_flagSparseGRM, t_isFastTest, t_pval_cutoff_for_fastTest, t_locationMat, t_valueVec, t_dimNum, t_isCondition, t_condition_genoIndex, t_is_Firth_beta, t_pCutoffforFirth, t_offset, t_resout))
 }
@@ -117,6 +121,18 @@ move_forward_iterator_Vcf <- function(i) {
     invisible(.Call('_SAIGE_move_forward_iterator_Vcf', PACKAGE = 'SAIGE', i))
 }
 
+set_iterator_inVcz <- function(chrom, beg_pd, end_pd) {
+    invisible(.Call('_SAIGE_set_iterator_inVcz', PACKAGE = 'SAIGE', chrom, beg_pd, end_pd))
+}
+
+check_Vcz_end <- function() {
+    .Call('_SAIGE_check_Vcz_end', PACKAGE = 'SAIGE')
+}
+
+move_forward_iterator_Vcz <- function(i) {
+    invisible(.Call('_SAIGE_move_forward_iterator_Vcz', PACKAGE = 'SAIGE', i))
+}
+
 fast_logistf_fit <- function(x, y, weight, offset, firth, col_fit, init, maxit, maxstep, maxhs, lconv, gconv, xconv, isfirthconverge) {
     .Call('_SAIGE_fast_logistf_fit', PACKAGE = 'SAIGE', x, y, weight, offset, firth, col_fit, init, maxit, maxstep, maxhs, lconv, gconv, xconv, isfirthconverge)
 }

R/SAIGE_SPATest_Marker.R

@@ -36,7 +36,54 @@ SAIGE.Marker = function(traitType,
   }
 
   ## set up an object for genotype
-  if(genoType != "vcf"){
+  if(genoType == "vcf"){
+   if(!isAnyInclude){	  
+    if(chrom == ""){
+      stop("chrom needs to be specified for single-variant assoc tests when using VCF as input\n")
+    }else{
+      set_iterator_inVcf("", chrom, 1, 250000000)
+    }
+   }
+    if(outIndex > 1){
+	move_forward_iterator_Vcf(outIndex*nMarkersEachChunk)    
+    }
+    isVcfEnd =  check_Vcf_end()
+    if(!isVcfEnd){
+    	#outIndex = 1
+    	genoIndex = rep("0", nMarkersEachChunk) 
+    	genoIndex_prev = rep("0", nMarkersEachChunk) 
+	#nChunks = outIndex + 1
+	is_marker_test = TRUE
+        i = outIndex
+    }else{
+	is_marker_test = FALSE    
+	stop("No markers are left in VCF")
+    }
+  } else if (genoType == "vcz") {
+    if (!isAnyInclude) {
+      if (chrom == "") {
+        stop("chrom needs to be specified for single-variant assoc tests when using VCZ as input\n")
+      } else {
+        set_iterator_inVcz(chrom, 1, 250000000)
+      }
+    }
+
+    if (outIndex > 1) {
+      move_forward_iterator_Vcz(outIndex * nMarkersEachChunk)
+    }
+
+    isVczEnd <- check_Vcz_end()
+
+    if (!isVczEnd) {
+      genoIndex <- rep("0", nMarkersEachChunk)
+      genoIndex_prev <- rep("0", nMarkersEachChunk)
+      is_marker_test <- TRUE
+      i <- outIndex
+    } else {
+      is_marker_test <- FALSE
+	    stop("No markers are left in VCZ")
+    }
+  }else{
       #markerInfo = objGeno$markerInfo
     if(LOCO){
       genoIndex = genoIndex[which(CHROM == chrom)]
@@ -71,30 +118,6 @@ SAIGE.Marker = function(traitType,
       is_marker_test = TRUE
       i = outIndex    
     }
-
-  }else{
-   if(!isAnyInclude){	  
-    if(chrom == ""){
-      stop("chrom needs to be specified for single-variant assoc tests when using VCF as input\n")
-    }else{
-      set_iterator_inVcf("", chrom, 1, 250000000)
-    }
-   }
-    if(outIndex > 1){
-	move_forward_iterator_Vcf(outIndex*nMarkersEachChunk)    
-    }
-    isVcfEnd =  check_Vcf_end()
-    if(!isVcfEnd){
-    	#outIndex = 1
-    	genoIndex = rep("0", nMarkersEachChunk) 
-    	genoIndex_prev = rep("0", nMarkersEachChunk) 
-	#nChunks = outIndex + 1
-	is_marker_test = TRUE
-        i = outIndex
-    }else{
-	is_marker_test = FALSE    
-	stop("No markers are left in VCF")
-    }
   }
 
   chrom = "InitialChunk"
@@ -105,7 +128,7 @@ SAIGE.Marker = function(traitType,
 #time_left = system.time({
 
 
-    if(genoType != "vcf"){	
+    if(genoType != "vcf" && genoType != "vcz"){	
       tempList = genoIndexList[[i]]
       genoIndex = as.character(format(tempList$genoIndex, scientific = FALSE))
       tempChrom = tempList$chrom
@@ -133,7 +156,7 @@ SAIGE.Marker = function(traitType,
     #print(ptm)
     #print("gc()")
     #print(gc())
-    if(genoType != "vcf"){
+    if(genoType != "vcf" && genoType != "vcz"){
       cat(paste0("(",Sys.time(),") ---- Analyzing Chunk ", i, "/", nChunks, ": chrom ", chrom," ---- \n"))
     }else{
       cat(paste0("(",Sys.time(),") ---- Analyzing Chunk ", i, " :  chrom ", chrom," ---- \n"))
@@ -150,6 +173,8 @@ SAIGE.Marker = function(traitType,
 
   if(genoType == "vcf"){
     isEnd_Output =  check_Vcf_end()
+  }else if (genoType == "vcz") {
+    isEnd_Output = check_Vcz_end()
   }else{
     isEnd_Output = (i==nChunks)	
   }
@@ -189,6 +214,12 @@ SAIGE.Marker = function(traitType,
     if(isVcfEnd){
 	is_marker_test = FALSE	     
     }
+  }else if(genoType == "vcz") {
+    isVczEnd = check_Vcz_end()
+    cat("isVczEnd ", isVczEnd, "\n")
+    if (isVczEnd) {
+      is_marker_test = FALSE
+    }
   }else{
     if(i > nChunks){
       is_marker_test = FALSE

R/SAIGE_Test_main.R

@@ -6,6 +6,7 @@
 #' @param vcfFile character. Path to vcf file
 #' @param vcfFileIndex character. Path to vcf index file. Indexed by tabix. Path to index for vcf file by tabix, .csi file using 'tabix --csi -p vcf file.vcf.gz'
 #' @param vcfField character. genotype field in vcf file to use. "DS" for dosages or "GT" for genotypes. By default, "DS".
+#' @param vczFile character. Path to VCZ group.
 #' @param savFile character. Path to sav file
 #' @param savFileIndex character. Path to index for sav file .s1r
 #' @param bedFile character. Path to bed file (PLINK)
@@ -64,6 +65,7 @@ SPAGMMATtest = function(bgenFile = "",
                  vcfFile = "",
                  vcfFileIndex = "",
                  vcfField = "DS",
+                 vczFile = "",
                  savFile = "",
                  savFileIndex = "",
                  bedFile="",
@@ -334,14 +336,15 @@ SPAGMMATtest = function(bgenFile = "",
     }
 
     nsample = length(obj.model$y)
-    cateVarRatioMaxMACVecInclude = c(cateVarRatioMaxMACVecInclude, nsample)	
-    
+    cateVarRatioMaxMACVecInclude = c(cateVarRatioMaxMACVecInclude, nsample)
+
     #in Geno.R
     objGeno = setGenoInput(bgenFile = bgenFile,
                  bgenFileIndex = bgenFileIndex,
                  vcfFile = vcfFile,   #not activate yet
                  vcfFileIndex = vcfFileIndex,
                  vcfField = vcfField,
+                 vczFile = vczFile,
                  savFile = savFile,
                  savFileIndex = savFileIndex,
                  sampleFile = sampleFile,

conda_env/createCondaEnvSAIGE_steps.txt

@@ -1,4 +1,4 @@
-conda create -n RSAIGE r-essentials r-base=4.1.2 python=3.10
+conda create -n RSAIGE r-essentials r-base=4.3.0 python=3.10
 conda activate RSAIGE
 conda install -c anaconda cmake
 conda install -c conda-forge gettext lapack r-matrix

extdata/input/genotype_100markers.vcz/.zgroup

@@ -0,0 +1,3 @@
+{
+    "zarr_format": 2
+}

extdata/input/genotype_100markers.vcz/call_genotype/.zarray

@@ -0,0 +1,25 @@
+{
+    "chunks": [
+        10,
+        1000,
+        2
+    ],
+    "compressor": {
+        "blocksize": 0,
+        "clevel": 7,
+        "cname": "zstd",
+        "id": "blosc",
+        "shuffle": 2
+    },
+    "dimension_separator": "/",
+    "dtype": "|i1",
+    "fill_value": null,
+    "filters": null,
+    "order": "C",
+    "shape": [
+        100,
+        10000,
+        2
+    ],
+    "zarr_format": 2
+}

extdata/input/genotype_100markers.vcz/call_genotype/.zattrs

@@ -0,0 +1,8 @@
+{
+    "_ARRAY_DIMENSIONS": [
+        "variants",
+        "samples",
+        "ploidy"
+    ],
+    "description": ""
+}

extdata/input/genotype_100markers.vcz/call_genotype_mask/.zarray

@@ -0,0 +1,25 @@
+{
+    "chunks": [
+        10,
+        1000,
+        2
+    ],
+    "compressor": {
+        "blocksize": 0,
+        "clevel": 7,
+        "cname": "zstd",
+        "id": "blosc",
+        "shuffle": 2
+    },
+    "dimension_separator": "/",
+    "dtype": "|b1",
+    "fill_value": null,
+    "filters": null,
+    "order": "C",
+    "shape": [
+        100,
+        10000,
+        2
+    ],
+    "zarr_format": 2
+}

extdata/input/genotype_100markers.vcz/call_genotype_mask/.zattrs

@@ -0,0 +1,8 @@
+{
+    "_ARRAY_DIMENSIONS": [
+        "variants",
+        "samples",
+        "ploidy"
+    ],
+    "description": ""
+}

extdata/input/genotype_100markers.vcz/call_genotype_phased/.zarray

@@ -0,0 +1,23 @@
+{
+    "chunks": [
+        10,
+        1000
+    ],
+    "compressor": {
+        "blocksize": 0,
+        "clevel": 7,
+        "cname": "zstd",
+        "id": "blosc",
+        "shuffle": 2
+    },
+    "dimension_separator": "/",
+    "dtype": "|b1",
+    "fill_value": null,
+    "filters": null,
+    "order": "C",
+    "shape": [
+        100,
+        10000
+    ],
+    "zarr_format": 2
+}

extdata/input/genotype_100markers.vcz/call_genotype_phased/.zattrs

@@ -0,0 +1,7 @@
+{
+    "_ARRAY_DIMENSIONS": [
+        "variants",
+        "samples"
+    ],
+    "description": ""
+}