使用smartpca做种群遗传学中常用的PCA分析实例
使用smartpca做种群遗传学中常用的PCA分析实例
用户7010445
发布于 2020-03-03 14:54:01
发布于 2020-03-03 14:54:01
本文中使用到的数据是 文献笔记三十五：水稻细胞器基因组数据做群体遗传学分析 文章中提到的水稻叶绿体的那篇论文中提供的 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文件大家可以到原论文中找到下载链接
本文参与 腾讯云自媒体同步曝光计划，分享自微信公众号。
原始发表：2019-08-28，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读