📚 课程目录第2章原始数据处理基础

原始数据处理基础：从BCL到表达矩阵

咱们直接进入正题。这一章咱们聊聊单细胞测序最底层的活儿——怎么把测序仪吐出来的BCL文件，折腾成咱们能用的表达矩阵。

很多人一上来就对着表达矩阵做质控、降维、聚类，但说实话，基础数据处理要是埋了雷，后面所有分析都是虚假繁荣。我一直强调：上游流程比下游推演更考验工程能力。今天就把三个主流工具——Cell Ranger、STARsolo、alevin——掰开揉碎讲一遍。

BCL → FASTQ：谁也绕不开的第一步

测序仪产出的是BCL文件（bcl格式），这是二进制原始数据。想读它得先转成FASTQ。这一步通常用Illumina官方的bcl2fastq或更新的bcl-convert。

我在一个大规模10X项目里见过有人直接跳过了抽检环节，结果跑完Cell Ranger发现barcode错误率奇高——后来查出来是sample sheet写错了。唉，白算三天。

关键参数：

--sample-sheet：决定哪些样本、哪些lane合并。格式错一个逗号直接报错。
--use-bases-mask：定义哪几个碱基是barcode、哪几个是UMI、哪几个是转录本序列。10X v2/v3的mask不同。
--no-lane-splitting：建议打开，避免一个样本拆成多个lane目录，后面整理麻烦。

曾经有次我处理完才发现某个样本的lane1和lane2深度差了三倍——因为没加--no-lane-splitting，每个lane单独输出，后面的合并脚本还得自己写。记住，单细胞测序的lane差异通常无需拆分。

命令行示例：

bcl2fastq --runfolder /data/raw/run_20240801 \
  --sample-sheet /data/meta/SampleSheet.csv \
  --use-bases-mask Y26,I8,Y91 \
  --no-lane-splitting \
  --output-dir /data/fastq

三大流程对比：选错工具，事倍功半

工具	核心策略	最快场景	典型坑
Cell Ranger	内嵌STAR + 自定义barcode校正	10X标准数据	非10X数据强改参数容易崩
STARsolo	利用STAR比对 + 内置umi/barcode处理	需要比对灵活性的场景	多sample合并要靠脚本
alevin	基于salmon的轻量化比对 + 计数	大样本、内存有限、速度优先	双端数据对齐模式需注意

我个人的选择策略很简单：如果跑10X官方的数据，无脑Cell Ranger；如果用的是非10X的商业试剂盒或自建protocol，STARsolo更灵活；如果样本量上百个、服务器内存还吃紧，alevin能省一半时间。

Cell Ranger：开箱即用，但别太信任默认

Cell Ranger的流程是串起来的：从FASTQ到gene expression matrix一条命令搞定。但许多人不知道它内部做了什么。我用一个项目举例：

cellranger count --id=Sample1 \
  --transcriptome=/ref/refdata-gex-GRCh38-2024-A \
  --fastqs=/data/fastq \
  --sample=Sample1 \
  --expect-cells=5000 \
  --localcores=16 \
  --localmem=64

--expect-cells这个参数影响barcode过滤的阈值。我遇到过有用户设了10000，但实际只有2000个细胞，结果大部分barcode被当作background保留，后面去doublet都救不回来。老老实实按经验估一个值，或跑完cellranger count后看那个filtered_feature_bc_matrix里的barcode数，再微调参数重跑。

千万不要直接用raw_feature_bc_matrix做下游分析！那里面包含了所有没有被过滤掉的背景barcode。Cell Ranger的过滤逻辑已经去掉了一大批低质量的barcode，但还不够。我们后面章节会讲二次质控。

Cell Ranger还会输出一个metrics_summary.csv，里面有几个关键指标：Fraction Reads in Cells（低于70%通常说明背景高）、Median Genes per Cell（低于500要考虑是否组织解离有问题）。

STARsolo：比Cell Ranger更通透

STARsolo本质是STAR比对器的扩展模式。它用一行STAR命令完成比对、barcode校正、UMI去重、表达计数。我最喜欢它的地方是——输出完全透明，每一步的barcode列表、UMI counts都能自己查。

命令示例：

STAR --genomeDir /ref/STAR_GRCh38 \
  --readFilesIn R1.fastq.gz R2.fastq.gz \
  --readFilesCommand zcat \
  --soloType CB_UMI_Simple \
  --soloCBwhitelist barcode_whitelist.txt \
  --soloCBstart 1 --soloCBlen 16 \
  --soloUMIstart 17 --soloUMIlen 12 \
  --outSAMtype BAM Unsorted \
  --outBAMcompression 6

说一下--soloType。默认的CB_UMI_Simple适用于大多数10X风格的barcode+UMI设计。如果是其他类型（比如带UMI只读、双端barcode），得改成CB_UMI_Complement或CB_UMI_SimpleCR。我踩过坑：一个客户用自制barcode（15bp + 10bp UMI），我没改--soloCBlen，结果STARsolo报错“barcode长度不对”。

STARsolo的输出目录结构：

Solo.out/Gene/raw：未过滤的barcode x gene矩阵（类似Cell Ranger的raw）
Solo.out/Gene/filtered：经过barcode数量阈值过滤后的矩阵
Solo.out/Gene/umi：UMI count原始表（可以用来验证去重逻辑）

alevin：轻量级的“快枪手”

alevin是salmon套件的一部分，采用伪比对（quasi-mapping）而不是完全匹配。好处是超快，坏处是对于高度相似的基因家族（如HLA、TCR）会有点模糊。不过对于常规的单细胞转录组，精度跟STAR几乎没差别。

我用alevin处理过一批384个样本的smart-seq2-like数据，10个样本并行跑，不到2小时全部完成。换成Cell Ranger大概要4小时。

典型调用：

salmon alevin -l ISR \
  -1 R1.fastq.gz \
  -2 R2.fastq.gz \
  --tgMap txp2gene.tsv \
  --chromium \
  --sketch \
  -p 16 \
  -o alevin_output

--chromium代表10X Chromium平台的barcode设计，如果是其他平台（如Drop-seq），得用--dropseq或自己写--barcodeLength和--umiLength。还有个参数--sketch，能在保证准确率的同时降低内存，我建议大样本必开。

表达矩阵的结构与关键核对

无论用哪个工具，最后得到的表达矩阵都是barcode × gene的稀疏矩阵。通常包含三个文件（Market Exchange Format，MEX）：

matrix.mtx：非零表达数据
barcodes.tsv：所有通过过滤的cell barcode
features.tsv：基因名/ID列表

我在做整合分析前一定会做的一件事：核对barcodes.tsv第一列的格式。不同工具输出的barcode格式可能不同——Cell Ranger是AAACCCAAGAAACTCT-1这样带后缀，STARsolo可能不带后缀，alevin默认也没有。后续整合时如果不统一，就全乱了。

我建议在进入下游分析前，先用一个简单脚本把barcode后缀统一加上或去掉。比如统一截掉-1后缀，或者追加-sample1等标识，避免合并不同样本时barcode冲突。

# 在R中快速统一格式（去掉后缀）
barcodes <- readLines("barcodes.tsv")
barcodes <- gsub("-1$", "", barcodes)
writeLines(barcodes, "barcodes_unified.tsv")

SVG知识框架：原始数据处理全流程

总结一下这一步的避坑清单

BCL转换阶段：检查sample sheet里barcode序列跟试剂盒是否一致（10X v2 vs v3 sequence不同）。
工具选择：不要用Cell Ranger处理非10X数据，不要用alevin处理超高深度的全长scRNA-seq。
参数调优：STARsolo的--soloCBwhitelist必须用对版本的白名单文件；Cell Ranger的--expect-cells影响很大。
输出检查：每次跑完都看一眼matrix.mtx的行数和列数是否合理，barcode数量是否在预期范围内。

我个人习惯把每个样本的barcode计数画一个简单的rank vs UMI count图（俗称“膝图”），不需要等下游，直接在终端用awk和gnuplot就能快速判断：有没有明显的拐点？拐点之前是不是有大量barcode背景？这点我们在下一章正式讲质控时还会展开。

好了，BCL到表达矩阵的路径清楚了。从原始数据到可分析的数据，这一步走稳了，后面的整合分析才能站得住。接下来我们就要拿着这些表达矩阵，开始做多样本整合和质控参数的优化了。