GWAS实战之制作PLINK格式的文件（上）
GWAS实战之制作PLINK格式的文件（上）
在上一期内容中，小陈让大家下载了一些数据。从今天开始，我们就要用这些数据干点酷酷的事了。
library(data.table)
library(GEOquery)
gset <- getGEO("GSE148812",getGPL=T) #使用GEOquery包获取数据矩阵
length(gset)
# 1
gset <- gset[[1]] # 提取矩阵数据
save(gset, file = "gset.Rdata") #保存
load("./gset.Rdata") # 加载
pdata <- pData(gset) # 提取样本的表型信息
pdata <- pdata[,c("title","age:ch1", "gender:ch1", "smoking_status:ch1")] #选取需要的列
head(pdata)
从上图中我们可以看到，”title”这一列是由两部分组成，一个是原始数据中样本的ID，一个是样本的表型（GSM开头的ID只是样本在GEO中的ID，不是原始数据里样本的ID），我们需要将原始的样本ID提取出来，代码如下：
id <-unlist(strsplit(as.character(pdata$title),split=":"))[seq(1,dim(pdata)[1]*2,2)]# 分割并提取原始样本ID
pdata$id <- id
head(pdata)
fwrite(pdata,"./sample_info.csv") #保存样本信息至此，我们已经把表型数据做好了。
在PLINK软件中，我们通常需要两个文件，一个是以.map为后缀的文件，另一个是以.ped为后缀的文件。其中，.map文件存储的是SNP位点的信息，主要由四列构成，第一列是染色体位置，第二个是SNP的ID（通常是rsID，但也可以是其它ID，只要保证唯一性即可），第三个是摩尔根位置，通常都可以设为0，第四个是碱基对位置（base position）。
而.ped文件存储了样本信息，包括表型和基因型，其列数在6列以上，前六列数据和.fam文件的前六列一致，在往期推文中可以查到------初探PLINK文件格式（bed，bim，fam）。从第七列开始就是基因型信息，基因型用A/T/C/G表示，如果基因型信息缺失，则用0表示。如下图所示，第一个人的第一个SNP的基因型信息就是缺失的（第七列和第八列为0），而其第二个SNP的基因型就是AA。
接下来我先简单介绍一下如何制作.map文件。
我们把之前下载的SOFT文件解压缩后打开（GWAS实战教程前言与示例数据），将前后的注释信息全部删除，最后的到如下所示的数据：
annot <-fread("GSE148812_family.soft", header=T) # 读取SOFT文件
mapping <-annot[,c("ID","Name")] # 将ID 和Name这两列选好用于后续的mapping
fwrite(mapping, 'mapping_exm.tsv',sep='\t') # 保存mapping文件
head(mapping)
这里可能有小伙伴会问，为什么要提取这两列呢？因为在这套数据里的ID这一列是作为突变的marker在基因型文件中使用的，它是用来和基因型文件匹配用的，而Name这一类又包含rsID的信息，是后续注释用的，因此我们需要把这两列提取出来。
colnames(annot) # 查看annot的列名并确定制作map文件所要提取的列
annot <- annot[,c("Chr","ID","AddressB_ID", "MapInfo")] # 提取map文件所需的4列信息
annot$AddressB_ID <- 0 #将第三列数据全部变成0
head(annot)
fwrite(annot, 'myWES.map', sep = '\t', col.names=F)# 保存数据并去掉列名，使用\t分割 关于表型文件和.map文件的制作就先讲到这里，下期我将介绍如何制作.ped文件，敬请期待！
腾讯云开发者
