拿到一份肿瘤样本的MAF文件很多人的第一反应是看VAF列——0.2、0.3、0.45……然后习惯性地给低频突变贴上“亚克隆”标签把接近0.5的位点当成交互验证后的“克隆性驱动事件”。这个判断方向没错但如果样本纯度只有40%一个VAF0.2的位点很可能其实是所有肿瘤细胞都携带的克隆性突变反过来纯度90%的样本里一个VAF0.35的位点也可能只是某条亚克隆分支上的事件。VAF、MAF、CCF这三个词在肿瘤基因组分析里几乎天天见真正能把它们算清楚、讲明白的人却没那么多。这篇文章就从我的实际分析经验出发把这三个指标的计算逻辑、Python实现和应用场景一次说透。1. 三个指标到底在描述什么先弄清楚VAF、MAF、CCF的角色1.1 从一条测序read开始理解VAFVAF的全称是Variant Allele Frequency变异等位基因频率。它描述的是在某个基因组位点上覆盖到该位点的所有测序read里有多少比例支持变异等位基因。举个例子某一位点一共被测序read覆盖了100次其中35条read携带了变异碱基65条read是参考碱基那这个位点的VAF就是35/1000.35。这个数字本身是一个非常底层的观测值它受到三方面因素的共同影响样本中肿瘤细胞的比例纯度purity、该位点在肿瘤细胞中的拷贝数状态CN、以及真正携带该突变的肿瘤细胞比例CCF。正因为VAF是这三者的混合体所以单纯看VAF高低做克隆性判断在低纯度样本或拷贝数扩增区域会得出很离谱的结论。1.2 MAF一表揽尽所有变异信息的格式MAF在这里不是Minor Allele Frequency而是Mutation Annotation Format是一种表格化的变异注释文件格式。TCGA当年为了让各个中心的变异检出结果可以横向比较把VCF中杂乱的INFO字段、注释字段抽出来整理成一行一个变异位点的MAF文件。MAF文件的意义在于它把“这个样本有哪些突变、每个突变有多少read支持、注释到了什么基因、产生了什么氨基酸改变”这些信息全部标准化了。只要会读表格就能对一份肿瘤样本的突变全貌有一个直观认识。后面我们会专门花一章讲MAF的字段和解析方法。1.3 CCF回答“肿瘤异质性”的关键指标CCF全称Cancer Cell Fraction癌细胞分数。它回答的问题是在这份肿瘤样本的癌细胞群体中有多大比例的癌细胞携带某个特定突变。如果某个突变出现在肿瘤发生早期存在于所有癌细胞里那它的CCF理论上就是1我们称之为克隆性突变。如果某个突变只存在于一部分癌细胞中比如某个亚克隆分支上那它的CCF就小于1我们称之为亚克隆突变。CCF的计算需要用到VAF、肿瘤纯度、局部拷贝数、突变等位基因拷贝数四个参数。这就是为什么不能只靠VAF判断克隆性VAF是表面观测值CCF是校正了纯度和拷贝数之后的“细胞比例”。2. VAF计算从BAM到数字哪些细节决定了你的结果是否可靠2.1 基础公式与最小可用代码VAF的基础公式非常简单VAF alt_count / (ref_count alt_count)但到了真实项目里你会发现“ref_count alt_count”和“总深度”经常对不上。有些VCF的FORMAT里DP字段包含的是所有比对到该位点的read数而AD字段只包含allele-specific的支持数。中间差出来的部分可能是其他非ref/alt的碱基、比对质量不过滤的read、或indel附近的复杂比对。我的建议是既然要算VAF就用ref和alt的计数来算不要直接用DP做分母。下面是最小可用代码用Python直接解析VCF并计算每个突变位点的VAFimport gzip def parse_vcf_and_calc_vaf(vcf_path): results [] opener gzip.open if vcf_path.endswith(.gz) else open with opener(vcf_path, rt) as f: for line in f: if line.startswith(#): continue fields line.strip().split(\t) chrom, pos, ref, alt fields[0], int(fields[1]), fields[3], fields[4] # 只取FORMAT为AD:DP的情况或直接取AD字段 format_keys fields[8].split(:) sample_values fields[9].split(:) sample_dict dict(zip(format_keys, sample_values)) if AD not in sample_dict: continue # AD格式ref_count,alt_count ad_values sample_dict[AD].split(,) ref_count int(ad_values[0]) alt_count int(ad_values[1]) if len(ad_values) 1 else 0 total_count ref_count alt_count vaf alt_count / total_count if total_count 0 else float(nan) results.append({ chrom: chrom, pos: pos, ref: ref, alt: alt, ref_count: ref_count, alt_count: alt_count, vaf: round(vaf, 4) }) return results这段代码只处理单个样本多样本VCF要加一层循环。2.2 计算前的过滤没有质控的VAF就是噪声直接拿VCF的AD字段算VAF你会得到很多假信号。我自己的标准流程在计算VAF之前会做一轮过滤总深度不足的位点直接丢弃。点突变一般要求refalt≥20深度30x以下的WES数据里低于这个阈值的位点VAF波动非常大。变异支持读数太少。即使总深度达到100只有2条read支持变异这个位点也基本不可信。通常要求alt_count≥3或以上严格一点要求≥5。链偏向性过滤。如果绝大多数alt read都来自同一条链大概率是比对假象。可以用Fisher精确检验或者简单的链偏向比例SOR来判断。同义突变和已知胚系位点如果我们要做的是体细胞突变分析通常已经用配对正常样本过滤过胚系突变如果没有配对样本需要谨慎处理。这些过滤条件不是为了让你把数据“洗得更干净”而是避免后续CCF计算时把一个噪声点当成亚克隆事件去分析。2.3 深度和等位基因计数的定义差异最容易埋雷的地方不同caller输出的等位基因计数含义并不一致这是实际计算VAF时最大的坑。CallerAD字段含义注意事项Mutect2AD[0]为refAD[1]为alt可能有AD[2]以上总深度DP可能大于AD之和Strelka2没有AD用AU/CU/GU/TU四个字段需要取变异碱基对应的支持数VarScan2AD和DP都有注意某个等位基因计数为0的情况GATK UnifiedGenotyperADref,alt另有DP老工具但数据集中仍常见我遇到过最坑的情况是某个VCF文件用的是Mutect2流程但样本的FORMAT里除了AD还有一个叫DP的字段而DP的值和AD两个数加起来差了十几条read。一查文档才发现这些read被判定为“其他等位基因”或低质量read被AD排除但没有从DP中排除。用DP做分母算出来的VAF和用AD总和做分母算出来的VAF能差出5个百分点。我个人的习惯是解析VCF之后先用AD字段算出VAF再和VCF自带的VAF列如果有的话做交叉校验。一旦发现差异超过0.05就回去检查这个caller对深度的定义。3. MAF不是指标是标准化的“中间语言”解析与实战3.1 从VCF到MAF为什么需要一个标准化的中间层VCF格式本身信息量很大但缺点是每个变异检测工具输出的VCF在注释字段上风格差异巨大。有的在INFO里塞了一大堆自定义flag有的把群体频率、蛋白功能预测散在不同标签里。这时候拿MAF做统一承接就非常省事。MAF的每一行代表一个样本中的一个变异位点关键字段固定后续无论是统计VAF分布、画oncoplot还是做CCF计算都可以直接读入pandas操作不需要再跟VCF的FORMAT较劲。3.2 核心字段逐项说明MAF文件通常有很多列但真正在VAF/CCF分析里需要关心的核心字段就十几个。字段含义注意事项Chromosome染色体有的文件带chr前缀有的不带Start_Position起始坐标注意是1-based坐标End_Position终止坐标用于indel时需和Start一起看Reference_Allele参考等位基因Tumor_Seq_Allele2肿瘤样本中的变异等位基因t_ref_count肿瘤样本中参考等位基因read数VAF计算的分子分母来源t_alt_count肿瘤样本中变异等位基因read数n_ref_count配对正常样本中参考等位基因read数用于判断位点在正常样本中是否存在n_alt_count配对正常样本中变异等位基因read数Variant_Classification突变类型Missense_Mutation、Nonsense_Mutation等Variant_TypeSNP/INS/DEL等indel做VAF/CCF时要特别小心Tumor_Sample_Barcode样本识别号多样本MAF中不要串了HGVSp_Short蛋白变更注释如p.E746_A750delTCGA的MAF规范里还有几十个其他列包括dbSNP ID、COSMIC ID、ExAC_AF等但核心分析用上面这些字段就够了。3.3 用Python批量读取MAF并计算VAF读取MAF文件用pandas即可但要注意读的时候把t_ref_count、t_alt_count这些列转成数值类型否则后面做除法会报错或得到object类型。import pandas as pd def load_maf_and_calc_vaf(maf_path, sample_colTumor_Sample_Barcode): maf pd.read_csv(maf_path, sep\t, comment#, low_memoryFalse) # 只保留分析必需字段 required_cols [ Chromosome, Start_Position, End_Position, Reference_Allele, Tumor_Seq_Allele2, t_ref_count, t_alt_count, n_ref_count, n_alt_count, Variant_Classification, Variant_Type, HGVSp_Short ] available_cols [c for c in required_cols if c in maf.columns] maf maf[available_cols [sample_col] if sample_col in maf.columns else available_cols] # 计算tumor VAF maf[t_depth] maf[t_ref_count] maf[t_alt_count] maf[VAF] maf[t_alt_count] / maf[t_depth] # 计算normal VAF用于检查是否可能是胚系突变 maf[n_depth] maf[n_ref_count] maf[n_alt_count] maf[normal_VAF] maf[n_alt_count] / maf[n_depth] return maf # 使用示例 maf_df load_maf_and_calc_vaf(sample.maf) print(maf_df[[Chromosome, Start_Position, HGVSp_Short, VAF, normal_VAF]].head())这段代码会把MAF里的t_ref_count和t_alt_count加起来作为深度。需要注意如果MAF文件是某个工具自产自销的可能用的不是t_ref_count/t_alt_count这种列名而是Tumor_Ref_Count/Tumor_Alt_Count之类的变体。读文件之前最好先打印一下列名确认。4. CCF推导到Python实现把VAF真正换算成癌细胞分数4.1 从期望VAF反推CCF公式推导全程CCF的核心思想是反推。我们先回答一个问题如果某个突变存在于所有癌细胞中克隆性突变在已知纯度purity、局部拷贝数CN_tumor、突变等位基因拷贝数m的情况下理论上应该观测到多大的VAF正常细胞部分贡献的等位基因来自正常二倍体假设CN_normal2那么总的等位基因数归一化为总等位基因数 purity × CN_tumor (1 - purity) × CN_normal突变等位基因的“总量”为突变等位基因数 purity × m所以如果所有肿瘤细胞都携带该突变期望的VAF为expected_VAF (purity × m) / (purity × CN_tumor (1 - purity) × CN_normal)而我们实际观测到的VAF可能比这个小说明只有部分肿瘤细胞携带或比这个大说明存在拷贝数变异或纯度估计偏差。CCF的定义就是实际VAF占期望VAF的比例CCF observed_VAF / expected_VAF observed_VAF × (purity × CN_tumor (1 - purity) × CN_normal) / (purity × m)这个公式中observed_VAF是直接由测序深度算出来的purity、CN_tumor、m需要从其他渠道获取。4.2 纯度、拷贝数、m值从哪里来这一步是整个分析里最容易被忽略、也最容易出问题的地方。纯度purity表示样本中肿瘤细胞所占比例。临床病理切片可以估计但是比较粗糙更常用的方式是用ABSOLUTE、ASCAT、sequenza等工具基于WES/WGS的拷贝数信号估算纯度和倍性。如果项目里用的是公共数据例如TCGA有些数据门户会直接给出纯度估计值。局部拷贝数CN_tumor指突变位点所在基因组区段在肿瘤细胞中的总拷贝数。注意是包含参考等位基因和变异等位基因的总拷贝数。这个数可以从ASCAT、sequenza、FACETS的输出中获得通常是分段水平的拷贝数不是每个位点单独估计。突变等位基因拷贝数m指突变重复的次数。如果一个位点是杂合突变拷贝数是2且没有发生LOHm1如果拷贝数是3其中2条拷贝带突变m2。这个参数通常通过“依据VAF估计m取值”的方式确定比如比较观测VAF和不同m假设下的期望VAF选择最接近的那一个。正常拷贝数CN_normal通常为2但如果突变在X或Y染色体上且样本来自男性这个值会变为1。很多人会忽略这一步导致性染色体位点的CCF全部偏高。4.3 完整Python代码从MAF到CCF的一站式实现下面是一个功能完整的CCF计算函数输入MAF文件和purity输出每个位点的CCF。拷贝数信息需要单独传入我这里用一个简单的示例说明import pandas as pd import numpy as np def calc_ccf_from_vaf(vaf, purity, cn_tumor, m1, cn_normal2): 计算单个位点的CCF。 参数 vaf: 观测到的VAF0~1之间 purity: 肿瘤纯度0~1之间 cn_tumor: 该位点在肿瘤细胞中的总拷贝数 m: 突变等位基因拷贝数默认1 cn_normal: 正常细胞拷贝数常染色体为2 返回 ccf: 癌细胞分数 denominator purity * m if denominator 0: return float(nan) numerator vaf * (purity * cn_tumor (1 - purity) * cn_normal) ccf numerator / denominator return ccf def add_ccf_to_maf(maf_df, purity, cn_tableNone, cn_normal2): 给MAF加上CCF列。 cn_table: DataFrame包含chrom、start、end、cn_tumor列可为空。 如果cn_table为空假设中性拷贝数cn_tumor cn_normal。 df maf_df.copy() # 默认中性拷贝数 df[cn_tumor] cn_normal df[m] 1 if cn_table is not None: for idx, row in df.iterrows(): chrom str(row[Chromosome]) pos row[Start_Position] # 找到该位点所属的拷贝数区段 matched cn_table[ (cn_table[chrom].astype(str) chrom) (cn_table[start] pos) (cn_table[end] pos) ] if len(matched) 0: df.at[idx, cn_tumor] matched.iloc[0][cn_tumor] # 简单估计m如果cn_tumor大于2则可能有多条拷贝带突变 # 更严谨的方式是比较不同m取值下的CCF取0~1范围内的最大m max_m int(row[cn_tumor]) best_m 1 for m_try in range(1, max_m 1): ccf_try calc_ccf_from_vaf( row[VAF], purity, row[cn_tumor], m_try, cn_normal ) if 0 ccf_try 1.2: best_m m_try break df.at[idx, m] best_m df[CCF] df.apply( lambda r: calc_ccf_from_vaf( r[VAF], purity, r[cn_tumor], r[m], cn_normal ), axis1 ) return df # 使用示例 purity 0.8 maf_with_ccf add_ccf_to_maf(maf_df, purity) print(maf_with_ccf[[HGVSp_Short, VAF, CCF]].head())这段代码为了演示做了很多简化特别是对m的估计逻辑。真实项目中不建议用这种“找到第一个落在0~1.2范围内的m”的方式而应该用概率模型或者启发式规则。但在快速探索分析中这个函数已经足够给你一个位点之间相对可比的CCF值。4.4 CCF1和0的处理策略CCF计算出来大于1并不是代码写错了而是意味着“在给定的纯度和拷贝数假设下观测到的VAF超过了所有肿瘤细胞都携带该突变时的理论期望值”。这种情况在真实数据中非常常见尤其集中在拷贝数发生杂合性缺失LOH的区域或者纯度估计偏低时。我处理CCF1的方式是不直接截断为1而是先检查一下这个位点的拷贝数、m假设是否合理。如果m被低估了比如实际上是m2但我们按m1算CCF就会虚高。把m修正之后CCF通常会回落到1附近。如果所有假设都合理但CCF还是大于1.5甚至2那就要警惕样本污染、纯度估计是否离谱、或者该位点处于一个未被识别的拷贝数扩增区。CCF0的情况基本只会在允许负值时出现实际计算时用max(0, ccf)处理即可。5. 实战复盘从TCGA公开数据到自测样本我踩过的坑和经验5.1 三个文件的对齐问题样本ID、染色体命名、坐标体系在你把MAF、拷贝数结果、纯度信息拼在一起之前请先用下面的顺序做一次全面体检否则后面算出来的CCF很可能免费送你一堆低质量亚克隆突变样本IDMAF里的Tumor_Sample_Barcode和拷贝数文件里的样本ID、纯度文件里的样本ID三者必须完全一致。很多公开数据集的样本ID带有重复后缀比如“-01”和“-01A”直接join的时候会错失所有匹配。染色体命名chromosome有没有chr前缀必须统一。拷贝数文件来自ASCAT时通常不带chr而某些MAF文件里带chr。我习惯把所有染色体名统一成不带chr的格式再做匹配。坐标体系MAF的Start_Position是1-based坐标而ASCAT、sequenza输出的区段坐标有的是0-based。一个不太起眼的off-by-one错误会让边界位点匹配到错误的拷贝数区段CCF计算自然就错了。5.2 单个位点的CCF只是起点亚克隆结构的聚类思路很多人拿到CCF之后就直接用阈值比如CCF0.8算亚克隆给每个位点定性。这样做虽然省事但忽略了肿瘤亚克隆结构的一个重要特点同一个亚克隆分支上的突变CCF会落在一个相似的区间内而不是恰好等于某个精确值。所以更严谨的做法是做CCF聚类。方法上可以用PyClone、SciClone这类专门工具也可以用简单的Gaussian Mixture Model对突变位点的CCF分布做聚类。如果你不想引入太多复杂依赖直接在Python里用sklearn的GaussianMixture做一个快速聚类也是可行的from sklearn.mixture import GaussianMixture import numpy as np ccf_values maf_with_ccf[CCF].dropna().values.reshape(-1, 1) # 这里组分数量可以根据需要调整一般有2-4个亚克隆群体 gmm GaussianMixture(n_components3, random_state42) labels gmm.fit_predict(ccf_values) # 查看每个聚类的CCF均值 for i in range(gmm.n_components): cluster_ccf ccf_values[labels i] print(fCluster {i}: mean CCF {cluster_ccf.mean():.3f}, size {len(cluster_ccf)})这个思路特别适合在正式分析之前快速了解样本的亚克隆结构。需要注意的是CCF聚类时最好只使用点突变SNV因为indel在拷贝数状态估计和突变拷贝数估计上的不确定性更大聚类结果容易带偏。5.3 实操心得做这步分析时的几个小技巧把这几条经验放在最后都是我在不同项目里真金白银换来的第一纯度不一致时不要用固定VAF阈值过滤位点。同一个克隆性突变在纯度70%的样本里VAF可能在0.35左右在纯度30%的样本里就只有0.15。如果你在第一步就用VAF0.1过滤掉所有突变低纯度样本里真正的克隆性驱动事件会被直接扔掉。低频位点过滤应该放在CCF计算之后用CCF阈值过滤而不是VAF阈值。第二做CCF计算时尽量优先使用配对正常样本已验证的体细胞突变。没有正常样本的情况下很多低频突变其实是残留的胚系杂合位点。这类位点在CCF计算时往往集中在某个固定VAF附近聚类时会形成一个人造的“假亚克隆群体”。第三当m值不确定时尽量查看突变所在区域的等位基因特异性拷贝数。比如SEG文件里如果Allele-Specific Copy Number显示某区域有LOH那么突变等位基因拷贝数m很可能等于总拷贝数CN_tumor。只是套一个m1的固定值很多LOH区域的突变CCF会系统性偏高。第四纯度最好从多个来源交叉验证。如果ASCAT、sequenza和病理估计三者的纯度差别超过15%建议检查样本是不是混合了太多正常细胞或者算法选择的倍性不对。纯度值是整个CCF公式里最敏感的参数它在分子里和分母里都出现对结果的扰动非常直接。第五公开数据集的MAF并不都是“干净”的。有的数据集把cosmic、dbNSFP、gnomAD等几十个注释列全部塞进去稍不留神就会把注释相关的字符串列当成数值列处理。处理之前把列名全部print出来检查一遍永远比想当然地按列名索引更稳妥。在真实项目中VAF到CCF只是肿瘤异质性分析的第一步。后面的亚克隆结构重建、驱动事件时序推断都是建立在这几个指标算得准、算得合理的前提之上。希望这篇实战笔记能帮你少踩几个坑。