使用smartpca做种群遗传学中常用的PCA分析实例

本文中使用到的数据是 文献笔记三十五：水稻细胞器基因组数据做群体遗传学分析 文章中提到的水稻叶绿体的那篇论文中提供的 vcf格式文件，下载链接可以在论文中找到。论文中的vcf文件总共包括412个样本，本次分析只挑选出其中的20个，根据论文的补充材料1选取 Indica 10个， Japonica 10个。样本名分别为：

 

  
RWG-024
  
RWG-025
  
RWG-026
  
RWG-027
  
RWG-028
  
RWG-029
  
RWG-033
  
RWG-035
  
RWG-050
  
RWG-052
  
RWG-118
  
RWG-119
  
RWG-120
  
RWG-121
  
RWG-122
  
RWG-123
  
RWG-124
  
RWG-125
  
RWG-126
  
RWG-131
 

• 使用 vcftools工具从所有样本的vcf文件中选取所需要的样本，将样本名字放到文本文件里，我命名为 Japonica_Indica_sample_name.txt

 

  
vcftools --vcf 412_all_cp.recode.eva.vcf --keep Japonica_Indica_sample_name.txt --recode --recode-INFO-all --out keep_20_indv
 

• 查看新生成的文件的样本名

 

  
bcftools query -l keep_20_indv.recode.vcf
 

• 接下来的分析只保留snp位点

 

  
vcftools --vcf keep_20_indv.recode.vcf --remove-indels --recode --recode-INFO-all --out keep_20_indv.snp
 

• 接下来内容参考 The Simple Fool's Guide to population genomics using RNA-Seq

• http://sfg.stanford.edu/

• 首先是为smartpca程序准备输入文件 使用到的脚本是 vcf2smartpca.py；下载链接https://github.com/DeWitP/SFG/blob/master/scripts/vcf2smartpca.py 使用到的命令

 

  
python2 vcf2smartpca.py keep_20_indv.snp.recode.vcf passing_snps_q30 30 Japonica_Indica_sample_name_1.txt
 

第一个位置的参数是输出结果文件的前缀；第二个位置的参数是设置snp的最低质量值 第三个位置的参数准备一个输入文件，第一列是样本名，第二列是样本来自哪个种群，空格分隔

 

  
RWG-024 Japonica
  
RWG-025 Japonica
  
RWG-026 Japonica
  
RWG-027 Japonica
  
RWG-028 Japonica
  
RWG-029 Japonica
  
RWG-033 Japonica
  
RWG-035 Japonica
  
RWG-050 Japonica
  
RWG-052 Japonica
  
RWG-118 Indica
  
RWG-119 Indica
  
RWG-120 Indica
  
RWG-121 Indica
  
RWG-122 Indica
  
RWG-123 Indica
  
RWG-124 Indica
  
RWG-125 Indica
  
RWG-126 Indica
  
RWG-131 Indica
 

• 使用smartpca程序计算 这里参考利用EIGENSOFT中的smartpca模块进行PCA分析 直接使用 conda进行安装 conda install eigensoft 运行

 

  
smartpca -p par.passing_snps_q30 > passing_snps_q30_logfile.txt
 

本次分析会产生几个输出文件，我们感兴趣的是以.evec结尾的文件

 

  
 #eigvals: 15.125 1.810 1.011 0.415 0.342 0.155 0.096 0.030
  
 RWG-024 0.2357 0.2091 -0.3794 -0.0485 -0.0019 -0.1553 -0.2832 -0.3160 Japonica
  
 RWG-025 0.2181 0.0415 0.3163 0.2661 0.0070 0.0007 -0.1023 -0.7693 Japonica
  
 RWG-026 0.2391 0.2318 -0.4315 0.0067 -0.0004 0.0992 0.0451 0.1263 Japonica
  
 RWG-027 0.2122 0.0282 0.1958 -0.0676 0.0001 -0.5280 0.2445 0.1501 Japonica
  
 RWG-028 0.2122 0.0282 0.1958 -0.0676 0.0001 -0.5280 0.2445 0.1501 Japonica
  
 RWG-029 0.2420 -0.9137 -0.2317 0.0194 -0.0000 0.0479 -0.0133 -0.0024 Japonica
  
 RWG-033 0.2391 0.2318 -0.4315 0.0067 -0.0004 0.0992 0.0451 0.1263 Japonica
  
 RWG-035 0.1925 0.0415 0.3297 0.3343 -0.0002 0.0680 -0.6865 0.4615 Japonica
  
 RWG-050 0.2177 0.0357 0.3142 -0.8147 -0.0005 0.3631 -0.0824 -0.0127 Japonica
  
 RWG-052 0.2249 0.0869 0.2121 0.3764 0.0100 0.5096 0.5542 0.1152 Japonica
  
 RWG-118 -0.2290 -0.0026 -0.0141 -0.0069 0.9472 -0.0024 -0.0019 0.0010 Indica
  
 RWG-119 -0.2233 -0.0021 -0.0088 -0.0006 -0.1063 0.0035 0.0056 -0.0115 Indica
  
 RWG-120 -0.2182 -0.0017 -0.0056 0.0004 -0.1107 -0.0022 -0.0090 0.0616 Indica
  
 RWG-121 -0.2233 -0.0021 -0.0088 -0.0006 -0.1063 0.0035 0.0056 -0.0115 Indica
  
 RWG-122 -0.2233 -0.0021 -0.0088 -0.0006 -0.1063 0.0035 0.0056 -0.0115 Indica
  
 RWG-123 -0.2233 -0.0021 -0.0088 -0.0006 -0.1063 0.0035 0.0056 -0.0115 Indica
 

• 将这个文件导出，使用ggplot2做一个散点图 先删掉第一行，使用nodpad++打开，按住ALT键删除前面空白行。

 

  
df<-read.table("../../vcf_handling/passing_snps_q30_3645loci.evec",
  
 header=F)
  
dim(df)
  
head(df)
  
library(ggplot2)
  
install.packages("ggalt")
  
library(ggalt)
  
ggplot()+
  
 geom_point(data=df,aes(x=V2,y=V3,color=V10,shape=V10))+
  
 scale_color_manual(values=c("red","blue"))+
  
 theme_bw()+
  
 theme(legend.title = element_blank())+
  
 labs(x="PC1",y="PC2")+ylim(-0.1,0.3)+
  
 geom_encircle(data=df[c(1:5,7:10),],aes(x=V2,y=V3),
  
 color="blue",size=1,expand=0.05)
 

本文中用到的vcf文件大家可以到原论文中找到下载链接

或者直接关注我的公众号

小明的数据分析笔记本

后台回复 PCA_VCF即可