Remove_asint (#183)

* remove genotype_asint

* add Arrow writing instead of CSV

Author: olivierlabayle

Manifest.toml

@@ -5,6 +5,7 @@ version = "0.1.0"
 
 [deps]
 ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
+Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45"
 BGEN = "6db4b851-9beb-4b83-9d64-eb1cfb37721d"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 CategoricalArrays = "324d7699-5711-5eae-9e2f-1d82baa6b597"
@@ -23,5 +24,5 @@ CategoricalArrays = "0.10"
 Combinatorics = "1.0"
 DataFrames = "1.2"
 SnpArrays = "0.3"
-YAML = "0.4"
 TMLE = "0.11"
+YAML = "0.4"

Project.toml

@@ -82,10 +82,6 @@ function parse_commandline()
             required = false
             help = "Interaction orders to be estimated"
             default = "1,2"
-
-        "--genotypes-as-int"
-            action = :store_true
-            help = "If true, genotypes are encoded as the number of minor alleles."
     end
 
     @add_arg_table s["from-param-file"] begin

bin/tmle_inputs.jl

@@ -8,6 +8,7 @@ using Mmap
 using YAML
 using Combinatorics
 using TMLE
+using Arrow
 
 
 ###############################################################################

src/TargeneCore.jl

@@ -194,7 +194,6 @@ function tmle_inputs_from_actors(parsed_args)
     orders = TargeneCore.parse_orders(parsed_args["from-actors"]["orders"])
     extraW = TargeneCore.read_txt_file(parsed_args["from-actors"]["extra-confounders"])
     extraC = TargeneCore.read_txt_file(parsed_args["from-actors"]["extra-covariates"])
-    genotypes_asint = parsed_args["from-actors"]["genotypes-as-int"]
 
     # Retrieve SNPs and environmental treatments
     bqtls, transactors, extraT = TargeneCore.treatments_from_actors(
@@ -204,7 +203,7 @@ function tmle_inputs_from_actors(parsed_args)
     )
     # Genotypes and final dataset
     variants = TargeneCore.all_variants(bqtls, transactors)
-    genotypes = TargeneCore.call_genotypes(bgen_prefix, variants, call_threshold; asint=genotypes_asint)
+    genotypes = TargeneCore.call_genotypes(bgen_prefix, variants, call_threshold)
     data = TargeneCore.merge(traits, pcs, genotypes)
 
     # Parameter files
@@ -215,5 +214,5 @@ function tmle_inputs_from_actors(parsed_args)
     )
 
     # write data
-    CSV.write(string(outprefix, ".data.csv"), data)
+    Arrow.write(string(outprefix, ".data.arrow"), data)
 end

src/tmle_inputs/from_actors.jl

@@ -26,33 +26,6 @@ end
 check_genotypes_encoding(val::T, type) where T = 
     T <: type || throw(MismatchedCaseControlEncodingError())
 
-get_genotype_encoding(val::NamedTuple) = typeof(val.case) <: Real ? Real : String
-get_genotype_encoding(val::T) where {T<:String} = String
-get_genotype_encoding(val::T) where {T<:Real} = Real
-get_genotype_encoding(val) = throw(ArgumentError(string("Genotype value(s):", val, " not allowed")))
-
-"""
-We enforce that all genotypes are encoded in the same way.
-That is either with an integer representing the number of minor alleles 
-or the string representation.
-"""
-function get_genotype_encoding(parameters, variants::Set{Symbol})
-    genotypes_encoding = nothing
-    for Ψ in parameters
-        for (T, val) in zip(keys(Ψ.treatment), Ψ.treatment)
-            if T in variants
-                if genotypes_encoding === nothing 
-                    genotypes_encoding = get_genotype_encoding(val)
-                else
-                    check_genotypes_encoding(val, genotypes_encoding)
-                end
-            end
-        end
-    end
-    return genotypes_encoding <: Real
-end
-
-
 
 function get_variables(parameters, traits, pcs)
     variants = Set{Symbol}()
@@ -75,9 +48,6 @@ end
 fix_mismatch(variant, allele::T, actual_alleles) where T = 
     throw(ArgumentError(string("Can't deal with ", variant, "'s allele ", allele, " of type: ", T)))
 
-fix_mismatch(variant, allele::Real, actual_alleles) = 
-    allele ∈ actual_alleles || throw(AbsentAlleleError(variant, allele))
-
 function fix_mismatch(variant, allele::String, actual_alleles)
     length(allele) == 2 || throw(NotSNPAndWontFixError(variant, allele))
     rev_allele = reverse(allele)
@@ -203,12 +173,10 @@ function tmle_inputs_from_param_files(parsed_args)
 
     # Genotypes and full data
     variables = TargeneCore.get_variables(parameters, traits, pcs)
-    genotypes_asint = TargeneCore.get_genotype_encoding(parameters, variables.variants)
     genotypes = TargeneCore.call_genotypes(
         bgen_prefix, 
         Set(string.(variables.variants)), 
-        call_threshold; 
-        asint=genotypes_asint
+        call_threshold
     )
     data = TargeneCore.merge(traits, pcs, genotypes)
 

src/tmle_inputs/from_param_files.jl

@@ -12,7 +12,7 @@ read_data(filepath) = CSV.read(filepath, DataFrame, types=Dict(:SAMPLE_ID => Str
 
 function write_tmle_inputs(outprefix, final_dataset, parameters; batch_size=nothing)
     # Write final_dataset
-    CSV.write(string(outprefix, ".data.csv"), final_dataset)
+    Arrow.write(string(outprefix, ".data.arrow"), final_dataset)
     # Write param_files
     if batch_size !== nothing
         for (batch_id, batch) in enumerate(Iterators.partition(parameters, batch_size))
@@ -56,22 +56,13 @@ all_snps_called(found_variants::Set{<:AbstractString}, variants::Set{<:AbstractS
     variants == found_variants
 
 """
-    genotypes_encoding(variant; asint=true)
+    genotypes_encoding(variant)
 
-If asint is true then the number of minor alleles is reported, otherwise string genotypes are reported.
+String genotypes are reported.
 """
-function genotypes_encoding(variant; asint=true)
-    minor = minor_allele(variant)
+function genotypes_encoding(variant)
     all₁, all₂ = alleles(variant)
-    if asint
-        if all₁ == minor
-            return [2, 1, 0]
-        else
-            return [0, 1, 2]
-        end
-    else
-        return [all₁*all₁, all₁*all₂, all₂*all₂]
-    end
+    return [all₁*all₁, all₁*all₂, all₂*all₂]
 end
 
 NotAllVariantsFoundError(found_snps, snp_list) = 
@@ -83,7 +74,7 @@ NotBiAllelicOrUnphasedVariantError(rsid) = ArgumentError(string("Variant: ", rsi
 
 This function assumes the UK-Biobank structure
 """
-function call_genotypes(bgen_prefix::String, variants::Set{<:AbstractString}, threshold::Real; asint=true)
+function call_genotypes(bgen_prefix::String, variants::Set{<:AbstractString}, threshold::Real)
     chr_dir_, prefix_ = splitdir(bgen_prefix)
     chr_dir = chr_dir_ == "" ? "." : chr_dir_
     genotypes = nothing
@@ -102,7 +93,7 @@ function call_genotypes(bgen_prefix::String, variants::Set{<:AbstractString}, th
                         continue
                     end
                     minor_allele_dosage!(bgenfile, variant)
-                    variant_genotypes = genotypes_encoding(variant; asint=asint)
+                    variant_genotypes = genotypes_encoding(variant)
                     probabilities = probabilities!(bgenfile, variant)
                     size(probabilities, 1) != 3 && throw(NotBiAllelicOrUnphasedVariantError(rsid_))
                     chr_genotypes[!, rsid_] = call_genotypes(probabilities, variant_genotypes, threshold)

src/tmle_inputs/tmle_inputs.jl

@@ -1,4 +1,5 @@
 [deps]
+Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45"
 BGEN = "6db4b851-9beb-4b83-9d64-eb1cfb37721d"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 CategoricalArrays = "324d7699-5711-5eae-9e2f-1d82baa6b597"

test/Project.toml

@@ -1,15 +1,15 @@
 Parameters:
   - type: IATE
     target: "BINARY_1"
-    treatment: (RSID_2 = (case = 1, control = 0), TREAT_1 = (case = 1, control = 0))
+    treatment: (RSID_2 = (case = "AA", control = "GG"), TREAT_1 = (case = 1, control = 0))
     confounders: []
   - type: ATE
     target: "CONTINUOUS_2"
-    treatment: (RSID_2 = (case = 1, control = 0),)
+    treatment: (RSID_2 = (case = "AA", control = "GG"),)
     confounders: [22001]
     covariates: [COV_1, 21003]
   - type: CM
     target: "*"
-    treatment: (RSID_2 = 1, )
+    treatment: (RSID_2 = "AA", )
     confounders: [22001]
     covariates: [COV_1, 21003]