SeekSoulTools 1.3.0 源码解读

寻因单细胞分析管线的 scRNA 与 V(D)J 流程

single-cell
python
bioinformatics
Author

Liripo

Published

June 16, 2026

寻因（SeekGene）的 SeekSoulTools（seekgenebio/seeksoultools，本文基于 v1.3.0）是一套开源单细胞管线，从 FASTQ 一路到表达矩阵和 V(D)J 组装。它整体用 Python（click）串流程，重计算的部分交给 Rust（simpleqc、search-pattern）和成熟的外部工具（STAR、featureCounts、trust4、qualimap）。

下面按数据流记一下我读源码时觉得值得留的点，分 scRNA 和 V(D)J 两块。

1 scRNA 流程

scRNA 的核心命令是 rna run，一条龙跑完 step1→step4。四个 step 分别做：

  1. step1：barcode 提取——交给 Rust 写的 simpleqc
  2. step2：比对+计数到基因——STAR、featureCounts 等
  3. step3：UMI 折叠+矩阵+cell calling
  4. step4：Seurat 降维聚类

每个 step 也可以单独跑（rna step1/step2/step3/step4）。

1.1 step1：把 barcode/UMI 提到 FASTQ 里

这步交给 Rust 写的 simpleqc。它做两件事：质控统计 + 把 read1 里的 cell barcode、UMI 校正后写到新的 R1，把 cDNA read（R2）原样输出。后续 STAR 直接拿这个 R2 比对。

输出到 Analysis/qc/ 目录下，文件大概是这些：

Analysis/qc/
├── demo_r1.fq.gz          # 校正后的 barcode+UMI
├── demo_r1.multi.fq.gz    # 多匹配的 barcode（默认丢弃）
├── demo_r2.fq.gz          # cDNA read
├── demo_r1_base_stat.csv
├── demo_r1_qual_stat.csv
├── demo_r2_base_stat.csv
├── demo_r2_qual_stat.csv
└── demo_summary.txt

simpleqc 的 read 结构用「Read Structure」描述（fgbio 规范），比如 17C12M+T 表示 17 个碱基的 Cell barcode（C）、12 个碱基的 UMI（M）、然后是模板（T）。上层用更可读的 bc:0:16,um:17:28，再由 readFormat2rs 函数翻译：

def readFormat2rs(readFormat:str):
    # "bc:0:15,um:16:25" -> "17C25M+T"
    rs = ""
    tmp = sorted(
        [p.split(":") for p in readFormat.split(",")],
        key=lambda x: int(x[1])
    )
    for (p, start, end) in tmp:
        start = int(start); end = int(end)
        if p == "bc":
            rs += f"{end-start+1}C"
        elif p == "um":
            rs += f"{end-start+1}M"
    rs += "+T"
    return rs

按 start 排序是为了保证 bc 永远在 um 前面——即便用户写反了。

1.2 step2：比对与基因计数

align 函数串起 STAR → 排序 → qualimap → featureCounts → 按名字排序，每一步都可由 steps 参数跳过：

if ("steps" not in kwargs) or (not kwargs["steps"]):
    kwargs["steps"] = ["STAR", "SortByPos", "qualimap", "FeatureCounts", "SortByName"]

整个 step2 的关键不是「怎么比对」，而是怎么从几个工具的输出里凑出 10x 风格的 QC 指标。mapping_summary 把 STAR 的 Log.final.out 和 qualimap 的 rnaseq_qc_results.txt 拆开读，重新算了一套：

指标 计算
Reads Mapped to Genome (unique + multi + too many) / Total
Reads Mapped Confidently to Genome (aligned to genes + no feature) / Total
Reads Mapped to Exonic Regions exonic / Total
Reads Mapped to Intronic Regions intronic / Total
Reads Mapped to Intergenic Regions intergenic / Total
Note

注意 “Mapped Confidently to Genome” 的算法：是 aligned to genes + no feature assigned。qualimap 里 “no feature” 指没落到任何注释区间，看似该排除，但这里把它算进 confident——因为落不到 feature 上不代表比对错，只是注释没覆盖到。比对的正确性由 STAR 的 unique mapping 决定，feature 是否注释上是另一回事。

region 参数控制计数时是否含内含子（--include-introns）。featureCounts 的链向通过 sc5p 推断：5’ 数据反向（-s 2）， otherwise 正向（-s 1）。

最后一步 SortByName 很关键：featureCounts 出来的 bam 要按 read name 排序，因为 step3 的计数逻辑是 groupby(sam_file, key=lambda x: x.qname.split("_",1)[0])——按 barcode 名分组遍历 bam。坐标排序没法高效做到这件事。

1.3 step3：UMI 折叠与 cell calling

这是整个 scRNA 最绕的一步。输入是按 barcode 分组好的 bam，输出是 raw matrix 和 filtered matrix。

1.3.1 按 (barcode, gene) 分桶 + UMI 校正

bam2table 外层按 barcode 分组，内层对每个 barcode 再按 read name（barcode_umi）二次分组：

for barcode, g in groupby(sam_file, key=lambda x: x.qname.split("_", 1)[0]):
    counts_dict, geneid_umi_dict = umi_count(g, ...)

umi_count 里每个 read 拿 featureCounts 打的 XT tag（基因 id），落到 (gene_id, umi) 桶里。这里有个多比对处理：一个 read 的 XT 可能是逗号分隔的多个基因（XT:Z:A,B），要调 multiple_alignment 用 IntervalTree 判断它到底真落哪个基因——只匹配到一个才采纳，否则丢弃。

UMI 校正用的 UMIClusterer，看实现就是 umi_tools 那套「adjacency + 连通图」聚类：汉明距离 ≤ 1 的 UMI 连边，每个连通分量里 count 最高的留下、其余并过去。umi_correct_method 支持三种（cluster/adjacency/directional），默认 adjacency。

校正完每个 (barcode, gene) 就是确定的 UMI 数，写到 counts.xls（cellID / geneID / UMINum / ReadsNum），这文件直接能转 raw matrix。

1.3.2 从 counts.xls 到稀疏矩阵

write_raw_matrix 按 barcode 分块读 counts.xls，攒成 COO 三元组再压成稀疏矩阵：

mat = coo_matrix((data, (row, col)), shape=(name_df.shape[0], n))

输出 10x 三件套：matrix.mtx.gz、features.tsv.gz、barcodes.tsv.gz。

1.3.3 cell calling：EmptyDrops 的实现

判定哪些 barcode 是真细胞，走的是 cell_identify.R 这个 R 脚本。它不是调 DropletUtils 的现成函数，而是把 10x cellranger 的 EmptyDrops 算法用 R 重写了一遍——分位数法找 top barcode、估背景 profile、多项式模拟算 p 值。

流程是：

  1. 找 top barcode：find_top_umi_barcode，按预期细胞数 recovered_cells 的 0.99 分位定基线，基线的 10% 作为 cutoff，bootstrap 100 次取均值。这一批是确定细胞的”种子”。
  2. 找 ambient candidate：UMI ≥ max(500, median*0.01)、且不是 top 的 barcode。这是”可能是细胞也可能是空滴”的待判集合，限制最多 20000 个。
  3. 估背景 profile：find_ambient_background 取 UMI 排名靠后的那批（假设是纯空滴），goodTuringProportions 做 Good-Turing 平滑得到每个基因的背景比例 ambient_prop。
  4. 算每个 candidate 的 p 值：核心假设是——空滴的表达应该服从背景的多项式分布。compute_simulate_pvalues 模拟 10000 次多项式采样，算出每个总 UMI 数 n 对应的 log-likelihood 分布；compute_ambient_candidate_pvalues 把真实 candidate 的 likelihood 跟模拟分布比，(1 + 比它小的次数)/(1 + 模拟次数) 就是 p 值。
  5. BH 校正 + cutoff：p.adjust(method="BH")，小于 pvalue（默认 0.01）的 candidate 也算细胞，跟 top 合起来输出。
Note

forceCell 模式：如果不信任自动 calling，--forceCell N 直接取 UMI 最多的前 N 个 barcode，跳过整个模拟。cell_calling 里 if kwargs["forceCell"] != None 就走这条分支。

模拟那块有个性能优化：相邻的 n（比如 n=100 和 n=101）不用各采一次样，而是在 n=100 的 curr_counts 基础上再加一个基因增量更新 likelihood（curr_loglk + log(p_j) + log(n/curr_counts[j])），省掉大量重算。

1.4 step4：Seurat

do_seurat 调 R 脚本跑标准 Seurat 流程（NormalizeData → FindVariableFeatures → ScaleData → PCA → clustering → FindMarkers），输出降维聚类结果，没什么特别的。

2 V(D)J 流程

V(D)J 入口是 run_trust4 函数。它没有从头造轮子，而是在 trust4 基础上包了一层 QC、UMI 校正、contig 重分配。整个流程：

simpleqc → enrichment_qc → fastq-extractor → subset_reads
        → filter_umis → trust4 组装 → re_assign → report

2.1 simpleqc：跟 scRNA 共用

跟 step1 一样，输出 _r1.fq.gz / _r2.fq.gz。summary_simpleqc 从 _summary.txt 读 total_num / valid_num，从 qual stat csv 算 Q30。注意：后续所有比例的分母都是 Number of Valid Reads（校正后的有效 read），不是原始 read 数。

2.2 enrichment_qc：Reads Mapped to Any V(D)J Gene

这步回答”数据里到底有没有 V(D)J”。做法是：对 bcrtcr.fa（所有 BCR/TCR 的 V/D/J/C 序列）建 STAR 索引，把 R2 比上去，数每条链的 read 数。

count_chain 有个细节值得说：

with pysam.AlignmentFile(bam) as bamfh:
    for r in bamfh:
        chain = r.reference_name[:3]   # "TRBV1" -> "TRB"
        if r.get_tag("NH")==1:         # 唯一比对
            d[chain] += 1
        else:
            tmp[r.qname][chain] += 1   # 多比对先攒着

NH tag = 1 表示唯一比对，直接计。多比对的 read（比如同时软比对到 TRBV 和 TRBJ）按 read name 攒到 dict，最后仲裁：只比对到一个链就计入；比对到多个链就比 count，count 相同判 unknown。这样避免一条 read 被多个链重复计。

summary_enrichment_qc 算链比例时按 chain 类型筛：

  • TR 数据：只看 TRA + TRB，忽略 IG 链（IGH/IGK/IGL）。TCR 文库混入 IG 是噪声，不该计入 “Any V(D)J”。
  • IG 数据：看 IGH + IGK + IGL。

而且有个硬卡：Reads Mapped to Any V(D)J Gene < 0.01 直接 logger.error，提示数据里没这个链，大概率是 chain 选错了或者文库搞错了。

2.3 fastq-extractor 与 subset_reads

trust4 的 fastq-extractor 负责：用 bcrtcr.fa 抓 R2 里命中 V(D)J 的 read，把对应的 barcode/UMI 从 R1 抽出来，输出三件套（_preassemble_bc.fa / _preassemble_umi.fa / _preassemble.fq），同时统计每个 barcode 的 UMI 数和 read 数（_bc_count.tsv）。

然后是 subset_reads：reads 数过高的 barcode 截断到 reads_limit（默认 80000）。原因是一些 barcode 可能是 doublet 或气泡污染，read 数异常高，直接组装会拖慢又影响克隆定量。做法是把每个 barcode 的 read 数封顶，多余的 read 砍掉。原始全量计数存到 _{samplename}_bc_count.tsv（带下划线前缀），截断后的覆盖 {samplename}_bc_count.tsv——这样后面算 “Fraction Reads in Cells” 时还能拿到原始数。

2.4 filter_umis：VDJ 的 UMI 校正

trust4 原版对 VDJ UMI 校正比较粗，seeksoultools 自己加了一步 filter_umis。两阶段：

1. 按阅读量过滤

每个 barcode 内，先算 UMI count 的 N50：

def calculate_n50(umi_counts: dict):
    counts = sorted(umi_counts.values(), reverse=True)
    total_sum = sum(counts)
    threshold = total_sum * 0.5
    current_sum = 0
    for count in counts:
        current_sum += count
        if current_sum >= threshold:
            return count
    return 0

N50 是组装/基因组学里的指标：把 UMI 按支持 read 数从高到低排，累计到总量 50% 时的那个 read 数。然后阈值定在 0.1 * N50，count 低于它的 UMI 丢掉。

Note

为什么用 N50 而不是固定阈值？不同 barcode 的测序深度差异大，固定阈值（比如丢掉 count<3 的）在深测序样本里会误杀真 UMI、在浅测序样本里又滤不干净。N50 是个自适应的”该 barcode 自身规模”的度量，乘 0.1 当下限。

2. adjacency 聚类校正

跟 scRNA 那步一样，UMIClusterer("adjacency") 把汉明距离 1 的 UMI 合并到 count 高的那个。校正过程写到 _{samplename}_filter.detail.xls，长这样：

barcode    original_umi  original_count  corrected_umi  corrected_count
TGATCAGATGGATTTGA  ACCGGCGTTTTC  11   ACCGGCGTTTTA  33
TGATCAGATGGATTTGA  ACCGGCGTTTTG   6   ACCGGCGTTTTA  39

最后重写出新的 barcode/UMI/R2 fastq，UMI 改了之后连 read name 都同步改掉（name_parts[1] = corrected_umi）。

2.5 trust4 组装 + re_assign

trust4 那几个子命令（trust4 组装、annotator 注释、trust-barcoderep/trust-simplerep/trust-airr 报告）seeksoultools 都用了，但加了自己的料。

Warning

trust4_wrapper 里组装那步硬编码了 -t 8，跟上层传进来的 threads 脱节：

cmd1 = (f"trust4  -t 8 -f {fa} -o {wd}/{samplename} ...")

源码里也标了注释（参考笔记：“脚本代码固定了 8 个线程，应该修改下”）。如果你跑大批量样本，这里会成为瓶颈，建议改成 -t {threads}。

annotator 加了 --readAssignment：输出 _assign.out，两列 read_name \t contig，记录每条 read 归属的 contig。这是为了下一步的 re_assign。

2.5.1 re_assign：按 UMI 重新分配 contig

这是 seeksoultools 最有自己想法的一步（utils/reassign.py）。trust4 默认按 read 分配 contig，但 VDJ 里真正该数的是 UMI——一个 UMI 对应一个原始 mRNA 分子。re_assign 把分配单位从 read 换成 UMI。

逻辑是：对每个 (barcode, UMI)，统计它下面的 read 落到各 contig 的次数，然后按 ratio 仲裁。process_assign_file 用多进程分块统计：

for barcode_umi, contig_counts in global_barcode_umi_dict.items():
    sorted_counts = sorted(contig_counts.items(), key=lambda x: x[1], reverse=True)
    if len(sorted_counts) == 1:
        final_contig_barcodes[sorted_counts[0][0]].add(barcode_umi)
    else:
        ratio = sorted_counts[1][1] / sorted_counts[0][1]
        if ratio < 0.15:
            final_contig_barcodes[sorted_counts[0][0]].add(barcode_umi)
        else:
            for contig, _ in sorted_counts:
                final_contig_barcodes[contig].add(barcode_umi)

举例，barcode bc1 的某个 UMI 下有：

contig read 数
contigA 20
contigB 1

ratio = 1/20 = 0.05 < 0.15，这个 UMI 归 contigA。

但如果：

contig read 数
contigA 20
contigB 10

ratio = 10/20 = 0.5 ≥ 0.15，说明 contigB 也有相当支持，UMI 同时归 contigA 和 contigB。

分配完，update_cdr3_file 用结果改写 cdr3.out 的第 11 列（UMI count），让每个 contig 的 UMI 数反映重新分配后的结果。这一步直接影响克隆大小的定量。

2.6 报告生成

报告那块调了 dandelion：dandelion.tl.find_clones() 按 CDR3 序列相似度在 VDJ 数据里识别克隆，生成 clone_id。Estimated Number of Cells、Number of Cells with Productive V-J Spanning Pair 这些指标都是从 dandelion 出来的 airr_rearrangement.tsv 算的（get_clone_matrix 函数）。

最后还有两个输出值得说：

_filtered_contig_annotations.tsv：用 Rust 写的 search-pattern 生成（src/search-pattern）。它基于 LEADER 序列的匹配位置，把每条 contig 切成 FWR1/CDR1/FWR2/CDR2/FWR3/CDR3/FWR4 这些功能区。FWR1 的左边界是 LEADER 匹配终点、右边界是 CDR1 匹配起点，依此类推。10x 的 all_contig_annotations 也是这个结构。

barcode rank 图：读 _bc_count.tsv，按 UMI 数降序画。关键是它把曲线分成三段画——cell / 混合 / background：

# BBBCC CCBBCBCC BBBBBBBBB
#       |         |
#    idx_first  idx_last

每个逻辑分段单独处理（参考笔记：df_bg.shape[0]==0 全是细胞、df_cells.shape[0]==0 全是背景、普通情况分”第一个 cell 前的混合段” + “中间指数窗口” + “尾段” + “background”）。窗口按 idx_s + 1000*4**n 指数扩张，保证前端密集采样、后端稀疏，画出来曲线平滑。