Skip to content

Commit 8bf6cf8

Browse files
gaowclaude
andcommitted
Fix nSample parameter bug and enforce exact N-in = N-out row counts
Root cause: the R implementation was using the GWAS sample size (e.g. 5000) for the nSample parameter, but the original DENTIST binary always uses the LD reference panel size (ref.N from the .fam file, e.g. 165). This parameter controls the SVD truncation rank K = min(idx.size(), nSample) * propSVD, so using the wrong value produces completely different imputation results (imputed_z correlation of only 0.40 instead of 1.0). With nSample correctly set to the reference panel size (NULL default in dentist_from_files), R now matches the binary to machine precision on both the 2500-variant subset and full 17421-variant toy dataset: - imputed_z max diff: 5e-6, correlation: 1.0 - rsq max diff: 5e-7, correlation: 1.0 - duplicate detection: 100% agreement Changes: - Update nSample docs in dentist(), dentist_single_window(), and dentist_from_files() to clarify it is the reference panel size, not the GWAS sample size - Remove incorrect filter in merge_windows() that dropped variants with imputed_z == 0 & rsq == 0, which caused output to have fewer rows than input - Rewrite unit tests (25 -> 43) with exact row count assertions (expect_equal(nrow(res), N)) instead of allowing ranges - Add new tests: windowed mode exact row counts, outlier_stat formula verification, X matrix input, fill region coverage, column checks, original_z preservation Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
1 parent 0446c35 commit 8bf6cf8

5 files changed

Lines changed: 174 additions & 109 deletions

File tree

R/dentist_qc.R

Lines changed: 10 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -46,8 +46,9 @@ resolve_LD_input <- function(R = NULL, X = NULL, nSample = NULL, need_nSample =
4646
#' @param R Square LD correlation matrix. Provide either \code{R} or \code{X}.
4747
#' @param X Genotype matrix (samples x SNPs). If provided, LD is computed via
4848
#' \code{compute_LD(X)} and \code{nSample} defaults to \code{nrow(X)}.
49-
#' @param nSample The number of samples. Required when \code{R} is provided; inferred
50-
#' from \code{X} when \code{X} is provided.
49+
#' @param nSample The number of samples in the LD reference panel (NOT the GWAS sample
50+
#' size). This controls the SVD truncation rank K = min(idx_size, nSample) * propSVD.
51+
#' Required when \code{R} is provided; inferred from \code{X} when \code{X} is provided.
5152
#' @param window_size The size of the window for dividing the genomic region. Default is 2000000.
5253
#' @param pValueThreshold The p-value threshold for significance. Default is 5e-8.
5354
#' @param propSVD The proportion of singular value decomposition (SVD) to use. Default is 0.4.
@@ -154,8 +155,9 @@ dentist <- function(sum_stat, R = NULL, X = NULL, nSample = NULL,
154155
#' @param R Square LD correlation matrix. Provide either \code{R} or \code{X}.
155156
#' @param X Genotype matrix (samples x SNPs). If provided, LD is computed via
156157
#' \code{compute_LD(X)} and \code{nSample} defaults to \code{nrow(X)}.
157-
#' @param nSample Number of samples. Required when \code{R} is provided; inferred
158-
#' from \code{X} when \code{X} is provided.
158+
#' @param nSample Number of samples in the LD reference panel (NOT the GWAS sample
159+
#' size). Controls the SVD truncation rank. Required when \code{R} is provided;
160+
#' inferred from \code{X} when \code{X} is provided.
159161
#' @param pValueThreshold P-value threshold for outlier detection. Default is 5e-8.
160162
#' @param propSVD SVD truncation proportion. Default is 0.4.
161163
#' @param gcControl Logical; apply genomic control. Default is FALSE.
@@ -531,10 +533,6 @@ merge_windows <- function(dentist_result_by_window, window_divided_res) {
531533
filter(index_global >= window_divided_res$fillStartIdx[k] & index_global < window_divided_res$fillEndIdx[k])
532534
merged_results <- rbind(merged_results, extracted_results)
533535
}
534-
# Filter out un-imputed variants (imputed_z == 0 and rsq == 0) after merging,
535-
# matching the original DENTIST binary which skips these when writing output.
536-
merged_results <- merged_results %>%
537-
filter(!(imputed_z == 0 & rsq == 0))
538536
return(merged_results)
539537
}
540538

@@ -593,8 +591,10 @@ read_dentist_sumstat <- function(gwas_summary) {
593591
#' @param gwas_summary Path to the GWAS summary statistics file (DENTIST format:
594592
#' tab-separated, 8 columns with header: SNP A1 A2 freq beta se p N). May be gzipped.
595593
#' @param bfile PLINK binary file prefix (expects .bed/.bim/.fam files).
596-
#' @param nSample Sample size. If NULL, taken from the N column of the summary file
597-
#' (using the median value). Default is NULL.
594+
#' @param nSample Number of samples in the LD reference panel. If NULL (recommended),
595+
#' uses the reference panel size from the genotype matrix. This controls the SVD
596+
#' truncation rank K = min(idx_size, nSample) * propSVD. Note: this is the reference
597+
#' panel size, NOT the GWAS sample size. Default is NULL.
598598
#' @param window_size Window size in base pairs. Default is 2000000.
599599
#' @param pValueThreshold P-value threshold for outlier detection. Default is 5.0369e-8.
600600
#' @param propSVD SVD truncation proportion. Default is 0.4.

man/dentist.Rd

Lines changed: 11 additions & 5 deletions
Some generated files are not rendered by default. Learn more about customizing how changed files appear on GitHub.

man/dentist_from_files.Rd

Lines changed: 4 additions & 2 deletions
Some generated files are not rendered by default. Learn more about customizing how changed files appear on GitHub.

man/dentist_single_window.Rd

Lines changed: 3 additions & 2 deletions
Some generated files are not rendered by default. Learn more about customizing how changed files appear on GitHub.

0 commit comments

Comments
 (0)