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

这个代码同样是一个处理基因型数据的入门级别的 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.pedsplit_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.pedsplit_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-2022 Vincere Zhou
  • 访问人数: | 浏览次数:

请我喝杯茶吧~

支付宝
微信