PCAone Tutorial

For most modern CPUs and Linux systems, download the one named with avx2.

pkg=https://github.com/Zilong-Li/PCAone/releases/latest/download/PCAone-avx2-Linux.zip
wget $pkg || curl -LO $pkg
unzip -o PCAone-avx2-Linux.zip

If the server is too old to support avx2 instruction, download the alternative version.

pkg=https://github.com/Zilong-Li/PCAone/releases/latest/download/PCAone-x64-Linux.zip
wget $pkg || curl -LO $pkg
unzip -o PCAone-x64-Linux.zip

pkg=http://popgen.dk/zilong/datahub/pca/example.tar.gz
wget $pkg || curl -LO $pkg
tar -xzf example.tar.gz && rm -f example.tar.gz

You should find a fold named example with some example data.

There are compiled binaries provided for both Linux and Mac platform. Check the releases page to download one.

PCAone is also available from bioconda.

conda config --add channels bioconda
conda install pcaone
PCAone --help

PCAone can be running on a normal computer/laptop with x86-64 instruction set architecture. PCAone has been tested on both Linux and MacOS system. To build PCAone from the source code, the following dependencies are required:

• GCC/Clang compiler with C++11 support
• GNU make
• zlib

We recommend building the software from source with MKL as backend to maximize the performance. For MacOS users, we recommend using llvm by brew install llvm instead of the default clang shipped with MacOS. Check out the mac workflow.

Run ./PCAone --help to see all options including hidden advanced options. The below are common useful options.

Main options:
  -h, --help                     print all options including hidden advanced options
  -d, --svd arg (=2)             SVD method to be applied. default 2 is recommended for big data.
                                 0: the Implicitly Restarted Arnoldi Method (IRAM)
                                 1: the Yu's single-pass Randomized SVD with power iterations
                                 2: the accurate window-based Randomized SVD method (PCAone)
                                 3: the full Singular Value Decomposition.
  -b, --bfile arg                prefix to PLINK .bed/.bim/.fam files
  -B, --binary arg               path of binary file
  -c, --csv arg                  path of comma seperated CSV file compressed by zstd
  -g, --bgen arg                 path of BGEN file compressed by gzip/zstd
  -G, --beagle arg               path of BEAGLE file compressed by gzip
  --USV arg                      prefix to PCAone .eigvecs/.eigvals/.loadings/.mbim
  --read-U arg                   path of file with left singular vectors (.eigvecs)
  --read-V arg                   path of file with right singular vectors (.loadings)
  --read-S arg                   path of file with eigen values (.eigvals)
  -k, --pc arg (=10)             top k principal components (PCs) to be calculated
  -m, --memory arg (=0)          RAM usage in GB unit for out-of-core mode. default is in-core mode
  -n, --threads arg (=12)        the number of threads to be used
  -o, --out arg (=pcaone)        prefix to output files. default [pcaone]
  -p, --maxp arg (=40)           maximum number of power iterations for RSVD algorithm
  -S, --no-shuffle               do not shuffle columns of data for --svd 2 (if not locally correlated)
  -v, --verbose                  verbose message output
  -V, --printv                   output the right eigenvectors with suffix .loadings
  -C, --scale arg (=0)           do scaling for input file.
                                 0: do just centering
                                 1: do log transformation eg. log(x+0.01) for RNA-seq data
                                 2: do count per median log transformation (CPMED) for scRNAs
  --emu                          use EMU algorithm for genotype input with missingness
  --pcangsd                      use PCAngsd algorithm for genotype likelihood input
  --maf arg (=0)                 exclude variants with MAF lower than this value
  --match-bim arg                the .mbim file to be matched, where the 7th column is allele frequency
  --project arg (=0)             project the new samples onto the existing PCs.
                                 0: disabled
                                 1: by multiplying the loadings with mean imputation for missing genotypes
                                 2: by solving the least squares system Vx=g. skip sites with missingness
                                 3: by Augmentation, Decomposition and Procrusters transformation
  --inbreed arg (=0)             compute the inbreeding coefficient accounting for population structure.
                                 0: disabled
                                 1: compute per-site inbreeding coefficient and HWE test
  --ld                           output a binary matrix for downstream LD related analysis
  --ld-r2 arg (=0)               r2 cutoff for LD-based pruning. (usually 0.2)
  --ld-bp arg (=0)               physical distance threshold in bases for LD. (usually 1000000)
  --ld-stats arg (=0)            statistics to calculate LD r2 for pairwise SNPs.
                                 0: the ancestry adjusted, i.e. correlation between residuals
                                 1: the standard, i.e. correlation between two alleles
  --print-r2                     print LD r2 to *.ld.gz file for pairwise SNPs within a window
  --clump arg                    assoc-like file with target variants and pvalues for clumping
  --clump-names arg (=CHR,BP,P)  column names in assoc-like file for locating chr, pos and pvalue
  --clump-p1 arg (=0.0001)       significance threshold for index SNPs
  --clump-p2 arg (=0.01)         secondary significance threshold for clumped SNPs
  --clump-r2 arg (=0.5)          r2 cutoff for LD-based clumping
  --clump-bp arg (=250000)       physical distance threshold in bases for clumping

\newpage

This depends on your datasets, particularlly the relationship between number of samples (N) and the number of variants / features (M) and the top PCs (k). Here is an overview and the recommendation.

PCAone is designed to be extensible to accept many different formats. Currently, PCAone can work with SNP major genetic formats to study population structure. such as PLINK, BGEN and Beagle. Also, PCAone supports a comma delimited CSV format compressed by zstd, which is useful for other datasets requiring specific normalization such as single cell RNAs data.

PCAone has both in-core and out-of-core mode for 3 different partial SVD algorithms, which are IRAM (--svd 0), Yu+Halko RSVD (--svd 1) and PCAone window-based RSVD (--svd 2). Also, PCAone supports full SVD (--svd 3) but with only in-core mode. Therefore, there are 7 ways for doing PCA in PCAone. In default PCAone uses in-core mode with --memory 0, which is the fastest way to do calculation. However, in case the server runs out of memory with in-core mode, the user can trigger out-of-core mode by specifying the amount of memory using --memory option with a value greater than 0.

PCAone will automatically apply the standard normalization for genetic data. Additionally, there are 3 different normalization method implemented with --scale option.

• 0: do just centering by substracting the mean
• 1: do log transformation (usually for count data, such as bulk RNA-seq data)
• 2: do count per median log transformation (usually for single cell RNA-seq data)

One should choose proper normalization method for specific type of data.

Project new samples onto existing PCs is supported with --project option. First, we run PCAone on a set of reference samples and output the loadings:

PCAone -b ref_samples -k 10 --printv -o ref

Then, we need to read in the SNPs loadings from the ref set (--read-V) and its scaling factors (--read-S) as well as the allele frequencies form the .mbim file via --match-bim. Note: one can use the --USV option instead to simplify the usage since v0.4.8 Here is the example command to project new target samples and get the PCs coordinates.

PCAone -b new_samples \
       --USV ref \  ## prefix to .eigvecs, .eigvals, .loadings
       --project 2 \  ## check the manual on projection methods
       -o new

To test Hardy-Weinberg equilibrium in presence of population structure, we need to work on the so-called individual allele frequencies matrix \(\Pi\), which can be reconstructed via the output of PCAone, i.e the .eigvecs,.eigvals,.loadings and .mbim files, generated by

PCAone -b example/plink -k 3 -V -o pcaone

Then we apply --inbreed 1 option to obtain the P value of HWE and inbreeding coefficient per-site. The output file is named with suffix .hwe.

PCAone -b example/plink \
       --USV pcaone \  
       --inbreed 1 \ 
       -o inbreed

LD patterns vary across diverse ancestry and structured groups, and conventional LD statistics, e.g. the implementation in plink --ld, failed to model the LD in admixed populations. Thus, we can use the so-called ancestry-adjusted LD statistics to account for population structure in LD. See our paper for more details.

To calculate the ancestry-adjusted LD matrix, we first figure out the number of principal components (-k/--pc) that capture population structure. In this example, assuming that 3 PCs can accout for population structure, we enable --ld option to calculate and output the ancestry adjusted LD matrix in a file with suffix .residuals.

./PCAone -b example/plink -k 3 --ld -o adj

Currently, the LD R2 for pairwise SNPs within a window can be outputted via --print-r2 option.

./PCAone -B adj.residuals \
         --match-bim adj.mbim \
         --ld-bp 1000000 \
         --print-r2 \
         -o adj

We provide the calc_decaybin.R script to parse the output file .ld.gz and calculate the average R2 for each bin as well as plotting. We also provide the nextflow ld.nf for benchmarking the LD statistics.

Given the LD binary file .residuals and its associated variant file .mbim, we can do pruning based on user-defined thresholds and windows

./PCAone -B adj.residuals \
         --match-bim adj.mbim \
         --ld-r2 0.8 \
         --ld-bp 1000000 \
         -o adj

Likewise, we can do clumping based on the Ancestry-Adjusted LD matrix and user-defined association results

./PCAone -B adj_ld.residuals \
         --match-bim adj.mbim \
         --clump example/plink.pheno0.assoc,example/plink.pheno1.assoc  \
         --clump-p1 0.05 \
         --clump-p2 0.01 \
         --clump-r2 0.1 \
         --clump-bp 10000000 \
         -o adj

We want to compute the top 40 PCs for this genotype dataset using 20 threads and only 2 GBs memory. We will use the proposed window-based RSVD algorithm with default setting --svd 2.

./PCAone --bfile example/plink -k 40 -m 2 -n 20

Then, we can make a PCA plot in R.

pcs <- read.table("pcaone.eigvecs2",h=F)
plot(pcs[,1:2+2], col=factor(pcs[,1]), xlab = "PC1", ylab = "PC2")
legend("topright", legend=unique(pcs[,1]), col=1:4, pch = 21, cex=1.2)

Imputation tools usually generate the genotype probabilities or dosages in BGEN format. To do PCA with the imputed genotype probabilities, we can work on BGEN file with --bgen option instead.

./PCAone --bgen example/test.bgen -k 10 -n 4 -m 2

Then, we can make a PCA plot in R.

pcs <- read.table("pcaone.eigvecs",h=F)
pop <- read.table("example/plink.fam",h=F)[,1]
plot(pcs[,1:2], col=factor(pop), xlab = "PC1", ylab = "PC2")
legend("topright", legend=unique(pop), col=factor(unique(pop)), pch = 21, cex=1.2)

In this example, we run PCA for the single cell RNAs-seq data using a different input format with a normalization method called count per median log transformation (CPMED).

./PCAone --csv example/BrainSpinalCord.csv.zst -k 10 -n 20 -m 4 --scale 2 --svd 1

It should take around 5 minutes.