Add RareEpistasisMeta package for rare variant epistasis meta-analysis - #16
Add RareEpistasisMeta package for rare variant epistasis meta-analysis#16meibujun wants to merge 3 commits into
Conversation
There was a problem hiding this comment.
💡 Codex Review
Here are some automated review suggestions for this pull request.
ℹ️ About Codex in GitHub
Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you
- Open a pull request for review
- Mark a draft as ready
- Comment "@codex review".
If Codex has suggestions, it will comment; otherwise it will react with 👍.
Codex can also answer questions or update the PR. Try commenting
| function load_vcf(path::AbstractString; impute::Symbol=:mean) | ||
| open_stream(path) do io | ||
| samples = String[] | ||
| variants = DataFrame() | ||
| genotype_rows = Vector{Vector{Float64}}() | ||
| for line in eachline(io) | ||
| startswith(line, "##") && continue | ||
| if startswith(line, "#CHROM") | ||
| header = split(line, '\t') | ||
| samples = header[10:end] | ||
| variants = DataFrame(CHROM = String[], POS = Int[], ID = String[], REF = String[], ALT = String[]) | ||
| continue | ||
| end | ||
| fields = split(line, '\t') | ||
| push!(variants, (fields[1], parse(Int, fields[2]), fields[3], fields[4], fields[5])) | ||
| format_fields = split(fields[9], ':') | ||
| gt_index = findfirst(==("GT"), format_fields) | ||
| gt_index = isnothing(gt_index) ? 1 : gt_index | ||
| gts = map(x -> parse_gt(x, gt_index), fields[10:end]) | ||
| push!(genotype_rows, gts) | ||
| end | ||
| G = reduce(hcat, genotype_rows)' | ||
| impute_missing!(G; method = impute) | ||
| return G, variants, DataFrame(Sample = samples) |
There was a problem hiding this comment.
Avoid transposing VCF genotypes before imputation
In load_vcf the collected genotype vectors are n_samples × n_variants once reduce(hcat, genotype_rows) completes, but the subsequent transpose (') flips the matrix to n_variants × n_samples. As a result, impute_missing! iterates over sample columns rather than variant columns and every downstream routine that assumes rows are samples (e.g. QC and statistical methods) receives a matrix with axes swapped. This mis‑orientation will lead to dimension mismatches and incorrect imputation whenever a VCF file is read.
Useful? React with 👍 / 👎.
| for j in axes(G, 2) | ||
| column = view(G, :, j) | ||
| miss = count(isnan, column) / length(column) | ||
| if miss > missing_rate | ||
| keep[j] = false | ||
| continue | ||
| end | ||
| valid = column[.!isnan.(column)] | ||
| if isempty(valid) | ||
| keep[j] = false | ||
| continue | ||
| end | ||
| maf = mean(valid) / 2 | ||
| if maf < maf_threshold | ||
| keep[j] = false |
There was a problem hiding this comment.
Compute minor allele frequency instead of mean allele count
The variant filter treats maf = mean(valid) / 2, which corresponds to the major allele frequency when the alternate allele is common. For monomorphic variants with all alternate alleles (frequency ≈1), this code returns maf≈1 and the variant passes the maf_threshold despite having zero minor‑allele variation. The filter therefore retains sites that should be discarded and can distort downstream analyses. maf should use min(p, 1-p) where p = mean(valid)/2 to reflect the true minor allele frequency.
Useful? React with 👍 / 👎.
Summary
Testing
julia --project=. -e 'using Pkg; Pkg.instantiate(); Pkg.test()'(fails: Julia binary unavailable in container)https://chatgpt.com/codex/tasks/task_e_68ee1dd547dc83219041127ab9e76482