Excel格式的SNP数据怎么变为plink格式

2022-05-19 09:25:53 浏览数 (1)

有时候,我们会遇到Excel格式的基因型数据,这篇博文介绍一下如何手动转为plink格式。

可以在Excel中整理,也可以在R语言中整理。数据量少的话,就在Excel中整理,数据量大的话,就在R语言中整理就行。

主要思路是根据plink的格式特点,针对性的满足,然后导出,就可以了。

1. Excel中的基因型数据格式

第一列是snpID,第二列是染色体,第三列是物理位置,第四列是参考基因组分型,第五列以后是每个样本的具体分型。

整体而言,每一行是一个snp,第五列以后每一列是一个样本。

2. plink的格式

「.map格式」

格式说明链接: http://zzz.bwh.harvard.edu/plink/data.shtml#map

❝map格式的文件, 主要是图谱文件信息, 主要包括染色体名称, 所在的染色体和所在染色体的坐标.❞

1, map文件没有行头

2, map文件包括四列: 染色体, SNP名称, SNP位置, 碱基对坐标

  • 染色体编号为数字, 未知为0
  • SNP名称为字符或数字, 如果不重要, 可以从1编号, 注意要和bed文件SNP列一一对应
  • 染色体的摩尔位置(可选项, 可以用0)
  • SNP物理坐标

3, 如果只有SNP名称, 可以手动构建map文件, 第二列为SNP名称, 其它三列为0即可.

Example:

代码语言:javascript复制
1 snp1 0 1
1 snp2 0 2
1 snp3 0 3
  • 这里有3个SNP, 分别名为snp1, snp3, snp3 「(第二列)」
  • 这三个SNP在第一个染色体上 「(第一列)」
  • 第三列为0
  • 第四列为SNP所在染色体的坐标

「.ped格式」格式说明链接:http://zzz.bwh.harvard.edu/plink/data.shtml#ped

❝bed格式的文件, 主要包括SNP的信息, 包括个体ID, 系谱信息, 表型和SNP的分型信息.❞

1, 数据没有行头, 空格或者tab隔开的文件 2, 必须要有六列, 包括系谱信息, 表型信息

  • 第一列: Family ID # 如果没有, 可以用个体ID代替
  • 第二列: Individual ID # 个体ID编号
  • 第三列: Paternal ID # 父本编号
  • 第四列: Maternal ID # 母本编号
  • 第五列: Sex (1=male; 2=female; other=unknown) # 性别, 如果未知, 用0表示
  • 第六列: Phenotype # 表型数据, 如果未知, 用0表示
  • 第七列以后: 为SNP分型数据, 可以是AT CG或11 12, 或者A T C G或1 1 2 2

3, 上面六列, 必须要有, 如果没有相关数据, 用0表示.

所以,下面的任务就是把Excel的格式,变为plink的ped和map格式。

3. R语言操作

3.1 读取数据

代码语言:javascript复制
library(openxlsx)
library(tidyverse)

dat = read.xlsx("SNP-excel.xlsx")
dat[1:10,1:10]

3.2 map数据整理

代码语言:javascript复制
map = dat %>% select(2,1,x = 3,p = 3)
head(map)

3.3 ped数据整理

下面这个代码复杂一点,主要的逻辑:

  • 去除中间的及列
  • 然后进行转置
  • 变为plink的格式
代码语言:javascript复制
ped = dat %>% select(-c(1:4)) %>% t() %>% as.data.frame() %>% 
  mutate(ID = rownames(ped)) %>% 
  mutate(x3=0,x4=0,x5=0,x6=0) %>% 
  select(FID=ID,IID=ID,x3,x4,x5,x6,everything())
ped[1:10,1:10]

3.4 导出数据

这里有个坑:默认plink的缺失为0,如果AT连写的格式,是00,但是如果一列都是00的时候,就变为了一个0,就会map和ped不匹配。

比较靠谱的方式是,将缺失变为##,然后将其变为00.

代码语言:javascript复制
library(data.table)
fwrite(map, "file.map",col.names = F,quote = F,sep = " ")
fwrite(ped, "file.ped",col.names = F,quote = F,sep = " ",na = "##")
代码语言:javascript复制
 sed -i 's/##/00/g' file.ped

4. 测试plink

代码语言:javascript复制
 plink --file file --missing

搞定!


大家好,我是邓飞,一个持续分享的农业数据分析师

0 人点赞