规整某类基因型数据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

1
Rscript mass_to_plink.R rawdata.xlsx raw.map

参数说明

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

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

运行python代码

python 代码文件见:mass_to_plink.py

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

1
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)等等。

1
2
3
4
5
6
7
8
$ sha1sum new*.ped
5150663442126f2f4e5b57e9d2d58118bf1f90f3 new1.ped
5150663442126f2f4e5b57e9d2d58118bf1f90f3 new2.ped

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

  • 版权声明: 本博客所有文章除特别声明外,著作权归作者所有。转载请注明出处!
  • Copyrights © 2019-2022 Vincere Zhou
  • 访问人数: | 浏览次数:

请我喝杯茶吧~

支付宝
微信