全基因组关联分析（GWAS）大家都不陌生，今天我们给大家介绍下各种格式之间转化在R语言是怎么实现的。首先我们来看下GWAS都有哪些数据格式：
1. HapMap格式，这也是当时全基因组计划的简称，自此这个也成为了其主要的一种文件格式。
其变量构成：
| 名称 | 描述 |
|---|---|
| rs# | SNP的标识符 |
| alleles | 基于NCBI数据库的SNP等位基因 |
| chrom | SNP所在的染色体 |
| pos | SNP在染色体上的位置 |
| strand | 相对于参考序列的方向，（+）前，（-）反 |
| assembly# | NCBI参考序列的版本 |
| center | 基因分型 |
| protLSID | HapMap protocol |
| assayLSID | HapMap assay used for genotyping |
| panelLSID | panel of individuals genotyped |
| QCcode | 质量控制 |
| …… | 样本数据 |
image.png
数据实例：
[图片上传中...(image.png-77c038-1608883012663-0)]
2. plink数据 ped/map。这个数据格式需要两个文件共同保存数据一个map文件一个ped数据文件。
Ped文件的结构：
image.png
image.png
Map文件结构：
image.png
image.png
3. BED/BIM/FAM文件
BED文件结构主要是二进制文件，它的具体内容我们估计不好看，就以网页的数据为例，给大家看下长啥样子：
image.png
BIM文件结构：
image.png
image.png
FAM文件的结构：
image.png
以上就是GWAS主要的文件结构，在R语言中还有另外一个结构就是GDS结构，此结构由R包gdsfmt进行创建编辑。今天我们主要讲下在包SNPRelate中如何实现这些数据结构之间的转化。
首先看下包的安装，还是需要bioconductor环境进行启动安装：
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install(version = "3.12")
BiocManager::install(c("SNPRelate"))
library("SNPRelate")
vcf.fn<-"zhili2.gvcf"
You can use VCFtools to make a PED and MAP file from VCF. This is PLINK format. Many PCA programs take PLINK input or offer conversion scripts.
I ended up using SNPRelate. After some silly errors here is how I got it to work:
setwd("/xxx/pca")
library("SNPRelate")
vcf.fn<-"~/xxx/tmp.vcf"
snpgdsVCF2GDS(vcf.fn, "ccm.gds", method="biallelic.only")
genofile <- openfn.gds("ccm.gds")
ccm_pca<-snpgdsPCA(genofile)
plot(ccm_pca$eigenvect[,1],ccm_pca$eigenvect[,2] ,col=as.numeric(substr(ccm_pca$sample, 1,3) == 'CCM')+3, pch=2)
接下来看下里面数据转化的主要函数：
| 名称 | 功能 |
|---|---|
| snpgdsPED2GDS | 将ped/map文件转化为GDS |
| snpgdsBED2GDS | 将BED/BIM/FAM文件转化为GDS |
| snpgdsGDS2PED | 将GDS文件转化为PED/MAP文件 |
| snpgdsGDS2BED | GDS转化为BED/BIM/FAM文件 |
| snpgdsVCF2GDS | VCF文件转化为GDS文件 |
