肿瘤驱动基因鉴定¶

# === 使用 MutSigCV 鉴定驱动基因 ===

# 1. 准备输入文件（MAF 格式）
# MAF 文件包含所有样本的突变信息
# 通常从 Mutect2 → VEP → maf_tools 转换

# 2. 运行 MutSigCV（MATLAB 版本）
# MutSigCV 需要三个输入：突变文件、覆盖度文件、协变量文件
matlab -nodisplay -r \
  "MutSigCV('mutations.maf', 'coverage.txt', \
   'covariates.txt', 'mutsig_output')"

# === 使用 dNdScv（R包，更现代的方法） ===

# === dNdScv：基于 dN/dS 的驱动基因鉴定 ===
library(dndscv)                        # 加载dNdScv包

# 1. 准备输入数据（最少4列）
# sampleID, chr, pos, ref, alt
mutations <- read.table(               # 读取突变文件
  "somatic_mutations.txt",             # 突变文件
  header = TRUE                        # 有表头
)

# 2. 运行 dNdScv
result <- dndscv(                      # 运行分析
  mutations,                           # 突变数据
  refdb = "hg38",                      # 参考基因组版本
  max_muts_per_gene_per_sample = 3,    # 每样本每基因最多3个突变
  max_coding_muts_per_sample = 5000    # 每样本最多编码突变数
)

# 3. 查看显著驱动基因（q值<0.1）
sig_genes <- result$sel_cv             # 基因选择结果
sig_drivers <- sig_genes[
  sig_genes$qallsubs_cv < 0.1,        # 全局q值<0.1
]
print(sig_drivers[, c("gene_name",     # 打印基因名
                       "n_syn",         # 同义突变数
                       "n_mis",         # 错义突变数
                       "n_non",         # 无义突变数
                       "wmis_cv",       # 错义dN/dS
                       "wnon_cv",       # 无义dN/dS
                       "qallsubs_cv")]) # q值

# 4. 判断驱动类型（癌基因 vs 抑癌基因）
for (i in 1:nrow(sig_drivers)) {       # 遍历每个驱动基因
  gene <- sig_drivers$gene_name[i]     # 基因名
  wmis <- sig_drivers$wmis_cv[i]       # 错义dN/dS
  wnon <- sig_drivers$wnon_cv[i]      # 无义dN/dS

  if (wmis > 1 & wnon <= 1) {         # 错义高、无义不高
    cat(gene, "→ 可能是癌基因\n")      # 功能获得型
  } else if (wnon > 1) {              # 无义突变显著
    cat(gene, "→ 可能是抑癌基因\n")    # 功能丧失型
  }
}

# 使用 OncoKB API 注释驱动突变的临床意义
import requests                        # HTTP请求库
import pandas as pd                    # 数据处理

def annotate_oncokb(mutations_df, token):
    """用 OncoKB 注释突变的临床可操作性"""
    url = "https://www.oncokb.org/api/v1/annotate/mutations/byProteinChange"
    headers = {"Authorization": f"Bearer {token}"}  # API密钥

    results = []                       # 存储注释结果
    for _, row in mutations_df.iterrows():  # 遍历突变
        params = {
            "hugoSymbol": row["gene"],     # 基因名
            "alteration": row["protein_change"],  # 蛋白变化
            "tumorType": row["cancer_type"],      # 癌种
        }
        resp = requests.get(url, headers=headers, params=params)
        if resp.status_code == 200:    # 请求成功
            data = resp.json()
            results.append({
                "gene": row["gene"],
                "mutation": row["protein_change"],
                "oncogenic": data.get("oncogenic", ""),      # 是否致癌
                "level": data.get("highestSensitiveLevel", ""),  # 最高证据等级
                "drugs": data.get("treatments", [])           # 对应药物
            })
    return pd.DataFrame(results)       # 返回注释表

# 使用示例
# result = annotate_oncokb(mutations, "YOUR_ONCOKB_TOKEN")