前期准备

给标记加上ID

SNP data通常都是以VCF格式文件呈现，拿到VCF文件的第一件事情就是添加各个SNP位点的ID。
先看一下最开始生成的VCF文件：

image

可以看到，ID列都是"."，需要我们自己加上去。我用的是某不知名大神写好的perl脚本，可以去我的github上下载，用法：

perl path2file/VCF_add_id.pl YourDataName.vcf YourDataName-id.vcf`

当然也可以用excel手工添加。添加后的文件如下图所示（格式：CHROMID__POS）：

image

SNP位点过滤（Missing rate and maf filtering）

SNP位点过滤前需要问自己一个问题，我的数据需要过滤吗？

一般要看后期是否做关联分析（GWAS）；如果只是单纯研究群体结构建议不过滤，因为过滤掉低频位点可能会改变某些样本之间的关系；如果需要和表型联系其来做关联分析，那么建议过滤，因为在后期分析中低频位点是不在考虑范围内的，需要保持前后一致。

如果过滤，此处用到强大的plink软件，用法：

plink --vcf YourDataName-id.vcf --maf 0.05 --geno 0.2 --recode vcf-iid -out YourDataName-id-maf0.05 --allow-extra-chr

参数解释：--maf 0.05：过滤掉次等位基因频率低于0.05的位点；--geno 0.2：过滤掉有20%的样品缺失的SNP位点；--allow-extra-chr：我的参考数据是Contig级别的，个数比常见分析所用的染色体多太多，所以需要加上此参数。

格式转换

将vcf文件转换为bed格式文件。
这里注意一点！！！！：应该是软件的问题，需要把染色体/contig名称变成连续的数字（1 to n），不然会报错无法算出结果！（坑）

plink --vcf YourDataName-id-maf0.05.vcf --make-bed --out snp --chr-set 29 no-xy

参数解释：--chr-set 给出染色体/contig的数目；no-xy 没有xy染色体。

用gcta做PCA分析

gcta输出grm阵列（genetic relationship matrix）

gcta64 --make-grm --out snp.gcta --bfile snp --autosome-num 29

参数解释：--autosome-num常染色体数目。

gcta计算PCA

gcta64 --grm snp.gcta --pca 20 --out snp.gcta

参数解读：--pca 20 保留前20个PCA。

特征值结果储存在snp.gcta.eigenval中，特征向量储存在snp.gcta.eigenvec中。

结果处理

将特征值结果和特征向量结果用R处理为可读性结果。写好的R包我放在了Github中：PCA2normal_format.R，大家自行下载使用。

如果不想下载，直接复制如下代码：

eigvec <- read.table("snp.gcta.eigenvec", header = F, stringsAsFactors = F)
write.table(eigvec[2:ncol(eigvec)], file = "gcta.eigenvector.xls", sep = "\t", row.names = F, col.names = T, quote = F)

eigval <- read.table("snp.gcta.eigenval", header = F)
pcs <- paste0("PC", 1:nrow(eigval))
eigval[nrow(eigval),1] <- 0
percentage <- eigval$V1/sum(eigval$V1)*100
eigval_df <- as.data.frame(cbind(pcs, eigval[,1], percentage), stringsAsFactors = F)
names(eigval_df) <- c("PCs", "variance", "proportion")
eigval_df$variance <- as.numeric(eigval_df$variance)
eigval_df$proportion <- as.numeric(eigval_df$proportion)
write.table(eigval_df, file = "gcta.eigenvalue.xls", sep = "\t", quote = F, row.names = F, col.names = T)

转换前snp.gcta.eigenvec

image

转换后gcta.eigenvector.xls

image

转换前snp.gcta.eigenval（PCA方差）

image

转换后gcta.eigenvalue.xls（PCA方差+方差解释率）

image

用LDAK做PCA分析

相比gcta，能用LD对结果进行校正，具体来说，就是先用LD计算每个SNP位点的权重，根据权重再计算Kinship，这样的结果更接近真实情况。

LDAK输出grm阵列（genetic relationship matrix）

在不考虑权重的情况下，方法如下：

ldak5.linux --calc-kins-direct snp.ldak --bfile snp --ignore-weights YES --kinship-gz YES --power -0.25

用LD计算每个SNP位点的权重，根据权重再计算Kinship

#切割
ldak5.linux --cut-weights snp.sections --bfile snp
#查看有多少个section
cat snp.sections/section.number
#根据自己的section个数分别计算权重（我这里是31个）
for section in {1..31}; do ldak5.linux --calc-weights snp.sections --bfile snp --section $section; done
#weight文件整合，给SNP赋权重值
ldak5.linux --join-weights snp.sections --bfile snp
#输出grm阵列
ldak5.linux --calc-kins-direct snp.ldak.weight --bfile snp --weights snp.sections/weights.all --kinship-gz YES --power -0.25

LDAK计算PCA（calculate PCA）

不考虑权重

ldak5.linux --pca snp.ldak --grm snp.ldak --axes 222

参数解释：--axes 样本数量，这里需要准确填写，不然无法用特征值计算方差解释率。

考虑权重

ldak5.linux --pca snp.ldak.weight --grm snp.ldak.weight --axes 222

特征值结果储存在snp.ldak.weight.values中，特征向量储存在snp.ldak.weight.vect中。

结果处理

和gcta方法一样（只是输入文件名称不同），用同一个R包将特征值结果和特征向量结果用R处理为可读性结果，这里不再赘述。

同样，如果不想下载，直接复制如下代码：

eigvec <- read.table("snp.ldak.weight.vect", header = F, stringsAsFactors = F)
colnames(eigvec) <- c("FID", "Sample", paste0("PC", 1:20))
write.table(eigvec[2:ncol(eigvec)], file = "ldak.eigenvector.xls", sep = "\t", row.names = F, col.names = T, quote = F)

eigval <- read.table("snp.ldak.weight.values", header = F)
pcs <- paste0("PC", 1:nrow(eigval))
eigval[nrow(eigval),1] <- 0
percentage <- eigval$V1/sum(eigval$V1)*100
eigval_df <- as.data.frame(cbind(pcs, eigval[,1], percentage), stringsAsFactors = F)
names(eigval_df) <- c("PCs", "variance", "proportion")
eigval_df$variance <- as.numeric(eigval_df$variance)
eigval_df$proportion <- as.numeric(eigval_df$proportion)
write.table(eigval_df, file = "ldak.eigenvalue.xls", sep = "\t", quote = F, row.names = F, col.names = T)

数据可视化

用R画散点图即可，散点图的画法由于篇幅原因，另外写个帖详细说明。这里直接分享一下我作图的方法：

数据准备

除了上面获得的两个文件：ldak.eigenvector.xls和ldak.eigenvalue.xls外，还需要准备一个命名为pca.pop.xls的文件，该文件包含如下内容：

image

注：第一列为排序；第二列为ID+vcf_id；第三列为vcf_id，需要和ldak.eigenvector.xls文件中SampleI ID一致；第四列为分组信息1；第五列为基于分组信息1给予的颜色信息；第六列为基于分组信息2给予的标记形状信息；第七列为分组信息2（例子中为地理来源）。

赋值

按照如下代码给各个参数赋值（赋予路径信息）：

eigvec <- "E:/Desktop/PCA/ldak.eigenvector.xls"
eigval <- "E:/Desktop/PCA/ldak.eigenvalue.xls"
popinfo <- "E:/Desktop/PCA/pca.pop.xls"
key <- "ldak_PCA"
od <- "E:/Desktop/PCA"

注：我在桌面建了一个名为PCA的文件夹，把ldak.eigenvector.xls、ldak.eigenvalue.xls和pca.pop.xls三个文件都放在里面；key是指定输出文件的文件名前缀；od是指定输出文件存放目录。

调用pca.plot2d.r包作图

poptable <- read.table(popinfo, header = T, comment.char = "")
pop <- unique(poptable[,4:7])
print(pop)
source("pca.plot2d.r")
pca_plot(eigenvector = eigvec, eigenvalue = eigval,
         group = popinfo, key = key, outdir = od,
         shape = T, shapes = pop$pch, border = T, border_size = 2.5,
         line0 = T, line0_size = 1)

执行后会返回，分别以PC1-PC2、PC1-PC3和PC2-PC3为坐标轴的PCA图，包括pdf和png两种格式的文件。