当我们进行群体分析时，获得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
