Enable priors from Constrained Genomic Regression, and update version number to 1.0.0

Project.toml

@@ -1,19 +1,27 @@
 name = "LocalAncestry"
 uuid = "9c717bbe-74a8-4a1b-8702-c9261a6bba1a"
 authors = ["Bjarke G. Poulsen", "Emre Karaman"]
-version = "1.0.0-DEV"
+version = "1.0.0"
 
 [deps]
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+NLopt = "76087f3c-5699-56af-9a33-bf431cd00edd"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 VariantCallFormat = "28eba6e3-a997-4ad9-87c6-d933b8bca6c1"
 
 [compat]
 CSV = "0.10.15"
 DataFrames = "1.7.0"
+LinearAlgebra = "1.11.0"
+NLopt = "1.1.3"
 OrderedCollections = "1.8.0"
+Statistics = "1.11.1"
+StatsBase = "0.34.5"
 Tables = "1.12.0"
 VariantCallFormat = "0.5.6"
 julia = "1.6.7"

src/LocalAncestry.jl

@@ -1,5 +1,5 @@
 module LocalAncestry
-using VariantCallFormat, CSV, DataFrames, OrderedCollections, Tables
+using VariantCallFormat, CSV, DataFrames, OrderedCollections, Tables, Statistics, NLopt, StatsBase, LinearAlgebra
 
 # Write your package code here.
 include("assignMissing.jl")

src/misc.jl

@@ -109,6 +109,10 @@ function getPopulationDictionary(x)
     return y
 end
 
+function getPopulations(x::Vector{String})::Vector{String}
+    return unique(x)
+end
+
 function alleleFrequencies(x, y)
     pops = getPopulations(y)
     p = zeros(Float32, size(x, 2), length(pops))

src/origins.jl

@@ -16,7 +16,7 @@ function origins(chromosome, reference_path, target_path, referenceOrigins, orig
     haplotypeLibrary, nHaplotypeBlocks = getHaploBlocks(minHaploSize, incHaploSize, haploCrit, referenceData, popDict, 1)
 
     # Get priors
-    priorProb, priorLevel = getPriors(referenceOriginsVector, referenceData, targetIndividuals, targetData, originPriors)
+    priorProb, priorLevel = getPriors(referenceOriginsVector, referenceData, targetIndividuals, targetData, originPriors, haplotypeLibrary, ploidity)
 
     # Get log-likelihoods
     LL = calculateBlockFrequencies(haplotypeLibrary, referenceData, popDict)
@@ -28,19 +28,20 @@ function origins(chromosome, reference_path, target_path, referenceOrigins, orig
     return postProb, postClass, haplotypeLibrary
 end
 
-function getPriors(referenceOriginsVector, referenceData, targetIndividuals, targetData, originPriors)
+function getPriors(referenceOriginsVector, referenceData, targetIndividuals, targetData, originPriors, haplotypeLibrary, ploidity)
 
     if originPriors == "flat"
         x, n = priorsFlat(referenceOriginsVector, targetIndividuals)
     elseif originPriors[1:3] == "CGR"
-        x, n = priorsCGR(referenceData, targetData, targetIndividuals, referenceOriginsVector, originPriors)
+        x, n = priorsCGR(referenceData, targetData, targetIndividuals, referenceOriginsVector, originPriors, haplotypeLibrary, ploidity)
     else
         throw(DomainError(originPriors, "Expected 'flat' or 'CGR'"))
     end
 
     return x, n
 end
 
+
 function getProbabilities(predictType, targetIndividuals, ploidity, LL, populations, nHaplotypeBlocks, priorProb, priorLevel, probStayState, targetData)
 
     # Instantiate output

src/predict.jl

@@ -6,7 +6,7 @@ function predictNaiveBayes!(postProb, targetData, targetIndividuals, ploidity, L
         for h in 1:ploidity
             idname = id * "_hap" * string(h)
             for (r, reg) in enumerate(keys(LL))
-                thisprior = returnPrior(priorProb, priorLevel, id, h, r)
+                thisprior = returnPrior(priorProb, priorLevel, id, h, reg)
                 z = targetData[2*i+h-2, reg]
                 if in(z, keys(LL[reg]))
                     tmpdict = naive_bayes_predict2(thisprior, LL[reg][z], populations)

src/priors.jl

@@ -59,7 +59,7 @@ function priorsFlat(C, I)
 end
 
 # Main function for constrained genomic regression (CGR)
-function priorsCGR(referenceData, targetData, targetIndividuals, referenceOriginsVector, certainty)
+function priorsCGR(referenceData, targetData, targetIndividuals, referenceOriginsVector, certainty, haplotypeLibrary, ploidity)
     n = "block"
     populations = getPopulations(referenceOriginsVector)
     p = alleleFrequencies(referenceData, referenceOriginsVector)
@@ -68,18 +68,19 @@ function priorsCGR(referenceData, targetData, targetIndividuals, referenceOrigin
     opt = optimizerCGR(populations)
     priors = repeat([1 ./ length(populations)], length(populations))
 
-    out = performCGR(targetData, targetIndividuals, ploidity, opt, priors, certainty)
+    out = performCGR(targetData, targetIndividuals, ploidity, opt, priors, certainty, haplotypeLibrary, p, populations)
 
     return out, n
 end
 
+
 # Support function for CGR
 function optimizerCGR(x::Vector{String})
     n = length(x)
     o = NLopt.Opt(:LN_COBYLA, n)
     lower_bounds!(o, zeros(Float64, n))
     NLopt.equality_constraint!(o, (x, g) -> constraintCGR(x, g), 1e-8)
-    maxeval!(o, 1000)
+    maxeval!(o, 3000)
     return o
 end
 
@@ -99,36 +100,64 @@ function constraintCGR(x::Vector, grad::Vector)
     return sum(x) - 1
 end
 
-function certaintyScaleCGR(certainty, min_x, min_f, populations, priors)
-
-    if certainty == "CGRfull"
-        min_x = max.(log.(min_x), repeat([-10^10], length(populations)))
-    elseif certainty == "CGR"
-        min_x = max.(log.(min_f .* priors .+ (1 - min_f) .* min_x), repeat([-10^10], length(populations)))
+function certaintyScaleCGR(certainty, min_x, min_f, populations, priors, p, r, x)
+        npop = size(p,2)
+        for i in 1:npop
+            x = x - p[r,i] ./ npop 
+        end
+        Freference = mean(x .* x)
+    if certainty == "CGR"
+        scalar = 1
+    elseif certainty == "CGRdet"
+        scalar = (det(cov2cor(transpose(p[r, :] .- 0.5) * (p[r, :] .- 0.5))))*(1-min_f)
+    elseif certainty == "CGRdetsqrt"
+        scalar = (det(cov2cor(transpose(p[r, :] .- 0.5) * (p[r, :] .- 0.5))))*sqrt(1-min_f)
+    elseif certainty == "CGRdetsqd"
+        scalar = (det(cov2cor(transpose(p[r, :] .- 0.5) * (p[r, :] .- 0.5))))*(1-min_f)^2
+    elseif certainty == "CGRFdet"
+        scalar = (det(cov2cor(transpose(p[r, :] .- 0.5) * (p[r, :] .- 0.5))))*(1-min_f/Freference)
+    elseif certainty == "CGRFdetsqrt"
+        scalar = (det(cov2cor(transpose(p[r, :] .- 0.5) * (p[r, :] .- 0.5))))*sqrt(1-min_f/Freference)
+    elseif certainty == "CGRFdetsqd"
+        scalar = (det(cov2cor(transpose(p[r, :] .- 0.5) * (p[r, :] .- 0.5))))*(1-min_f/Freference)^2
+    elseif certainty == "CGRfull"
+        scalar = 1-min_f
     elseif certainty == "CGRsqrt"
-        scalar = sqrt(min_f)
-        min_x = max.(log.(scalar .* priors .+ (1 - scalar) .* min_x), repeat([-10^10], length(populations)))
+        scalar = sqrt(1-min_f)
+    elseif certainty == "CGRF"
+        scalar = 1 - min_f/Freference
     else
-        throw(DomainError(certainty, "expects 'CGRfull' or 'CGR'"))
+        throw(DomainError(certainty, "expects 'CGR', 'CGRdet', 'CGRdetsqrt', 'CGRdetsqd', 'CGRfull', or 'CGRsqrt'"))
     end
 
-    return min_x
+    return max.(log.((1-scalar) .* priors .+ (scalar) .* min_x), repeat([-10^10], length(populations)))
 end
 
-function performCGR(targetData, targetIndividuals, ploidity, opt, priors, certainty)
+p = zeros(Float64, 6,3)
+
+
+function performCGR(targetData, targetIndividuals, ploidity, opt, priors, certainty, haplotypeLibrary, p, populations)
     out = Dict{String,Dict{String,Float64}}()
-    for (i, ind) in enumerate(targetIndividuals)
-        for h in 1:ploidity
-            outname = ind * "_hap" * string(h)
-            inindice = 2 * i + h - 2
-            y = (targetData')[:, inindice]
-            NLopt.min_objective!(opt, (x, g) -> objectiveCGR(x, g, y, p))
+    for r in keys(haplotypeLibrary)
+
+        blocks = unique(targetData[:,r], dims = 1)
+        for (ib, b) in enumerate(eachrow(blocks))
+            NLopt.min_objective!(opt, (x, g) -> objectiveCGR(x, g, blocks[ib,:], p[r, :]))
             min_f, min_x, _ = NLopt.optimize(opt, priors)
 
-            out[outname] = Dict(zip(populations, certaintyScaleCGR(certainty, min_x, min_f, populations, priors)))
+                blockmatch = [itb for (itb, tb) in enumerate(eachrow(targetData[:,r])) if all(b .== targetData[itb,r])] 
+                for bm in blockmatch
+                    ind = ceil( bm / ploidity) 
+                    h = bm - ind * ploidity + ploidity
+                outname = targetIndividuals[convert(Int64,ind)] * "_hap" * string(convert(Int8,h)) * "_reg" * string(r)
+
+                out[outname] = Dict(zip(populations, certaintyScaleCGR(certainty, min_x, min_f, populations, priors, p, r, blocks[ib,:])))
+                end
         end
     end
+
     return out
 end
 
 
+