手动写PCA代码的过程，使用参考群和候选群分开的形式，这样新样本无需进行特征分解。

目的

一般软件（如 plink, gcta等）计算 PCA 的过程的如下：

1. 构建 G 阵
2. 对 G 阵进行特征值分解

得到的特征向量就是需要的主成分。

当样本数较多时，PCA计算速度慢；而且当仅仅是新增部分样本时，均要使用全部样本重新计算一遍主成分，计算效率低。本脚本将做 PCA 的数据分为两个群体，参考群（旧数据）和候选群（新增数据）；对于参考群构建 SNP 的协方差矩阵并进行特征值分解；对于候选群，可以直接使用参考群的特征值分解结果，无需再次进行特征值分解，因此计算效率大大提高。因此本脚本适用于同一个群体不断新增样本的情况，如不更新参考群，则全部 PCA 分析只需要进行一次特征值分解。

过程

参考群

1. plink 转化为 012 格式（保存计数的碱基，后面还要用）

2. 遍历 raw 文件，缺失 (NA) 改为 5 。

   缺失改为5的好处在于，2p_j 的范围在 [0,2], 因此 M_ij - 2p_j 的范围在 [-2,2]，当缺失为 5时，其范围为 [3,5]，这样可以和正常值分开。

3. 构建 SNP 的协方差矩阵并分解，过程如下：

   • 生成基因型矩阵（一行一个样本，一列一个位点）

   • 计算 maf (去除缺失位点)

   • 每一列减去 2 * maf （中心化）

   • 将此时矩阵中大于 2 的元素设为 0 （缺失位点）

   • 计算每一列 SNP 的标准差（2 * maf * (1-maf)）

   • 每一列除以相应标准差，得到 X 矩阵（标准化）

   • SNP 的协方差矩阵计算公式如下 \[ G = \frac{1}{m} X^{T}X \]

   • 对 G 阵进行特征值分解，得到特征向量 u （可以提取前 n 个）和特征值 sigma

   • 通过计算 \(Xu\) 得到我们需要的主成分。

候选群

候选群必须包含参考群中的全部位点。

1. plink 转化为 012 格式（注意计数的碱基必须和参考群一致，使用 --recode-allele {ref_prefix}.recode_allele 参数设置这一点 ）

2. 遍历 raw 文件，缺失 (NA) 改为 5 。

3. 计算特征向量，过程如下：

   • 生成候选群的基因型矩阵（一行一个样本，一列一个位点）

   • 使用 参考群的 maf

   • 每一列减去 2 * maf （中心化）

   • 将此时矩阵中大于 2 的元素设为 0 （缺失位点）

   • 计算每一列 SNP 的标准差（2 * maf * (1-maf)）

   • 每一列除以相应标准差，得到 X 矩阵（标准化）

   • 通过计算 \(Xu\) （\(u\) 为参考群得到的特征向量）得到我们需要的主成分。

pca画图

可选

参考资料

1. https://vincere.fun/posts/fc532062/
2. https://vincere.fun/posts/1e12f93e/