Compute Principal Components from GDS Genotype Data

Computes principal components from per-chromosome SeqArray GDS files for population stratification adjustment in association analyses. Variants are filtered by MAF and missing rate, LD pruned, then PCA is performed via SNPRelate.

compute_pcs_gds(
  gds_files,
  n_pcs = 10L,
  maf_threshold = 0.05,
  missing_rate = 0.05,
  ld_threshold = 0.2,
  output_file = NULL,
  num_thread = 1L,
  algorithm = c("randomized", "exact"),
  seed = NULL,
  verbose = 1
)

Arguments

gds_files

Character vector. Paths to one or more SeqArray GDS files, typically one per chromosome. All files must exist and be readable.

n_pcs

Integer. Number of principal components to compute. Must be a positive integer (\(\geq 1\)). Default: 10.

maf_threshold

Numeric. Minimum minor allele frequency for variant inclusion. Must be in \([0, 0.5]\). Default: 0.05.

missing_rate

Numeric. Maximum per-variant missing rate. Variants with missing rate above this threshold are excluded. Must be in \([0, 1]\). Default: 0.05.

ld_threshold

Numeric. LD pruning \(r^2\) threshold. Variants with pairwise \(r^2\) above this threshold are pruned. Must be in (0, 1]. Default: 0.2.

output_file

Character or NULL. Optional path to save results. Extension determines format: .rds uses saveRDS() (preserves S3 class and all attributes), .csv writes a CSV with metadata comment header. Default: NULL (no file output).

num_thread

Integer. Number of threads for SNPRelate operations. Default: 1L.

algorithm

Character. PCA algorithm: "randomized" (default) uses randomized truncated SVD—fast and memory-efficient for large datasets (\(O(npk)\) time, \(O(nk)\) memory). "exact" computes the full \(n \times n\) GRM then extracts eigenvalues—requires \(O(n^2)\) memory (e.g., ~18 GB for 48K samples). Use "exact" only for small datasets or when exact eigendecomposition is needed. Default: "randomized".

seed

Integer or NULL. Random seed for reproducibility of LD pruning. Default: NULL (no seed set).

verbose

Integer. Verbosity level: 0 = silent, 1 = progress messages, 2 = detailed messages. Default: 1.

Value

An object of class c("glow_pcs", "data.frame") with columns:

sample.id

Character. Sample identifiers from the GDS file(s).

PC1, PC2, ..., PC{n_pcs}

Numeric. Principal component scores.

The following attributes are attached:

eigenvalues

Numeric vector of length n_pcs.

variance_proportion

Numeric vector. Proportion of variance explained by each PC.

total_variance_explained

Numeric scalar. Sum of variance proportions for all returned PCs.

n_variants_used

Integer. Number of variants after filtering and LD pruning.

n_variants_per_chr

Named integer vector. Pruned variant count per chromosome.

call_args

List. Arguments used for this call (for reproducibility).

Details

Processing Pipeline:

1. Validate inputs and check SNPRelate/SeqArray availability

2. Merge chromosome GDS files via SeqArray::seqMerge() (if multiple files; skipped for single file)

3. Convert SeqArray GDS to SNP GDS via SeqArray::seqGDS2SNP()

4. LD prune via SNPRelate::snpgdsLDpruning() with MAF and missing rate filtering

5. Compute PCA via SNPRelate::snpgdsPCA() on pruned variants

6. Build glow_pcs S3 object with eigenvalue diagnostics

7. Save to disk if output_file is specified

8. Clean up temporary files (via on.exit())

GDS Format Compatibility:

SNPRelate functions (snpgdsPCA, snpgdsLDpruning) require SNP GDS format (SNPGDSFileClass), not SeqArray format (SeqVarGDSClass). This function handles the conversion automatically via seqGDS2SNP().

Computational Complexity:

For \(n\) samples and \(p\) variants after pruning, PCA computation is \(O(n^2 p)\) for the covariance matrix. LD pruning is \(O(p \cdot w)\) where \(w\) is the sliding window size.

Dependencies:

SNPRelate and SeqArray are optional dependencies (in Suggests). The function checks availability at runtime and provides installation instructions if missing.

References

Zheng X, Levine D, Shen J, Gogarten SM, Laurie C, Weir BS (2012). A high-performance computing toolset for relatedness and principal component analysis of SNP data. Bioinformatics, 28(24), 3326-3328. DOI: 10.1093/bioinformatics/bts606

See also

plink_to_gds for converting PLINK files to GDS format

Examples

if (FALSE) { # \dontrun{
# Example 1: Basic PC computation from per-chromosome GDS files
gds_files <- sprintf("data/gds/chr%d.gds", 1:22)
pcs <- compute_pcs_gds(gds_files, n_pcs = 10)
print(pcs)

# Example 2: Compute and save with seed for reproducibility
pcs <- compute_pcs_gds(
  gds_files = gds_files,
  n_pcs = 20,
  output_file = "results/pcs.rds",
  seed = 42L,
  num_thread = 4L,
  verbose = 1
)

# Example 3: Scree plot to determine optimal number of PCs
plot(pcs, type = "scree")

# Example 4: PC scatter matrix for population structure
plot(pcs, type = "pairs", n_pairs = 6)

# Example 5: Single GDS file (skips merge step)
pcs <- compute_pcs_gds("data/merged.gds", n_pcs = 10)
} # }