根据芯片号拆分plink文件的代码实现

数据分析 / python
字数统计: 2.8k   |   阅读时长≈ 11 分钟

这个代码同样是一个处理基因型数据的入门级别的 python 脚本。

目的是根据芯片号，将一对plink文件，提取相应的基因型个体，拆分为一对或多对子文件。

使用软件

  • python 3.8 及以上版本

输入文件格式

plink格式的基因型文件包含两个文件：基因型文件（以.ped结尾）和图谱文件（以.map结尾），格式如下：

（1）基因型文件（ped文件）

FAM001 IND001 0 0 1 2 A A G G A C …

FAM001 IND002 0 0 1 2 A A A G 0 0 …

……

格式说明：前6列分别是：家系编号，个体编号，父亲编号，母亲编号，性别，表型值；性别编码规则：‘1’代表雄性，‘2’代表雌性，其他字符代表性别未知；从第7列开始是基因型（只支持二等位位点），每个位点2列，默认‘0’代表基因型缺失。

（2）图谱文件（MAP文件）

1 rs234567 0 1237793

1 rs224534 0 1237697

1 rs233556 0 1337456

……

格式说明：包含4列：染色体编号，标记名称，标记遗传图谱，标记物理位置；对于性染色体和线粒体的‘染色体’编号规则：总的常染色体对数加一 表示X染色体，Y染色体、假常染色体和线粒体在此基础上依次加一。

更加详细说明见官网

拆分文件

一列或两列，无标题，以空格或制表符分隔。

如果是一列，则为一列芯片号，程序会将这一列芯片号的基因型提取出来，生成一对子文件，命名为 split_plink.ped 和 split_plink.map

如果是两列，第一列是芯片号，第二列是芯片分类名（例如品种）。此时程序会将不同分类的芯片提取出来，各生成一对plink文件，plink文件前缀即为分类名。

输出文件格式

一对或多对 plink 文件。

运行代码

代码文件见：split_plink.py

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

1
python split_plink.py plink_prefix split_plink.txt

参数说明：

plink_prefix : 第一个参数为总的plink前缀

split_plink.txt：第二个参数为拆分文件

例子

举个简单的例子，方便理解。

总的plink为 1.ped 1.map。

1.ped 内容如下：

1
2
3
4
5
221 30800300 0 0 0 -9 C C T T T T A A C C
222 30800299 0 0 0 -9 T C T T T T G A T C
223 30800298 0 0 0 -9 T C T T T T A A C C
224 30800297 0 0 0 -9 C C T T T T G A T C
225 30800296 0 0 0 -9 T C T T T T A A C C

1.map 内容如下：

1
2
3
4
5
1	1_242598	0	242598
1	1_10673082	16.99376	10673082
1	1_10723065	17.17102	10723065
1	1_11407894	19.59961	11407894
1	1_11426075	19.66408	11426075

第一个拆分文件为 split_plink1.py ，内容只有一列，为3个芯片号。

1
2
3
4
30800300
30800299
30800298

运行命令，提取这三个个体的基因型，第一个参数为 1，为总的plink文件前缀；第二个参数为拆分文件名。

1
python split_plink.py 1 split_plink1.txt

生成 split_plink.ped 和 split_plink.map ，生成的 split_plink.map文件与1.map 一致， split_plink.ped 文件内容如下，为需要的三个芯片号的基因型。

1
2
3
4
221 30800300 0 0 0 -9 C C T T T T A A C C
222 30800299 0 0 0 -9 T C T T T T G A T C
223 30800298 0 0 0 -9 T C T T T T A A C C

第二个拆分文件为 split_plink2.py ，内容为两列，第二列模拟为 3 个物种，分别命名为 DD, LL, YY

1
2
3
4
5
6
30800300	DD
30800299	DD
30800298	LL
30800297	YY
30800296	YY

运行命令，提取这三个个体的基因型，第一个参数为 1，为总的plink文件前缀；第二个参数为拆分文件名。

1
python split_plink.py 1 split_plink2.txt

生成了三对 plink 文件，DD.ped/map ，LL.ped/map , YY.ped/map ，三个map 文件均与 1.map 一致。三个 ped 文件内容如下：

DD.ped 内容为

1
2
3
221 30800300 0 0 0 -9 C C T T T T A A C C
222 30800299 0 0 0 -9 T C T T T T G A T C

LL.ped 内容为

1
2
223 30800298 0 0 0 -9 T C T T T T A A C C

YY.ped 内容为

1
2
3
224 30800297 0 0 0 -9 C C T T T T G A T C
225 30800296 0 0 0 -9 T C T T T T A A C C

实现了将不同分类拆分为不同plink文件的需求。

代码说明

注释说明

首先是两行注释，第一句以#!开头，后面接着的是我在linux系统中安装的python路径。

1
#!/mnt/data/zhouziwen/bin/miniConda/miniconda3/bin/python

添加这句话的目的是指定解释器，就是说用什么去执行这个脚本。添加这个命令后，在给这个程序添加执行权限后（如何添加执行权限见linux 的 chmod命令 ），可以将这个脚本视为正常的shell脚本执行，例子如下（需要先修改为你自己的python路径）

1
./pick0.py id.txt bigfile.txt out.txt

如果你没看懂，不要紧，这个不重要，你就还是按照上面的运行代码部分写的运行命令去运行即可。

第二行注释是指定字符编码为 utf-8（字符编码的内容可以看看这篇博客，字符编码那点事） ，不然写中文注释会报错。

1
#coding=utf-8

后面是我自己写的两行注释，说明代码功能

1
2
# 按照一个芯片号文件拆分plink文件
# 如果芯片号文件只有一列，默认输出文件为 split_plink.ped split_plink.map

导入模块和参数

这里只导入了三个内置模块。

导入 defaultdict 函数的作用是可以指定字典的值为一种默认的数据类型，比如列表，下面会用到。

1
2
import os,sys
from collections import defaultdict

接下来，导入参数并读取文件。这里将第一个参数 (sys.argv[1]) 命名为 plink_prefix ，即 plink 文件前缀。第二个参数为拆分文件，打开文件，命名为 split_file 。打开 ped 文件，命名为 plink_file 。 最后一句将map文件名命名为 map_file_name

1
2
3
4
5
plink_prefix = sys.argv[1] # plink文件前缀
split_file = open(sys.argv[2],'r') # 拆分信息文件，无标题，一列（芯片号）或两列（芯片号-分类名）
plink_file = open(f"{plink_prefix}.ped",'r')
map_file_name = f"{plink_prefix}.map"

获取拆分文件列数

这里是通过 readline() 读取拆分文件的第一行内容，然后使用 split() 函数对第一行内容按空白字符拆分为列表，最后用 len() 函数求列表元素的个数。通过这种方法，实际获取的是拆分文件第一行的列数，命名为 column_num。

1
column_num = len(split_file.readline().split())

之后记得通过 seek(0) 函数将指针复原，不然下面对文件进行循环就是从第二行开始了，相当于剔除了第一行。使用了 seek(0) ，下面对文件进行循环，就依然是从第一行开始。

1
split_file.seek(0)

读取拆分文件内容

先创建一个空列表和 一个值默认为列表的字典

注意，使用 defaultdict(list) 必须先从collections 包中导入defaultdict 函数。

1
2
set1 = set()
dick = defaultdict(list)

遍历拆分文件内容，形成 dick 字典的内容。这里分了两种情况，拆分文件为一列和两列的情况，但是内容大同小异。

此处代码对于拆分文件的每一行首先判断列数是否均符合要求（是不是都是一列或两列），存在不符合要求的行，则报错 (sys.exit(1)) ；之后判断不同行之间的芯片号是否出现重复（通过 set1 进行判断），如果重复，则报错；最后将文件内容放入 dick 字典中，字典的键为分类名（如果拆分文件只有一列，则所有芯片号的分类名均设置为 split_plink），值为属于该分类名的所有芯片号组成的列表。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
j=0
if column_num == 1: # 只有一列
    for i in split_file:
        j += 1
        f = i.split()
        if len(f) == 1:
            if f[0] not in set1:
                set1.add(f[0])
                dick["split_plink"].append(f[0]) # 品种：芯片号的列表
            else:
                print(f"Error: duplicated chipid {f[0]} in {sys.argv[1]}\n")
                sys.exit(1)
        else:
            print(f"Error:{j} row in {sys.argv[1]} is not 1 columns\n")
            sys.exit(1)
elif column_num == 2:
    for i in split_file:
        j += 1
        f = i.split()
        if len(f) == 2:
            if f[0] not in set1:
                set1.add(f[0])
                dick[f[1]].append(f[0]) # 品种：芯片号的列表
            else:
                print(f"Error: duplicated chipid {f[0]} in {sys.argv[1]}\n")
                sys.exit(1)
        else:
            print(f"Error:{j} row in {sys.argv[1]} is not 2 columns\n")
            sys.exit(1)
else:
    print(f"Error: first row in {sys.argv[1]} is not 1 or 2 columns\n")
    sys.exit(1)

读取 ped 文件内容

先创建一个空集合和空字典

1
2
set2=set()
dick2 = {}

遍历 ped 文件，只做是否存在重复芯片号的判断，如果存在则报错，程序中止。

生成键为芯片号，值为整行内容的字典 dick2

1
2
3
4
5
6
7
8
9
for i in plink_file:
    f = i.split()
    if f[1] not in set2:
        set2.add(f[1])
        dick2[f[1]] = i # 芯片号：芯片信息的行
    else:
        print(f"Error: duplicated chipid {f[1]} in {sys.argv[2]}.ped\n")
        sys.exit(1)

写入结果文件

这里我们遍历 dick 的所有键，即遍历所有的分类名。首先这里先将原始的map文件复制为当前的分类名的map文件；之后我们创建out_plink_list 这个列表，作为拆分的ped文件内容的临时存储的地方；然后遍历该分类名下的所有芯片号，将这些芯片号的基因型信息存放在out_plink_list 列表中（这里会判断这些芯片号是否均在 plink 文件中，如果不在则程序报错中）。最后创建以分类号为前缀的 ped 文件，将 out_plink_list 中的全部内容写入其中。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
for i in dick: # 所有品种
    # 复制map文件
    os.system(f"cp {map_file_name} {i}.map") 
    out_plink_list = []
    for i2 in dick[i]:
        # 如果芯片号不在plink文件中报错
        if i2 not in dick2:
            print(f"Error: {i2} in {sys.argv[1]} not in {sys.argv[2]}\n")
            sys.exit(1)
        out_plink_list.append(dick2[i2])
    # 生成plink 文件
    out_plink_file = open(f"{i}.ped",'w')
    out_plink_file.writelines(out_plink_list)
    out_plink_file.close()

这里，如果拆分文件只有一列，那么 dick 中只有一种分类号，即split_plink ，因此最终只会生成一对plink文件，即 split_plink.ped/map 。

关闭文件

所有打开的文件记得关闭，由于写入的文件已经在循环中关闭了，最后只需要关闭读取的两个文件。

1
2
split_file.close()
plink_file.close()

小结

这个脚本其实和上一篇根据ID提取文件中相应行的python脚本实现与详细讲解 有异曲同工之处，说白了都是根据某一列从提取需要的特定行，像这种需求用 python 很容易实现，而且可以再加一些判断，可以做到比较稳健。

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

请我喝杯茶吧~

支付宝
微信