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

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

前言

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

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

使用软件

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

输入文件格式

plink 文本文件

格式说明见根据芯片号拆分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

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

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 的基础修改的，前面的内容基本没变，只看后面修改的地方。

计算不一致率

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

#计算不一致率
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}"])

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

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数目和样本数目。

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")

闲话

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