群体遗传中基于SNP的PCA分析

基于群体遗传中变异信息文件VCF来分析PCA

第一种方法

• 可以使用plink软件直接进行分析

plink --vcf all_genotypegvcf_filter_remove.vcf --pca -out all_genotypegvcf_plink

会有三个结果文件，
all_genotypegvcf_plink_plink.log：这个是日志文件
all_genotypegvcf_plink.eigenval：这个是样品名称的文件
all_genotypegvcf_plink.eigenvec：这个是pca结果

第二种方法

• 利用vcftools软件进行格式转换

vcftools --vcf all_genotypegvcf_filter_remove.vcf --plink --out all_genotypegvcf_vcftools #输出文件只要前缀就可以了

此时会生成两个文件
all_genotypegvcf_vcftools.ped
all_genotypegvcf_vcftools.map

再利用plink软件进行数据格式转换

plink --noweb --file all_genotypegvcf_vcftools --make-bed --out all_genotypegvcf_plink #输入和输出文件都只要前缀

生成3个文件
all_genotypegvcf_plink.bed
all_genotypegvcf_plink.bim
all_genotypegvcf_plink.fam

最后利用gcta软件进行pca构建

gcta --bfile all_genotypegvcf_plink --make-grm --autosome --out all_genotypegvcf_plink #输入和输出文件都只要前缀

此时生成一个文件
all_genotypegvcf_plink.grm.gz

gcta --grm all_genotypegvcf_plink --pca 8 --out all_genotypegvcf_plink

生成两个文件
all_genotypegvcf_plink.eigenval
all_genotypegvcf_plink.eigenvec

在all_genotypegvcf_plink.eigenvec文本的第一行再加入一行（之间以空格隔开）

1 2 pc1 pc2 pc3 pc4 pc5 pc6 pc7 pc8

打开R软件

#读取文件
a <- read.table("D:/all_genotypegvcf_plink.eigenvec", header=TRUE,sep=" ") 
#画图
plot(a$pc1,a$pc2, pch=c(1,2,3,4,5,6,7,8), col=c(1,2,3,4,5,6,7,8) , main="pca",xlab="pc1",ylab="pc2")
#添加图例
legend("bottomleft", c("A","B","C","D","E","F","G","H"), pch=c(1,2,3,4,5,6,7,8), col=c(1,2,3,4,5,6,7,8))

第三种方法

使用SNPRelate分析，但是这个我还没有尝试过