如何为药物遗传学分析挑选匹配度最优的单倍型？（Python 2.7）

嘿，我之前做过类似的遗传分析工具开发，给你梳理下思路和具体的Python 2.7实现方案，应该能帮到你！

要解决这个问题，核心是对齐患者SNP位点与单倍型的SNP位点，然后计算每个单倍型的匹配度，最终选出匹配度最高的那个。具体分这几步：

1. 标准化单倍型数据 #

首先得把你从数据库拿到的原始单倍型数据整理成易处理的结构。比如你给出的示例：

['UGT1A116', 'T;A;T;C;G;A;G;C'] ['UGT1A117', 'T;A;T;C>T;G;A;G;C']

这里注意C>T这种格式，通常表示该位点的等位基因是T（相对于参考序列的C），我们需要把每个单倍型的等位基因列表标准化，统一提取实际的等位基因碱基（比如把C>T转成T），同时按SNP的位置排序（确保和患者的SNP位点顺序完全对应）。

2. 对齐患者的SNP数据 #

你需要确保患者的SNP数据是按单倍型定义中的SNP位置顺序排列的，每个位点对应一个等位基因（比如患者在某个染色体位置的碱基是T，就对应单倍型中同位置的条目）。如果患者数据是按位置存储的，就先按染色体、起始位排序，和单倍型的SNP列表顺序对齐。

3. 计算匹配度并筛选 #

对每个单倍型，统计它和患者SNP列表中匹配的位点数量，或者计算匹配率（匹配数/总位点数），然后找出匹配度最高的单倍型。如果有多个单倍型匹配度相同且最高，建议全部返回（避免漏检）。

下面是具体的代码示例，包含数据解析、匹配计算的完整流程：

# 1. 解析并标准化单倍型数据
def parse_haplotypes(raw_haplotype_data):
    haplotype_dict = {}
    for item in raw_haplotype_data:
        name, alleles_str = item
        # 处理等位基因字符串，拆分并标准化（比如把C>T转成T）
        alleles = []
        for allele in alleles_str.split(';'):
            # 处理变异格式，提取实际等位基因
            if '>' in allele:
                standardized = allele.split('>')[1]
            else:
                standardized = allele
            alleles.append(standardized)
        haplotype_dict[name] = alleles
    return haplotype_dict

# 2. 计算匹配度
def calculate_match_score(haplotype_alleles, patient_alleles):
    # 确保两个列表长度一致（位点数量相同）
    if len(haplotype_alleles) != len(patient_alleles):
        raise ValueError("单倍型和患者的SNP位点数量不匹配！")
    # 统计匹配的位点数量
    match_count = 0
    for h_allele, p_allele in zip(haplotype_alleles, patient_alleles):
        if h_allele == p_allele:
            match_count += 1
    # 返回匹配数和匹配率
    return match_count, float(match_count)/len(haplotype_alleles)

# 3. 筛选匹配度最高的单倍型
def find_best_match(haplotype_dict, patient_alleles):
    best_score = -1
    best_haplotypes = []
    for name, alleles in haplotype_dict.items():
        match_count, match_rate = calculate_match_score(alleles, patient_alleles)
        # 更新最佳匹配
        if match_count > best_score:
            best_score = match_count
            best_haplotypes = [(name, match_count, match_rate)]
        elif match_count == best_score:
            best_haplotypes.append((name, match_count, match_rate))
    return best_haplotypes

# 示例使用
if __name__ == "__main__":
    # 模拟从数据库获取的原始单倍型数据
    raw_data = [
        ['UGT1A1*16', 'T;A;T;C;G;A;G;C'],
        ['UGT1A1*17', 'T;A;T;C>T;G;A;G;C'],
        ['UGT1A1*28', 'T;A;T;C;G;A;G;C>T']
    ]
    # 解析单倍型
    haplotype_dict = parse_haplotypes(raw_data)
    # 模拟患者的SNP数据（按单倍型位点顺序排列）
    patient_alleles = ['T', 'A', 'T', 'T', 'G', 'A', 'G', 'C']
    # 找出最佳匹配
    best_matches = find_best_match(haplotype_dict, patient_alleles)
    
    # 输出结果
    print("匹配度最高的单倍型：")
    for name, count, rate in best_matches:
        print("%s - 匹配位点：%d，匹配率：%.2f%%" % (name, count, rate*100))

• 位点对齐是关键：一定要保证单倍型的SNP列表和患者的SNP列表是按相同的位置顺序排列的，否则匹配结果完全无效。如果你的单倍型定义里包含位置信息（染色体、起始位、终止位），可以先写一个排序函数，按染色体号、起始位从小到大排序，再生成等位基因列表。
• 处理杂合情况：如果患者的SNP是杂合的（比如A/T），你可能需要调整匹配逻辑——比如只要单倍型的等位基因是杂合中的一个就算匹配，这时候需要修改calculate_match_score函数，比如把患者的等位基因拆分成列表，然后判断单倍型的等位基因是否在其中。
• Python 2.7兼容性：代码里用了float()确保除法是浮点型，print语句是Python2原生格式，避免了版本兼容问题。

内容的提问来源于stack exchange，提问作者Nulnulzeven