update for new DataFrames syntax

OpenMendel · Oct 11, 2019 · 67dfe0e · 67dfe0e
1 parent 62e47a9
commit 67dfe0e
Show file tree

Hide file tree

Showing 3 changed files with 26 additions and 25 deletions.
diff --git a/Project.toml b/Project.toml
@@ -17,6 +17,7 @@ StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d"
 [compat]
 OrdinalMultinomialModels = ">=0.3.0"
 SnpArrays = ">=0.1.0"
+DataFrames = ">=0.19.0"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

diff --git a/src/gwas.jl b/src/gwas.jl
@@ -148,7 +148,7 @@ function ordinalgwas(
     genomat = SnpArrays.SnpArray(bedfile, bedn)
     # extra columns in design matrix to be tested
     testdf = DataFrame(fittednullmodel.mf.data) # TODO: not type stable here
-    testdf[:snp] = zeros(size(fittednullmodel.mm, 1))
+    testdf[!, :snp] = zeros(size(fittednullmodel.mm, 1))
     #mfalt = ModelFrame(testformula, testdf)
     #mfalt.terms.intercept = false # drop intercept
     #Z = similar(ModelMatrix(mfalt).m)
@@ -182,7 +182,7 @@ function ordinalgwas(
                         if snponly
                             copyto!(ts.Z, @view(genomat[bedrowinds, j]), impute=true, model=snpmodel)
                         else # snp + other terms
-                            copyto!(testdf[:snp], @view(genomat[bedrowinds, j]), impute=true, model=snpmodel)
+                            copyto!(testdf[!, :snp], @view(genomat[bedrowinds, j]), impute=true, model=snpmodel)
                             #mfalt = ModelFrame(testformula, testdf)
                             #mfalt.terms.intercept = false # drop intercept
                             #ts.Z[:] = ModelMatrix(mfalt).m
@@ -226,7 +226,7 @@ function ordinalgwas(
                                 @view(genomat[bedrowinds, j]), 
                                 impute=true, model=snpmodel)
                         else # snp + other terms
-                            copyto!(testdf[:snp], @view(genomat[bedrowinds, j]), 
+                            copyto!(testdf[!, :snp], @view(genomat[bedrowinds, j]), 
                                 impute=true, model=snpmodel)
                             #mfalt = ModelFrame(testformula, testdf)
                             #mfalt.terms.intercept = false # drop intercept
@@ -421,9 +421,9 @@ function ordinalsnpsetgwas(
         snpsetFile = CSV.read(snpset, header = [:snpset_id, :snp_id], delim = " ")
         #make sure it matches bim file
         biminfo = CSV.read(bimfile, header = [:chr, :snp_id, :c3, :bp, :c5, :c6], delim = "\t")
-        snpsetFile[:snp_id] == biminfo[:snp_id] || throw(ArgumentError("snp order in snpset file
+        snpsetFile[!, :snp_id] == biminfo[!, :snp_id] || throw(ArgumentError("snp order in snpset file
         must match (in the same order) bimfile")) 
-        snpset_ids = unique(snpsetFile[:snpset_id])
+        snpset_ids = unique(snpsetFile[!, :snpset_id])
         nSets = length(snpset_ids)
         setlength = 0
     elseif isa(snpset, Integer)
@@ -521,7 +521,7 @@ function ordinalsnpsetgwas(
                 "snpsetid,nsnps,l2normeffect,pval")
             for j in eachindex(snpset_ids)
                 snpset_id = snpset_ids[j]
-                snpinds = findall(snpsetFile[:snpset_id] .== snpset_id)
+                snpinds = findall(snpsetFile[!, :snpset_id] .== snpset_id)
                 q = length(snpinds)
                 Z = zeros(size(fittednullmodel.mm, 1), q)
                 γ̂ = Vector{Float64}(undef, q)
@@ -737,7 +737,7 @@ function ordinalgwasGxE(
     genomat = SnpArrays.SnpArray(bedfile, bedn)
     cc = SnpArrays.counts(genomat, dims=1) 
     mafs = SnpArrays.maf(genomat)
-    envvar = DataFrame(fittednullmodel.mf.data)[e]
+    envvar = DataFrame(fittednullmodel.mf.data)[!, e]
     testvec = zeros(size(fittednullmodel.mm, 1), 1)
     # create SNP mask vector
     if snpinds == nothing

diff --git a/test/runtests.jl b/test/runtests.jl
@@ -9,7 +9,7 @@ const snpsetfile = datadir * "/hapmap_snpsetfile.txt"
     @time ordinalgwas(@formula(trait ~ sex), covfile, plkfile, test=:score)
     @test isfile("ordinalgwas.null.txt")
     @test isfile("ordinalgwas.pval.txt")
-    scorepvals = open(CSV.read, "ordinalgwas.pval.txt")[6][1:5]
+    scorepvals = open(CSV.read, "ordinalgwas.pval.txt")[!, 6][1:5]
     @test isapprox(scorepvals, [1.0, 4.56531284e-3, 3.10828383e-5, 1.21686724e-5, 8.20686005e-3], rtol=1e-4)
     rm("ordinalgwas.null.txt", force=true)
     rm("ordinalgwas.pval.txt", force=true)
@@ -19,7 +19,7 @@ end
     @time ordinalgwas(@formula(trait ~ sex), covfile, plkfile, test=:LRT)
     @test isfile("ordinalgwas.null.txt")
     @test isfile("ordinalgwas.pval.txt")
-    lrtpvals = open(CSV.read, "ordinalgwas.pval.txt")[7][1:5]
+    lrtpvals = open(CSV.read, "ordinalgwas.pval.txt")[!, 7][1:5]
     @test isapprox(lrtpvals, [1.0, 1.91858366e-3, 1.80505056e-5, 5.87338471e-6, 8.08102258e-3], rtol=1e-4)
     rm("ordinalgwas.null.txt", force=true)
     rm("ordinalgwas.pval.txt", force=true)
@@ -30,15 +30,15 @@ end
     ordinalgwas(@formula(trait ~ sex), covfile, plkfile, test=:score, snpmodel=DOMINANT_MODEL)
     @test isfile("ordinalgwas.null.txt")
     @test isfile("ordinalgwas.pval.txt")
-    scorepvals = open(CSV.read, "ordinalgwas.pval.txt")[6][1:5]
+    scorepvals = open(CSV.read, "ordinalgwas.pval.txt")[!, 6][1:5]
     @test isapprox(scorepvals, [1.0, 0.14295, 0.000471942, 0.00555348, 0.000652844], rtol=1e-4)
     rm("ordinalgwas.null.txt", force=true)
     rm("ordinalgwas.pval.txt", force=true)
     # recessive model 
     ordinalgwas(@formula(trait ~ sex), covfile, plkfile, test=:score, snpmodel=RECESSIVE_MODEL)
     @test isfile("ordinalgwas.null.txt")
     @test isfile("ordinalgwas.pval.txt")
-    scorepvals = open(CSV.read, "ordinalgwas.pval.txt")[6][1:5]
+    scorepvals = open(CSV.read, "ordinalgwas.pval.txt")[!, 6][1:5]
     @test isapprox(scorepvals, [1.0, 0.00673612, 0.000279908, 4.15322e-5, 0.167642], rtol=1e-4)
     rm("ordinalgwas.null.txt", force=true)
     rm("ordinalgwas.pval.txt", force=true)
@@ -48,7 +48,7 @@ end
     ordinalgwas(@formula(trait ~ sex), covfile, plkfile, link=ProbitLink(), pvalfile="opm.pval.txt")
     @test isfile("ordinalgwas.null.txt")
     @test isfile("opm.pval.txt")
-    scorepvals = open(CSV.read, "opm.pval.txt")[6][1:5]
+    scorepvals = open(CSV.read, "opm.pval.txt")[!, 6][1:5]
     @test isapprox(scorepvals, [1.0, 1.00769167e-2, 2.62725649e-5, 1.08974849e-5, 5.10288399e-3], rtol=1e-4)
     rm("ordinalgwas.null.txt", force=true)
     rm("opm.pval.txt", force=true)
@@ -59,7 +59,7 @@ end
     @test isfile("ordinalgwas.null.txt")
     @test isfile("first5snps.pval.txt")
     @test countlines("first5snps.pval.txt") == 6
-    scorepvals = open(CSV.read, "first5snps.pval.txt")[6]
+    scorepvals = open(CSV.read, "first5snps.pval.txt")[!, 6]
     @test isapprox(scorepvals, [1.0, 4.56531284e-3, 3.10828383e-5, 1.21686724e-5, 8.20686005e-3], rtol=1e-4)
     rm("ordinalgwas.null.txt", force=true)
     rm("first5snps.pval.txt", force=true)
@@ -70,7 +70,7 @@ end
     @time ordinalgwas(@formula(trait ~ sex), covfile, plkfile, test=:score, covrowinds=1:300, bedrowinds=1:300)
     @test isfile("ordinalgwas.null.txt")
     @test isfile("ordinalgwas.pval.txt")
-    scorepvals = open(CSV.read, "ordinalgwas.pval.txt")[6][1:5]
+    scorepvals = open(CSV.read, "ordinalgwas.pval.txt")[!, 6][1:5]
     @test isapprox(scorepvals, [1.0, 0.00355969, 0.000123604, 5.2213e-6, 0.00758234], rtol=1e-4)
     rm("ordinalgwas.null.txt", force=true)
     rm("ordinalgwas.pval.txt", force=true)
@@ -82,7 +82,7 @@ end
     testformula=@formula(trait ~ snp + snp & sex))
     @test isfile("ordinalgwas.null.txt")
     @test isfile("GxE.pval.txt")
-    scorepvals = open(CSV.read, "GxE.pval.txt")[6][1:5]
+    scorepvals = open(CSV.read, "GxE.pval.txt")[!, 6][1:5]
     @test isapprox(scorepvals, [1.0, 1.74460104e-2, 1.66707324e-4, 4.76376246e-5, 2.91384712e-2], rtol=1e-4)
     rm("ordinalgwas.null.txt", force=true)
     rm("GxE.pval.txt", force=true)
@@ -91,7 +91,7 @@ end
     testformula=@formula(trait ~ snp + snp & sex), test=:LRT, snpinds=1:5)
     @test isfile("ordinalgwas.null.txt")
     @test isfile("GxE.pval.txt")
-    lrtpvals = open(CSV.read, "GxE.pval.txt")[end]
+    lrtpvals = open(CSV.read, "GxE.pval.txt")[!, end]
     @test isapprox(lrtpvals, [1.0, 7.22410973e-3, 1.01730983e-4, 1.88174211e-5, 2.88295705e-2], rtol=1e-4)
     rm("ordinalgwas.null.txt", force=true)
     rm("GxE.pval.txt", force=true)
@@ -104,7 +104,7 @@ end
     snpset=250)
     @test isfile("ordinalgwas.null.txt")
     @test isfile("snpset.pval.txt")
-    scorepvals = open(CSV.read, "snpset.pval.txt")[end][1:5]
+    scorepvals = open(CSV.read, "snpset.pval.txt")[!, end][1:5]
     #@test isapprox(scorepvals, [1.0, 1.74460104e-2, 1.66707324e-4, 4.76376246e-5, 2.91384712e-2], rtol=1e-4)
     rm("ordinalgwas.null.txt", force=true)
     rm("snpset.pval.txt", force=true)
@@ -113,7 +113,7 @@ end
     snpset=25, test=:LRT)
     @test isfile("ordinalgwas.null.txt")
     @test isfile("snpset.pval.txt")
-    lrtpvals = open(CSV.read, "snpset.pval.txt")[end][1:5]
+    lrtpvals = open(CSV.read, "snpset.pval.txt")[!, end][1:5]
     @test isapprox(lrtpvals, [2.1817554071810948e-13, 0.2865769729670889, 0.32507802233937966,
     0.3344823237332578, 0.42948375949508427], rtol=1e-4)
     rm("ordinalgwas.null.txt", force=true)
@@ -125,7 +125,7 @@ end
     snpset = snpsetfile)
     @test isfile("ordinalgwas.null.txt")
     @test isfile("snpset.pval.txt")
-    scorepvals = open(CSV.read, "snpset.pval.txt")[end][1:5]
+    scorepvals = open(CSV.read, "snpset.pval.txt")[!, end][1:5]
     @test isapprox(scorepvals, [1.72134e-5, 0.036925, 0.747855,
      0.0276508, 0.611958], rtol=1e-4)
     rm("ordinalgwas.null.txt", force=true)
@@ -135,7 +135,7 @@ end
     snpset = snpsetfile, test = :lrt)
     @test isfile("ordinalgwas.null.txt")
     @test isfile("snpset.pval.txt")
-    lrtpvals = open(CSV.read, "snpset.pval.txt")[end][1:5]
+    lrtpvals = open(CSV.read, "snpset.pval.txt")[!, end][1:5]
     @test isapprox(lrtpvals, [6.75377e-13, 0.000256566, 0.359382,
      0.000163268, 0.0867508], rtol=1e-4)
     rm("ordinalgwas.null.txt", force=true)
@@ -168,14 +168,14 @@ end
     ordinalgwasGxE(@formula(trait ~ sex), covfile, plkfile, :sex, pvalfile = "gxe_snp.pval.txt",
     snpinds=1:5, test=:score)
     @test isfile("gxe_snp.pval.txt")
-    scorepvals = open(CSV.read, "gxe_snp.pval.txt")[end][1:5]
+    scorepvals = open(CSV.read, "gxe_snp.pval.txt")[!, end][1:5]
     @test isapprox(scorepvals, [1.0, 0.637742242597749, 0.9667114198051628,
     0.26352674694121003, 0.7811133315582837], rtol=1e-4)
     rm("gxe_snp.pval.txt", force=true)
     ordinalgwasGxE(@formula(trait ~ sex), covfile, plkfile, "sex", pvalfile = "gxe_snp.pval.txt",
     snpinds=1:5, test=:LRT)
     @test isfile("gxe_snp.pval.txt")
-    lrtpvals = open(CSV.read, "gxe_snp.pval.txt")[end][1:5]
+    lrtpvals = open(CSV.read, "gxe_snp.pval.txt")[!, end][1:5]
     @test isapprox(lrtpvals, [1.0, 0.6279730133445315, 0.9671662821946985,
     0.26693502209463904, 0.7810214899265426], rtol=1e-4)
     rm("gxe_snp.pval.txt", force=true)
@@ -202,7 +202,7 @@ end
         ordinalgwas(@formula(trait ~ sex), covfile, plinkfile, pvalfile = pvalfile)
         @test isfile(pvalfile)
         if chr == 1
-            pvals_chr1 = open(CSV.read, pvalfile)[6][1:5]
+            pvals_chr1 = open(CSV.read, pvalfile)[!, 6][1:5]
             @test isapprox(pvals_chr1, [1.0, 4.56531284e-3, 3.10828383e-5, 1.21686724e-5, 8.20686005e-3], rtol=1e-4)    
         end
         rm(plinkfile * ".pval.txt", force=true)
@@ -214,7 +214,7 @@ end
         ordinalgwas(nm, plinkfile, pvalfile = pvalfile)
         @test isfile(pvalfile)
         if chr == 1
-            pvals_chr1 = open(CSV.read, pvalfile)[6][1:5]
+            pvals_chr1 = open(CSV.read, pvalfile)[!, 6][1:5]
             @test isapprox(pvals_chr1, [1.0, 4.56531284e-3, 3.10828383e-5, 1.21686724e-5, 8.20686005e-3], rtol=1e-4)    
         end
         rm(pvalfile, force=true)
@@ -227,7 +227,7 @@ end
         ordinalgwas(nm, bedfile, bimfile, 324; pvalfile = pvalfile)
         @test isfile(pvalfile)
         if chr == 1
-            pvals_chr1 = open(CSV.read, pvalfile)[6][1:5]
+            pvals_chr1 = open(CSV.read, pvalfile)[!, 6][1:5]
             @test isapprox(pvals_chr1, [1.0, 4.56531284e-3, 3.10828383e-5, 1.21686724e-5, 8.20686005e-3], rtol=1e-4)    
         end
         rm(pvalfile, force=true)