# 生物信息学视角下的IGH基因重排分析：用Python解析B细胞受体多样性 在免疫系统的精密舞台上，B细胞扮演着抗体工厂的角色，而决定其产品特异性的核心蓝图，就编码在免疫球蛋白重链（IGH）基因的重排过程中。对于从事癌症免疫治疗、自身免疫病研究或抗体药物开发的生物信息学分析人员而言，深入理解并能够从海量测序数据中解析IGH基因的重排模式，是一项极具价值且充满挑战的核心技能。这不仅仅是识别一段序列，更是解读B细胞克隆的“身份证”、追踪其演化谱系、乃至窥探疾病发生机制的关键窗口。 传统的分析方法往往依赖于商业软件或固定的分析流程，虽然便捷，但有时就像使用一个黑箱，限制了我们对数据底层逻辑的深度定制和灵活探索。而Python，凭借其强大的生态库和灵活性，为我们打开了一扇自主分析的大门。本文将从一个实践者的角度，带你用Python工具链，从原始的FASTQ或BAM文件出发，一步步实现IGH基因重排的分析、特征提取与可视化。我们将重点关注如何构建一个轻量但功能完整的分析流程，处理VDJ重组算法中的关键步骤，提取决定抗体特异性的互补决定区（CDR3）序列，并最终将抽象的序列数据转化为直观的克隆演化图谱。无论你是希望为自己的研究项目搭建定制化分析管线，还是渴望更透彻地理解免疫组库测序数据的本质，这里的思路和代码都将提供切实的参考。 ## 1. 从原始数据到VDJ注释：构建分析流水线 拿到免疫组库测序数据后，第一步往往不是急于运行某个“万能”脚本，而是理解数据的来源和格式。数据可能来自靶向IGH基因的扩增子测序，也可能是全转录组或全外显子组测序中捕获的免疫受体序列。明确这一点，关系到后续比对参考数据库的选择和质量控制策略的制定。 一个稳健的分析流程始于严格的质量控制。除了常规的测序质量、接头污染和重复序列检查，对于免疫组库数据，我们还需特别关注引物区域的完整性，因为引物区域覆盖了V基因片段的部分保守序列，其缺失或错配会严重影响后续的基因注释准确性。 > 注意：不同的实验建库方案（例如，多重PCR引物集或5‘ RACE）会直接影响数据的起始位置和可用序列长度，在设置质量控制参数时需要据此调整。 完成质控后，核心步骤是将每条测序读段比对到免疫球蛋白基因参考序列上，以确定其使用的V、D、J基因片段。虽然已有MiXCR、IgBLAST等成熟工具，但用Python实现一个简化版的理解性流程，能让我们更清晰地把握其中的算法逻辑。这里，我们可以利用`Biopython`和`pandas`来构建一个基于局部比对和得分矩阵的注释器。 首先，我们需要准备参考数据库。可以从IMGT（国际免疫遗传学信息系统）下载所有功能性的IGHV、IGHD、IGHJ基因的核苷酸序列。将其整理成FASTA格式，并存储在一个便于快速查询的数据结构中，例如字典。 ```python import pandas as pd from Bio import SeqIO from Bio import pairwise2 from Bio.SubsMat import MatrixInfo as matlist def load_imgt_database(fasta_path): """加载IMGT基因序列数据库到字典""" gene_db = {} for record in SeqIO.parse(fasta_path, "fasta"): # record.id 格式示例: IGHV1-2*01 gene_name = record.id.split('*')[0] # 去除等位基因信息，保留基因名 gene_db[gene_name] = str(record.seq).upper() return gene_db # 示例：加载V基因数据库 ighv_db = load_imgt_database("IMGT_IGHV.fasta") ``` 接下来，针对一条待注释的序列，我们需要分别与V、D、J基因库进行比对。由于D基因片段较短，且两侧会发生核苷酸的随机插入/删除（N-核苷酸添加），直接比对效果可能不佳，通常更依赖其在V-J比对后留下的中间序列特征进行推断。这里先展示V基因的比对注释： ```python def annotate_v_gene(query_seq, gene_db, matrix=matlist.blosum62, gap_open=-10, gap_extend=-0.5): """使用局部比对算法注释最可能的V基因""" best_score = -float('inf') best_gene = None best_alignment = None for gene_name, ref_seq in gene_db.items(): # 使用局部比对，寻找query_seq与ref_seq的最佳匹配区域 alignments = pairwise2.align.localds(query_seq, ref_seq, matrix, gap_open, gap_extend) if alignments: top_alignment = alignments[0] score = top_alignment[2] # 比对得分 if score > best_score: best_score = score best_gene = gene_name best_alignment = top_alignment # 可以设置一个得分阈值，低于阈值则认为未注释到可信的V基因 if best_score < 150: # 阈值需根据实际情况调整 best_gene = "unassigned" return best_gene, best_score, best_alignment # 对一条测序读段进行V基因注释示例 test_seq = "CAGGTGCAGCTGGTGCAGTCTGGAGCTGAGGTGAAGAAGCCTGGGGCCTCAGTGAAGGTCTCCTGCAAGGCTTCTGGATACACCTTCACCGACTACTATATGCACTGGGTGCGACAGGCCCCTGGACAAGGGCTTGAGTGGATGGGATGGATCAACCCTAACAGTGGTGGCACAAACTATGCACAGAAGTTTCAGGGCAGGGTCACCATGACCAGGGACACGTCCATCAGCACAGCCTACATGGAGCTGAGCAGGCTGAGATCTGACGACACGGCCGTGTATTACTGTGCGAGAGA" v_gene, v_score, v_align = annotate_v_gene(test_seq, ighv_db) print(f"注释的V基因: {v_gene}, 比对得分: {v_score}") ``` 通过类似的流程完成J基因注释后，位于V基因比对结束位置和J基因比对开始位置之间的序列，就是包含了D基因片段和N-核苷酸添加的高变CDR3区。提取这段序列是后续克隆型分析的基础。 ## 2. CDR3序列提取与克隆型聚类：定义B细胞的身份标识 互补决定区3（CDR3）是抗体可变区中多样性最高、直接参与抗原结合的关键区域。其核苷酸序列由V基因的3‘端、D基因（重链中）、J基因的5’端，以及连接处随机添加的N-核苷酸共同构成。因此，CDR3序列几乎是每个B细胞克隆独一无二的分子标签。 提取CDR3序列需要精确识别V和J基因的边界。IMGT数据库定义了保守的氨基酸锚点（如V基因末端的保守半胱氨酸C104，J基因起始的保守苯丙氨酸F118），我们可以根据比对结果，定位这些保守位点，从而框定CDR3区域。 ```python def extract_cdr3_nt(query_seq, v_alignment, j_alignment, v_gene_name): """ 根据V和J基因的比对结果，提取核苷酸水平的CDR3序列。 假设v_alignment和j_alignment是pairwise2.align.localds返回的最佳比对元组。 """ # 从V基因比对结果中，找到query序列中与IMGT位置104（保守Cys）相对应的位置。 # 这需要根据比对字符串进行解析。这里是一个简化逻辑示例： # 实际中需要更复杂的逻辑来解析比对字符串，找到特定参考位置对应的查询位置。 # 简化示例：假设我们已经通过其他方法获得了V基因比对在查询序列上的结束位置（v_end_in_query） # 和J基因比对在查询序列上的开始位置（j_start_in_query） v_end_in_query = 120 # 示例值，需实际计算 j_start_in_query = 180 # 示例值，需实际计算 if v_end_in_query < j_start_in_query: cdr3_nt = query_seq[v_end_in_query: j_start_in_query] return cdr3_nt else: return None # V和J区域重叠，无效重排 # 在实际流程中，我们需要解析比对字符串的细节，这是一个更复杂的函数示例框架 def find_junction_from_alignment(alignment, ref_seq, target_ref_pos): """ 在比对结果中，找到参考序列上特定位置（target_ref_pos）在查询序列中对应的位置。 alignment: pairwise2的比对结果元组 (seqA, seqB, score, start, end) """ aligned_seqA, aligned_seqB, score, begin, end = alignment ref_pos = 0 query_pos = 0 for a, b in zip(aligned_seqA, aligned_seqB): if a != '-': # 查询序列当前位置有碱基 query_pos_increment = 1 else: query_pos_increment = 0 if b != '-': # 参考序列当前位置有碱基 ref_pos += 1 if ref_pos == target_ref_pos: # 找到参考序列目标位置，返回此时查询序列的索引（从0开始） # 注意：这里返回的是在原始查询序列中的近似位置，需要根据begin和gap进行调整 # 此为简化逻辑，实际应用需更精确处理gap。 return begin + query_pos query_pos += query_pos_increment return None ``` 获得CDR3核苷酸序列后，我们可以将其翻译成氨基酸序列。**相同的CDR3氨基酸序列（允许轻微的测序误差）通常被视为同一个B细胞克隆**。因此，克隆型聚类就转化为对大量CDR3氨基酸序列进行相似性分组的问题。一个简单有效的方法是使用精确匹配或允许少量错配的聚类算法。 ```python from collections import defaultdict import Levenshtein # 需要安装 python-Levenshtein 包 def cluster_cdr3_aa(cdr3_aa_list, max_distance=1): """ 使用编辑距离对CDR3氨基酸序列进行聚类。 cdr3_aa_list: 列表，每个元素是一个CDR3氨基酸序列字符串。 max_distance: 归为同一克隆的最大编辑距离。 返回一个字典，键为代表性序列（簇中心），值为属于该簇的所有序列列表。 """ clusters = {} seq_to_cluster = {} # 记录每个序列属于哪个簇 for seq in cdr3_aa_list: assigned = False # 遍历现有簇的代表性序列 for rep_seq in clusters.keys(): if Levenshtein.distance(seq, rep_seq) <= max_distance: clusters[rep_seq].append(seq) seq_to_cluster[seq] = rep_seq assigned = True break # 如果未分配到任何现有簇，则创建一个新簇 if not assigned: clusters[seq] = [seq] seq_to_cluster[seq] = seq return clusters, seq_to_cluster # 示例：假设我们有一组提取的CDR3氨基酸序列 sample_cdr3s = ["CARGGNYGYDFWS", "CARGGNYGYDFWS", "CARDSSGYDFWS", "CARGGSYDFWS", "CARGGNYGYYFWS"] clusters, mapping = cluster_cdr3_aa(sample_cdr3s, max_distance=2) print(f"共形成 {len(clusters)} 个克隆簇") for rep, members in clusters.items(): print(f"代表序列: {rep}, 成员数: {len(members)}") ``` 通过克隆型聚类，我们可以将数百万条读段归结为几十到几千个克隆型，并计算每个克隆型的频率（克隆丰度），这是评估免疫组库多样性、识别优势克隆的基础。 ## 3. 多样性度量与克隆结构可视化：从数据到洞察 量化免疫组库的多样性是许多研究问题的核心。多样性并非单一概念，它至少包含以下几个维度： | 多样性维度 | 描述 | 常用指标 | | :--- | :--- | :--- | | **丰富度 (Richness)** | 克隆型的绝对数量 | 观测克隆型数 (S) | | **均匀度 (Evenness)** | 各克隆型丰度分布的均匀程度 | 香农熵 (Shannon Index)、辛普森指数 (Simpson Index)、Pielou均匀度 | | **克隆结构** | 优势克隆与稀有克隆的组成 | 克隆丰度分布曲线、Gini系数 | 我们可以用`scipy`和`numpy`轻松计算这些指标： ```python import numpy as np from scipy.stats import gini def calculate_diversity_metrics(clone_frequencies): """ clone_frequencies: 列表或数组，每个元素是一个克隆型的频率（比例或绝对数） """ freqs = np.array(clone_frequencies) total = freqs.sum() proportions = freqs / total # 1. 丰富度：克隆型数量 richness = len(freqs) # 2. 香农熵 shannon = -np.sum(proportions * np.log(proportions)) # 3. 辛普森指数 (概率论形式，值越大多样性越低) simpson = np.sum(proportions ** 2) # 4. Pielou均匀度 (J) if richness > 1: j = shannon / np.log(richness) else: j = 1.0 # 5. Gini系数 (衡量不平等性，0完全平等，1完全不平等) gini_coefficient = gini(freqs) metrics = { 'Richness': richness, 'Shannon_Index': shannon, 'Simpson_Index': simpson, 'Pielou_Evenness': j, 'Gini_Coefficient': gini_coefficient } return metrics # 示例：假设有5个克隆型，其测序读段数分别为 1000, 200, 50, 20, 5 freqs = [1000, 200, 50, 20, 5] metrics = calculate_diversity_metrics(freqs) for k, v in metrics.items(): print(f"{k}: {v:.4f}") ``` 可视化是让这些数字“说话”的关键。我们可以用`matplotlib`或`seaborn`绘制多种图形来揭示克隆结构。 **克隆丰度分布曲线（Rank-Abundance Curve）**：将克隆型按丰度从高到低排序后绘图，能直观展示优势克隆的“统治力”和长尾稀有克隆的存在。 ```python import matplotlib.pyplot as plt import seaborn as sns def plot_rank_abundance(clone_frequencies, sample_name=""): """ 绘制克隆丰度排序曲线 """ sorted_freqs = np.sort(clone_frequencies)[::-1] # 降序排列 ranks = np.arange(1, len(sorted_freqs) + 1) plt.figure(figsize=(10, 6)) plt.plot(ranks, sorted_freqs, 'o-', linewidth=2, markersize=5) plt.yscale('log') # Y轴常用对数刻度以看清稀有克隆 plt.xscale('log') # X轴有时也用对数 plt.xlabel('Clone Rank (log scale)') plt.ylabel('Clone Frequency (log scale)') plt.title(f'Rank-Abundance Curve - {sample_name}') plt.grid(True, which="both", ls="--", alpha=0.3) plt.show() # 生成模拟数据并绘图 np.random.seed(42) # 模拟一个包含少数优势克隆和大量稀有克隆的分布 sim_freqs = np.concatenate([np.array([5000, 1000, 500]), np.random.randint(1, 100, 97)]) plot_rank_abundance(sim_freqs, "Simulated Sample") ``` **克隆谱系树（Phylogenetic Tree）或网络图**：如果我们有来自同一患者不同时间点（如治疗前后）或多个部位（如肿瘤组织、外周血）的样本，可以追踪特定克隆型的动态变化。通过比较不同样本中相同CDR3序列的丰度变化，可以绘制热图或折线图来展示克隆的扩增或消退。更进一步，如果对同一克隆型的序列进行更细致的突变分析，可以构建其内部演化树，揭示体细胞超突变（SHM）的积累过程。 ```python import pandas as pd def plot_clonal_dynamics(longitudinal_data): """ longitudinal_data: DataFrame，索引为克隆型ID，列为不同时间点，值为丰度。 绘制前N个优势克隆的动态变化。 """ # 选取基线时（第一列）丰度最高的前10个克隆 top_clones = longitudinal_data.iloc[:, 0].nlargest(10).index top_data = longitudinal_data.loc[top_clones] plt.figure(figsize=(12, 8)) for clone_id in top_data.index: plt.plot(top_data.columns, top_data.loc[clone_id], 'o-', label=clone_id, linewidth=2) plt.xlabel('Time Point') plt.ylabel('Clone Frequency (Reads per Million)') plt.title('Longitudinal Dynamics of Top 10 Clones') plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left') plt.grid(True, alpha=0.3) plt.tight_layout() plt.show() # 示例数据框构建 time_points = ['Pre', 'Post_1m', 'Post_3m', 'Post_6m'] clone_ids = [f'Clone_{i}' for i in range(100)] # 模拟数据：每个克隆在不同时间点有一个丰度值 np.random.seed(123) sim_dynamics = pd.DataFrame( np.random.randn(100, 4).cumsum(axis=1) + 50, # 模拟有趋势的随机变化 index=clone_ids, columns=time_points ) sim_dynamics = sim_dynamics.abs() # 取绝对值确保非负 plot_clonal_dynamics(sim_dynamics) ``` ## 4. 高级分析：体细胞超突变与克隆演化 在B细胞发育的生发中心阶段，活化的B细胞会对其IGHV基因引入高频的点突变，这一过程称为体细胞超突变（SHM），是抗体亲和力成熟的核心机制。分析SHM的模式和程度，对于理解免疫应答的质量、区分记忆B细胞亚群以及研究某些淋巴瘤的发病机制至关重要。 分析SHM，首先需要将测序读段与推断出的胚系V基因序列进行精确比对，识别出所有突变位点。突变通常集中在互补决定区（CDR），尤其是CDR1和CDR2，但框架区（FR）也会发生。 ```python def calculate_shm_rate(query_seq, germline_v_seq): """ 计算给定序列相对于其胚系V基因序列的突变率。 假设query_seq和germline_v_seq已经过比对（长度相同，可能包含gap）。 返回突变总数和突变率（突变数/可比对长度）。 """ if len(query_seq) != len(germline_v_seq): raise ValueError("Sequences must be aligned and of equal length.") mutations = 0 comparable_positions = 0 for q, g in zip(query_seq, germline_v_seq): if g != '-': # 胚系序列在该位置有碱基，为可比对位置 comparable_positions += 1 if q != g and q != '-': # 查询序列碱基不同且非缺失，计为突变 mutations += 1 if comparable_positions > 0: mutation_rate = mutations / comparable_positions else: mutation_rate = 0.0 return mutations, mutation_rate, comparable_positions # 示例：假设我们已经有了比对后的序列（包含gap） aligned_query = "CAGGTGCAGCTGGTGCAGTCTGGAGCTGAGGTGAAGAAGCCTGGGGCCTCAGTGAAGGTCTCCTGCAAGGCTTCTGGATACACCTTCACCGACTACTATATGCACTGGGTGCGACAGGCCCCTGGACAAGGGCTTGAGTGGATGGGATGGATCAACCCTAACAGTGGTGGCACAAACTATGCACAGAAGTTTCAGGGCAGGGTCACCATGACCAGGGACACGTCCATCAGCACAGCCTACATGGAGCTGAGCAGGCTGAGATCTGACGACACGGCCGTGTATTACTGTGCGAGAGA" aligned_germline = "CAGGTGCAGCTGGTGCAGTCTGGGGCTGAGGTGAAGAAGCCTGGGGCCTCAGTGAAGGTCTCCTGCAAGGCTTCTGGATACACCTTCACCAGCTACTATATGCACTGGGTGCGACAGGCCCCTGGACAAGGGCTTGAGTGGATGGGATGGATCAACCCTAACAGTGGTGGCACAAACTATGCACAGAAGTTTCAGGGCAGGGTCACCATGACCAGGGACACGTCCATCAGCACAGCCTACATGGAGCTGAGCAGGCTGAGATCTGACGACACGGCCGTGTATTACTGTGCGAGAGA" # 注意：这里为了示例，人为修改了germline的一处（GGG -> GAG），实际应从数据库获取精确胚系序列。 mutations, rate, length = calculate_shm_rate(aligned_query, aligned_germline) print(f"突变数: {mutations}, 可比对长度: {length}, 突变率: {rate:.4f}") ``` 更进一步，我们可以分析突变的频谱，即不同碱基转换（如A>G， C>T）和颠换的比例。这有助于判断突变是否由特定的脱氨酶（如AID）活性驱动，其典型的突变模式是偏向于在WRC（W=A/T, R=A/G）热点基序上发生C>T或G>A突变。 ```python from collections import Counter def analyze_mutation_spectrum(query_seq, germline_seq): """ 分析突变类型频谱。 返回一个计数器，记录各种突变类型（如 'A>G', 'C>T'）的数量。 """ spectrum = Counter() for q, g in zip(query_seq, germline_seq): if g != '-' and q != '-' and q != g: mutation = f"{g}>{q}" spectrum[mutation] += 1 return spectrum spectrum = analyze_mutation_spectrum(aligned_query, aligned_germline) print("突变频谱:") for mut, count in spectrum.most_common(): print(f" {mut}: {count}") ``` 将SHM分析与克隆型聚类结合，我们可以构建**克隆内演化树**。同一个克隆（共享相同CDR3）的不同序列，可能因为积累了不同的SHM而分化。我们可以将这些序列的多态性位点提取出来，利用邻接法（Neighbor-Joining）或最大简约法构建系统发育树，可视化克隆内部的微演化路径。 ```python from Bio.Phylo.TreeConstruction import DistanceCalculator, DistanceTreeConstructor from Bio.Phylo import draw from Bio import AlignIO import io def build_intraclonal_tree(sequence_list, clone_id): """ 为同一克隆内的多条序列构建系统发育树。 sequence_list: 列表，每个元素是（已比对到同一胚系V基因的）序列字符串。 """ # 将序列列表转换为多序列比对的格式（例如，写入临时字符串） aligned_seqs = [] for i, seq in enumerate(sequence_list): # 这里假设所有序列已经与同一个胚系序列对齐，长度一致 aligned_seqs.append(f">Seq_{i}\n{seq}") fasta_data = "\n".join(aligned_seqs) # 使用Bio.AlignIO读取 alignment = AlignIO.read(io.StringIO(fasta_data), "fasta") # 计算距离矩阵（使用简单的核酸身份差异） calculator = DistanceCalculator('identity') dm = calculator.get_distance(alignment) # 使用邻接法构建树 constructor = DistanceTreeConstructor() tree = constructor.nj(dm) # 绘制树（简单文本显示或图形显示） print(f"Intra-clonal tree for {clone_id}:") # draw.ascii(tree) # 在控制台打印ASCII树 # 或者使用matplotlib绘制 # fig = plt.figure(figsize=(10, 8)) # axes = fig.add_subplot(1, 1, 1) # draw(tree, axes=axes) # plt.title(f'Intra-clonal Phylogeny: {clone_id}') # plt.show() return tree # 示例：假设我们有一个克隆内的3条略有差异的序列 clonal_variants = [ aligned_query, # 原始序列 aligned_query.replace('A', 'G', 1), # 模拟一个A>G突变 aligned_query.replace('C', 'T', 1), # 模拟一个C>T突变 ] tree = build_intraclonal_tree(clonal_variants, "Clone_X") ``` 通过这些高级分析，我们不仅能回答“有哪些克隆”，更能深入探究“这些克隆是如何进化而来的”，从而在癌症免疫治疗中，帮助识别那些经过亲和力成熟、具有潜在抗肿瘤活性的T细胞克隆，或在淋巴瘤研究中，揭示肿瘤克隆的演化历史和驱动突变。