Pangolin剪接位点有害性预测
#software
Pangolin是一个剪接位点有害性预测工具。在发表文章里，Pangolin比SpliceAI更强。
Broad研究所的SpliceAI在线工具，已经同步接入了Pangolin得分。
SpliceAI
另外，我建议大家都转到Pangolin来。因为，SpliceAI无论是源代码、抑或是Illumina提供的预计算的打分结果，都是禁止商业使用的。这意味着，就算你不使用Illumina的预计算结果，每次结果都使用SpliceAI程序来重新计算得分，在商业用途上，也是违反使用协议的。
Pangolin
但Pangolin不同，他的代码使用了GPLv3进行开源。虽然着意味着，如果你在流程里使用的是实时计算的模式，那么你的流程代码也需要被GPLv3传染而开源。但是，你完全可以使用预计算的方式计算出得分，然后作为数据库进行注释。这个数据库，是分离的，不需要开源。
我在Zenodo上找到一套基于hg38的SNV的Pangolin预计算打分，这套数据是CC 4.0协议，商用友好。
由于Pangolin，只能去计算SNV以及简单的InDel（单独的片段插入和缺失，即不能计算CGACGTACG -> CATG 这种复杂类型），同时在预计算时，很难去找到合适的InDel进行推算，因此这套数据只有SNV。
CoLab
我是打算直接使用上面Zenodo中找到的预计算数据（需要进一步处理），但是，我尝试使用Google大善人的CoLab，进行预计算。
我的Google Drive空间是2TB，足够放下大量的碎片数据，当前的流程是这样的。
- 在Google Drive中新建下面的结构
- Pangolin_Project
- inputs
- outputs
- reference
- 在CoLab里，首先运行下面代码安装Pangolin
# 1. 获取 Pangolin 的最新源码
!git clone https://github.com/tkzeng/Pangolin.git
# 2. 安装 Pangolin 运行所需的依赖库
# (注意：原版 PyVCF 在新版 Python 中容易报错，这里我们使用修复版的 PyVCF3)
!pip install PyVCF3 gffutils biopython pandas pyfastx -q
# 3. 安装 Pangolin 到系统环境变量中
!cd Pangolin && pip install . -q
# 4. 验证是否安装成功
!pangolin -h
- 在CoLab里，运行下面的代码，下载数据库和测试
import os
import shutil
import subprocess
from google.colab import drive
# ==========================================
# 第一步：挂载 Google Drive 与 初始化目录
# ==========================================
print("正在挂载 Google Drive...")
drive.mount('/content/drive')
# Drive 上的“仓库”路径
DRIVE_BASE = '/content/drive/MyDrive/Pangolin_Project'
INPUT_DIR = os.path.join(DRIVE_BASE, 'inputs')
OUTPUT_DIR = os.path.join(DRIVE_BASE, 'outputs')
REF_DIR = os.path.join(DRIVE_BASE, 'reference')
# 自动在云盘中创建所需文件夹
for d in [INPUT_DIR, OUTPUT_DIR, REF_DIR]:
os.makedirs(d, exist_ok=True)
# Colab 上的“高速车间”路径
LOCAL_WORKSPACE = '/content/local_workspace'
os.makedirs(LOCAL_WORKSPACE, exist_ok=True)
# ==========================================
# 第二步：基座文件自动检测与下载
# ==========================================
DB_NAME = 'gencode.v38.annotation.db'
FASTA_NAME = 'Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz'
DRIVE_DB = os.path.join(REF_DIR, DB_NAME)
DRIVE_FASTA = os.path.join(REF_DIR, FASTA_NAME)
# 强制下载的直链（Dropbox链接已添加 ?dl=1）
DB_URL = "https://www.dropbox.com/sh/6zo0aegoalvgd9f/AADOhGYJo8tbUhpscp3wSFj6a/gencode.v38.annotation.db?dl=1"
FASTA_URL = "https://ftp.ensembl.org/pub/release-115/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz"
print("\n--- 正在检查云盘基座文件 ---")
if not os.path.exists(DRIVE_DB):
print(f"☁️ 未检测到 {DB_NAME}，正在从 Dropbox 自动下载...")
# 使用 wget 下载并显示进度条
subprocess.run(["wget", "-q", "--show-progress", "-O", DRIVE_DB, DB_URL], check=True)
print("✅ DB 注释库下载完成！")
else:
print(f"✅ 云盘中已存在 {DB_NAME}")
if not os.path.exists(DRIVE_FASTA):
print(f"☁️ 未检测到 {FASTA_NAME}，正在从 Ensembl 自动下载...")
subprocess.run(["wget", "-q", "--show-progress", "-O", DRIVE_FASTA, FASTA_URL], check=True)
print("✅ 参考基因组下载完成！")
else:
print(f"✅ 云盘中已存在 {FASTA_NAME}")
# ==========================================
# 第三步：将基座文件缓存到 Colab 本地 (提升 I/O 速度)
# ==========================================
print("\n--- 正在将基座文件缓存至本地高速车间 ---")
LOCAL_DB = os.path.join(LOCAL_WORKSPACE, DB_NAME)
LOCAL_FASTA = os.path.join(LOCAL_WORKSPACE, FASTA_NAME)
if not os.path.exists(LOCAL_DB):
print("正在复制 DB...")
shutil.copy2(DRIVE_DB, LOCAL_DB)
if not os.path.exists(LOCAL_FASTA):
print("正在复制 FASTA...")
shutil.copy2(DRIVE_FASTA, LOCAL_FASTA)
print("✅ 本地缓存完毕！")
# ==========================================
# 第四步：核心循环 - 蚂蚁搬家与断点续传
# ==========================================
print("\n--- 开始执行 Pangolin 调度任务 ---")
vcf_files = [f for f in os.listdir(INPUT_DIR) if f.endswith('.csv')]
vcf_files.sort()
if len(vcf_files) == 0:
print("⚠️ 任务队列为空！请在云盘的 Pangolin_Project/inputs/ 目录下放入切片好的 VCF 文件。")
else:
print(f"共发现 {len(vcf_files)} 个任务切片。")
for vcf_name in vcf_files:
base_name = vcf_name.replace('.vcf', '').replace('.gz', '')
expected_output_marker = os.path.join(OUTPUT_DIR, f"{base_name}.done")
# 【断点续传逻辑】
if os.path.exists(expected_output_marker):
print(f"⏭️ [跳过] {vcf_name} 已计算完毕。")
continue
print(f"▶️ [开始处理] {vcf_name} ...")
# 从 Drive 拉取小 VCF 到本地
drive_vcf_path = os.path.join(INPUT_DIR, vcf_name)
local_vcf_path = os.path.join(LOCAL_WORKSPACE, vcf_name)
shutil.copy2(drive_vcf_path, local_vcf_path)
local_out_prefix = os.path.join(LOCAL_WORKSPACE, base_name)
# 运行 Pangolin
cmd = ["pangolin", local_vcf_path, LOCAL_FASTA, LOCAL_DB, local_out_prefix]
try:
subprocess.run(cmd, check=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
# 把结果传回 Drive
for f in os.listdir(LOCAL_WORKSPACE):
if f.startswith(base_name) and f != vcf_name:
shutil.copy2(os.path.join(LOCAL_WORKSPACE, f), os.path.join(OUTPUT_DIR, f))
# 生成 .done 标记
with open(expected_output_marker, 'w') as f:
f.write("Completed.")
print(f"✅ [完成] {vcf_name} 结果已回传。")
except subprocess.CalledProcessError as e:
print(f"❌ [错误] 处理 {vcf_name} 时发生异常：\n{e.stderr.decode('utf-8', errors='ignore')}")
finally:
# 清理本地 VCF 和结果，防爆盘
for f in os.listdir(LOCAL_WORKSPACE):
if f.startswith(base_name):
os.remove(os.path.join(LOCAL_WORKSPACE, f))
print("\n🎉 调度脚本运行结束。")
- 在CoLab中，运行下面的代码进行Vcf切片测试，测试的是区域
17:7661779-7687702
import os
import pyfastx
# 路径设置
DRIVE_BASE = '/content/drive/MyDrive/Pangolin_Project'
INPUT_DIR = os.path.join(DRIVE_BASE, 'inputs')
LOCAL_WORKSPACE = '/content/local_workspace'
FASTA_PATH = os.path.join(LOCAL_WORKSPACE, 'Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz')
CHROM = "17"
START = 7661779
END = 7687702
POSITIONS_PER_CHUNK = 10000
fasta = pyfastx.Fasta(FASTA_PATH)
ref_seq = fasta[CHROM][START-1:END].seq.upper()
bases = ['A', 'C', 'G', 'T']
chunk_index = 1
current_vars = []
def write_csv(vars_list, chunk_idx):
"""将突变列表写入极简的 CSV 文件"""
csv_name = f"target_region_chunk_{chunk_idx:03d}.csv"
csv_path = os.path.join(INPUT_DIR, csv_name)
with open(csv_path, 'w') as f:
# Pangolin 默认识别的四个表头
f.write("CHROM,POS,REF,ALT\n")
for var in vars_list:
f.write(f"{var[0]},{var[1]},{var[2]},{var[3]}\n")
print(f"✅ 生成任务切片: {csv_name} (包含 {len(vars_list)} 个 SNV)")
print("\n--- 开始生成 CSV 切片 ---")
for i, ref_base in enumerate(ref_seq):
pos = START + i
if ref_base not in bases:
continue
for alt_base in bases:
if alt_base != ref_base:
current_vars.append((CHROM, pos, ref_base, alt_base))
if len(current_vars) >= POSITIONS_PER_CHUNK * 3:
write_csv(current_vars, chunk_index)
chunk_index += 1
current_vars = []
if current_vars:
write_csv(current_vars, chunk_index)
print(f"\n🎉 CSV 切片生成完毕！")
- 在CoLab中，重新执行3的代码，进行注释测试
理论上，完全跑通后，就可以自动运行。但是，CoLab是有时间限制的，普通用户似乎是1小时，而Pro用户似乎是24小时。在超出这个运行时长后，重新打开时，就要重新授权Google Drive权限，重新安装软件等等（相当于启动了一台新的服务器）。
因此，上述的代码是可以识别哪些切片是已经处理过的，不会重复处理。但是仍然需要中断后，手动重新执行。
后记
本篇是一个Pangolin替代SpliceAI的使用建议，以及白嫖Google大善人的CoLab服务的初探。当前我的CoLab主要是作为往Google Drive里用命令行下载进数据的接口，以及推送第三方数据到我的抱抱脸的接口，这样节省了本地下载以及上传的大量流量费用。