规整某类基因型数据R代码与python代码实现

这个是规整某类格式”不规范“的基因型数据的脚本实现说明，同时采用了 R 和 python 进行实现。

这个需求实现起来说难不难，说易不易，仅是为了将下面我要提到的这类下机的基因型数据转变为常用的plink格式。这里的两个代码可能本身对其他人无法直接使用，但是还是有些借鉴作用。

使用软件

• python 3.8 及以上版本
• plink 软件 v1.90，并且提前加到环境变量中。

输入文件格式

下机基因型数据

示例数据见 rawdata.xlsx ，是一个 Excel 表，内容见下表。第一列是样本名称，后面每一列表示一个SNP，标题是SNP的名称。也就是说，一行表示一个样本，一列表示一个SNP。如果基因型为两个字符，说明是杂合子，例如GA；如果基因型是一个字符，说明是纯合子，例如T (正常应该记为TT，这里省略了一个字符)；如果单元格内容为空，说明此处基因型缺失。

id	seq-rs319899733	seq-rs320947633	seq-rs326926754	seq-rs327695512
24	T	A		CT
31	T	GA	C	C
23	CT	GA	C	C
29	T	GA		T

参考map文件

包含下机基因型所有SNP的 map 文件，根据这个文件来生成下机基因型数据的map文件，示例数据见 raw.map

map 文件格式说明见根据芯片号拆分plink文件的代码实现

输出文件格式

一对 plink 文件，文件格式说明见根据芯片号拆分plink文件的代码实现

运行代码

运行R代码

R 代码文件见：mass_to_plink.R

将两个输入文件和本程序放在同一文件夹下，运行命令示范如下。运行结束后，会生成规整后的 plink 文件 new.ped/map

Rscript mass_to_plink.R rawdata.xlsx raw.map

参数说明：

rawdata.xlsx: 第一个参数为下机基因型数据

raw.map：第二个参数为参考map文件。

运行python代码

python 代码文件见：mass_to_plink.py

将两个输入文件和本程序放在同一文件夹下，运行命令示范如下

python mass_to_plink.py rawdata.xlsx raw.map new2

参数说明：

rawdata.xlsx: 第一个参数为下机基因型数据

raw.map：第二个参数为参考map文件。

new2 : 输出plink文件前缀

运行结束后，生成以第三个参数为前缀的plink 文件。

示例

示例数据和代码我都放到了一个文件夹内 mass_to_plink ，有兴趣可以自己跑一下。

代码说明

这两个代码内部我都有比较详细的注释，这里就不再贴代码了。

两个代码的实现思路是一样的，首先要对基因型数据进行规整，缺失要替换为00，纯合子要改为两个字符，生成 compound 格式的 ped 文件（正常的ped文件一个SNP占两列，compound 格式则是一个SNP占一列，一个SNP的两个碱基之间无分隔符）。之后从参考map中提取相应的SNP的行，生成需要的map文件。最后通过 plink 软件将 compound 格式的文件转化为正常的plink格式。

后续

R 本身就有 dataframe 这种数据类型，而 python 是需要借助 pandas 模块进行处理。其实我第一次是用 R 实现的，后面又用 python 写了一遍。分别用 R 和 python 跑了一遍以后，我又用之前提到的比对基因型的代码，比较了两个生成的plink文件，基因型完全一致（具体参考 两个Plink文件样本间基因型比对代码实现 和 两个Plink文件同一SNP基因型比对代码实现）。

这里也有一个简单的办法，就是通过 sha1sum 或 md5sum 进行文件校验，你会发现两个 ped 文件 或 map 文件都是一样的。但是这种做法并不稳健，如果校验码不一致，不一定说明两个plink文件基因型不一致，还有很多原因导致校验码不同，比如 SNP 顺序不同，性染色体编码不同，杂合子显示顺序不同（如 A G 与 G A）等等。

$ sha1sum new*.ped
5150663442126f2f4e5b57e9d2d58118bf1f90f3  new1.ped
5150663442126f2f4e5b57e9d2d58118bf1f90f3  new2.ped

$ sha1sum new*.map
f899867a1d9f762b7a3d7fb231d5fc13fdc44aa4  new1.map
f899867a1d9f762b7a3d7fb231d5fc13fdc44aa4  new2.map