基于两个plink文件共同样本共同位点的基因型比对代码实现

比对两个plink文件共同样本共同位点的基因型,合并了 两个Plink文件样本间基因型比对代码实现两个Plink文件同一SNP基因型比对代码实现 的内容。

前言

整理文档的时候突然想到,既然比对两个plink文件时从共同样本角度和从共同SNP角度比对都需要进行,那为什么不合并成一个脚本呢?

是的,这里需要合并成一个脚本,因此就有了这个脚本。这个脚本是在compair_snps.py基础上新增了一些代码得到的,同时从共同样本和共同SNP角度去比较基因型。

使用软件

  • python 3.8 及以上版本
  • plink 软件 v1.90 版本

输入文件格式

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

输出文件说明

两个输出文件名称前面相同,后面分别为 _snp.txt_sample.txt ,分别为从SNP角度和样本角度的比对结果。

*_snp.txt : 第一列为SNP名称,第二列为两个plink文件中共同样本中不一致的样本数目,第三列为共同样本总数,第四列为前两列相除得到的不一致率。

示例如下:

8_118011832 3508 3529 0.994
5_39432167 3459 3524 0.9816
14_133523132 3428 3518 0.9744

*_sample.txt : 第一列为样本名称,第二列为两个plink文件中该样本的共同位点中不一致的样本数目,第三列为共同位点总数,第四列为前两列相除得到的不一致率。

202888800026_R03C02 0 51315 0
202888800026_R03C03 0 51315 0
202888800026_R03C04 0 51315 0

运行代码

代码文件见:compair_two_plinks.py

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

1
python compair_two_plinks.py --plink1 plink1 --plink2 plink2

可用参数说明

--plink1: 第一个plink文件前缀,必选。

--plink2: 第二个plink文件前缀,可选。

--exclude-miss : 基因型比对是否先剔除缺失值,可选,只能为yes和no,默认为yes。如果该选项设置为no,比对时不会剔除缺失样本对,如果存在同一个样本,在一个基因型文件中有基因型,另一个基因型文件中则为缺失,该样本同样会被视为不一致位点;如果该选项设置为 yes ,比对时会剔除至少一个样本为缺失的位点,只比对两个样本均不为缺失的位点。

--save-temp : 是否保留中间文件,可选,只能为yes和no,默认为 no。

--out-prefix : 结果文件前缀,可选,默认为 cmp_two_plinks

代码说明

由于这个脚本是基于 compair_snps.py 的基础修改的,前面的内容基本没变,只看后面修改的地方。

计算不一致率

这里新增了从样本角度比对的代码

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
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
#计算不一致率
out_sample_list = []
out_snp_list = []

if args.exclude_miss == 'yes':
for i in range(snp_num):
first_index = range(0,common_sample_num)
second_index = range(common_sample_num,total_sample_num)
v1 = array1[first_index,i]
v2 = array1[second_index,i]
non_miss_index = np.where((v1!=3) & (v2!=3))
total_num = non_miss_index[0].shape[0]
v1_new = v1[non_miss_index]
v2_new = v2[non_miss_index]
no_match_num = np.count_nonzero(v1_new!=v2_new)
if total_num > 0:
no_match_rate = no_match_num/total_num
else:
no_match_rate = 0
out_snp_list.append([snp_list[i],f"{no_match_num}",f"{total_num}",f"{no_match_rate:.4f}"])
for i in range(len(common_sample_list)):
first_index = i
second_index = i+common_sample_num
v1 = array1[first_index,:]
v2 = array1[second_index,:]
non_miss_index = np.where((v1!=3) & (v2!=3))
total_num = non_miss_index[0].shape[0]
v1_new = v1[non_miss_index]
v2_new = v2[non_miss_index]
no_match_num = np.count_nonzero(v1_new!=v2_new)
if total_num > 0:
no_match_rate = no_match_num/total_num
else:
no_match_rate = 0
out_sample_list.append([common_sample_list[i],f"{no_match_num}",f"{total_num}",f"{no_match_rate:.4f}"])

elif args.exclude_miss == 'no':
for i in range(snp_num):
first_index = range(0,common_sample_num)
second_index = range(common_sample_num,total_sample_num)
v1 = array1[first_index,i]
v2 = array1[second_index,i]
total_num = common_sample_num
no_match_num = np.count_nonzero(v1!=v2)
no_match_rate = no_match_num/total_num
out_snp_list.append([snp_list[i],f"{no_match_num}",f"{total_num}",f"{no_match_rate:.4f}"])
for i in range(len(common_sample_list)):
first_index = i
second_index = i+common_sample_num
v1 = array1[first_index,:]
v2 = array1[second_index,:]
total_num = snp_num
no_match_num = np.count_nonzero(v1!=v2)
no_match_rate = no_match_num/total_num
out_sample_list.append([common_sample_list[i],f"{no_match_num}",f"{total_num}",f"{no_match_rate:.4f}"])

将比对结果按从不一致率从大到小

1
2
sorted_out_snp_list = sorted(out_snp_list, key = lambda s:float(s[3]), reverse = True) #不一致率从大到小
sorted_out_sample_list = sorted(out_sample_list, key = lambda s:float(s[3]), reverse = True) #不一致率从大到小

写入结果文件,统计不一致的SNP数目和样本数目。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
unconsistent_snp_num = 0 # 不一致的SNP数目
unconsistent_sample_num = 0 # 不一致的样本数目

for i in sorted_out_snp_list:
out_snp_file.write('\t'.join(i)+'\n')
if int(i[1]) > 0:
unconsistent_snp_num += 1

for i in sorted_out_sample_list:
out_sample_file.write('\t'.join(i)+'\n')
if int(i[1]) > 0:
unconsistent_sample_num += 1

print(f"unconsistent sample number: {unconsistent_sample_num}\n")
print(f"unconsistent snp number: {unconsistent_snp_num}\n")

闲话

定期整理文档,梳理之前写的代码是有益的。

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

请我喝杯茶吧~

支付宝
微信