plink 进行PCA分析

当我们进行群体分析时，获得vcf文件后，可以根据变异位点对这些样本进行PCA分析，现简单介绍

1、软件安装

conda install plink

2、简单操作

本次使用train.vcf.gz 作为例子

• 第一步：将vcf转换为plink格式

plink --vcf F_M_trans.recode.vcf.gz --recode --out testacc --const-fid --allow-extra-chr
  
 # --vcf vcf 或者vcf.gz
 # --recode 输出格式
 # --out 输入前缀
 # --const-fid  添加群体信息
 # --allow-extra-chr 允许非标准染色体编号

得到3个以.map, .nosex, .ped结尾的文件。

• 基于.ped生成一个bed文件（二进制文件）

plink --allow-extra-chr --file testacc --noweb --make-bed --out testacc
 
# --file .ped + .map 文件前缀
# --make-bed 建立一个新的二进制文件

得到2个以.bim, .bed结尾的文件。

• 进行PCA分析

plink --allow-extra-chr --threads 20 -bfile testacc --pca 20 --out testacc

 
# --threads 线程数
# --pca 主成分

得到2个以.eigenval, .eigenvec结尾的文件;其中.eigenval代表每个PCA所占的比重，另外一个记录特征向量，用于坐标轴的绘制

• 可视化
  将各个样本所在的群体以及样本名字，PCA1，PCA2整理在一个文本中；如下：

pop     name    pca1    pca2
B201    B201-1  -0.621307       -0.0873675
B201    B201-2  -0.374507       -0.00563367
B201    B201-3  -0.422489       -0.00238888
B202    B202-1  0.228835        -0.172511
B202    B202-2  0.21213 -0.142105
B202    B202-3  0.247804        -0.202722
B203    B203-1  0.216946        -0.141285
B203    B203-3  0.204785        -0.0981921
B204    B204-1  0.136368        0.930181
B204    B204-2  0.169868        -0.0802494

(ggplot2)
pca <- read.table("pca_info.txt",sep = "\t",header = T)
ggplot(pca, aes(x=pca1,y=pca2)) +
        geom_point(aes(color=pop, shape=pop),size=1.5)+
        labs(x="PC1",y="PC2")+theme_bw()+theme(legend.title = element_blank())

结果

** 若想分析部分样本，则可以使用--remove参数，后接一个文件，其格式为： 第一列：群体编号， 第二列：样本名称，在这个例子中

echo '0\tSP23' > remove.txt
plink --remove remove.txt --allow-extra-chr -bfile testacc --pca 20 --out testacc_dele