加载中…[vcftools]基于全基因组snp数据如何进行主成分分析（PCA）
(2014-08-26 02:16:55)
标签：
ngsvcftools |
分类： 【转】NGS算法 |
http://www.bioask.net/question/238
经过“基因组-health
(213256700)”的热情指导，以及我本人的理解，现将如何基于全基因组的SNP数据进行PCA分析流程记录下来，不妥之处欢迎拍砖。
1）全基因组snp数据格式为 .vcf
2）利用vcftools软件进行格式转换：vcftools --vcf tmp.vcf --plink --out tmp
此时会生成两个文件：tmp.ped 和 tmp.map
3）利用plink软件进行数据格式转换：./plink --noweb --file tmp --make-bed --out tmp
注意，输入文件和输出文件都不需要文件名的后缀，此时生成3个文件：tmp.bed，tmp.bim 和 tmp.fam
4）利用gcta软件进行pca构建
4.1 ./gcta --bfile tmp --make-grm --autosome --out tmp
此时生成一个文件：tmp.grm.gz
4.2 ./gcta --grm tmp --pca 3 --out pcatmp
此时生成两个文件：pcatmp.eigenval 和 pcatmp.eigenvec
5）将生成的pcatmp.eigenvec用文本编辑器打开，在最上面加入一行：1 2 pc1 pc2 pc3（之间以空格隔开），保存
6）打开R软件
6.1 输入文件：a <- read.table("D:/pcatmp.eigenvec", header=TRUE)
6.2 绘散点图：plot(a$pc1,a$pc2, pch=c(1,2,3,4,5,6,7,8,9,10), col=c(1,2,3,4,5,6,7,8,9,10) , main="pca",xlab="pc1",ylab="pc2")
6.3 添加图例： legend("bottomleft", c("CL","IN","GZ","DA","PP","YN","DX","JY","NP","SL"), pch=c(1,2,3,4,5,6,7,8,9,10), col=c(1,2,3,4,5,6,7,8,9,10))
文件 > 另存为 > Jpeg or Tiff
That's all, Game over. 再次向基因组-health (213256700)予以致谢！
1）全基因组snp数据格式为 .vcf
2）利用vcftools软件进行格式转换：vcftools --vcf tmp.vcf --plink --out tmp
此时会生成两个文件：tmp.ped 和 tmp.map
3）利用plink软件进行数据格式转换：./plink --noweb --file tmp --make-bed --out tmp
注意，输入文件和输出文件都不需要文件名的后缀，此时生成3个文件：tmp.bed，tmp.bim 和 tmp.fam
4）利用gcta软件进行pca构建
4.1 ./gcta --bfile tmp --make-grm --autosome --out tmp
此时生成一个文件：tmp.grm.gz
4.2 ./gcta --grm tmp --pca 3 --out pcatmp
此时生成两个文件：pcatmp.eigenval 和 pcatmp.eigenvec
5）将生成的pcatmp.eigenvec用文本编辑器打开，在最上面加入一行：1 2 pc1 pc2 pc3（之间以空格隔开），保存
6）打开R软件
6.1 输入文件：a <- read.table("D:/pcatmp.eigenvec", header=TRUE)
6.2 绘散点图：plot(a$pc1,a$pc2, pch=c(1,2,3,4,5,6,7,8,9,10), col=c(1,2,3,4,5,6,7,8,9,10) , main="pca",xlab="pc1",ylab="pc2")
6.3 添加图例： legend("bottomleft", c("CL","IN","GZ","DA","PP","YN","DX","JY","NP","SL"), pch=c(1,2,3,4,5,6,7,8,9,10), col=c(1,2,3,4,5,6,7,8,9,10))
文件 > 另存为 > Jpeg or Tiff
That's all, Game over. 再次向基因组-health (213256700)予以致谢！
喜欢
0
赠金笔