数据预处理：原始测序数据质量控制、比对与基因表达矩阵生成

从测序仪拿到原始fastq文件的那一刻，我心里清楚——这堆ATCG看着漂亮，但里面掺了多少接头、多少低质量碱基，才是决定分析成败的关键。本章的任务就是把这堆“生数据”变成细胞里的基因表达计数矩阵，说白了，就是做三件事：质量评估、数据清洗、比对定量。每一步踩过的坑，我都会用实际项目中的例子讲清楚。

1. 原始数据质量评估——FastQC是最好的体检报告

拿到fastq第一件事，二话不说先跑FastQC。我习惯用它的batch模式一次跑完所有样本，省事。

fastqc -t 8 -o qc_report/ *.fastq.gz

跑完会生成一堆html报告。我最关注这几个指标：

• Per base sequence quality：看每个碱基位置的Q值箱线图。如果尾部（比如后100bp）掉到Q20以下，就得考虑截短。
• Per sequence GC content：理论上是正态分布，如果出现双峰，说明可能有污染或建库扩增偏好。我在一个空间转录组项目里见过GC含量在40%和60%各有一峰，后来发现是残留的rRNA没去除干净。
• Overrepresented sequences：这里会列出高频序列。最常见的是接头和polyA尾巴。我记得有一次发现某个样本里20%的read都是ployA，一看注释文件，原来是mRNA capture beads上的polyT序列没被去除，导致跟组织polyA尾巴互补后大量复制。

2. 数据清洗——Trimmomatic / cutadapt

评估完就该动手洗数据了。我常用Trimmomatic，参数看着多，但其实就四大块：剪接头、滑窗截短、截头去尾、过滤短read。

java -jar trimmomatic-0.39.jar PE \
  sample_R1.fastq.gz sample_R2.fastq.gz \
  sample_R1_clean.fastq.gz sample_R1_unpaired.fastq.gz \
  sample_R2_clean.fastq.gz sample_R2_unpaired.fastq.gz \
  ILLUMINACLIP:adapters.fa:2:30:10 \
  LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

讲讲几个容易踩坑的参数：

• SLIDINGWINDOW:4:15 —— 窗口4个碱基，平均质量低于15就截断。我曾经在项目里用SLIDINGWINDOW:4:20，结果大量read后半段全部被切掉，比对率反而下降了。后来改成15，数据利用率提高一大截。
• MINLEN:36 —— 截完后read长度不足36bp的直接丢掉。这个值要根据你后面比对工具的要求来设。比如STAR一般要求最短read长度>30bp。
• LEADING/TRAILING —— 切除5'和3'端低质量碱基。我一般设3，太严容易把有效序列误切。

3. 比对与基因表达定量

数据洗干净了，就该把read贴到参考基因组上了。空间转录组最常用的比对工具是STAR，搭配featureCounts或STAR自带的--quantMode GeneCounts定量。如果是10x Visium数据，官方工具spaceranger已经封装好了，但本质还是STAR+wk定量。

我这边以STAR为例：

STAR --genomeDir /path/to/genome_index \
     --readFilesIn sample_R1_clean.fastq.gz sample_R2_clean.fastq.gz \
     --readFilesCommand zcat \
     --outSAMtype BAM SortedByCoordinate \
     --outBAMsortingThreadN 8 \
     --quantMode GeneCounts \
     --outFileNamePrefix sample_

运行完会生成一个 sample_ReadsPerGene.out.tab 文件，里面四列：GeneID、unstranded、first strand、second strand。对于大多数空间转录组实验（如Visium用的是polyA捕获），我们取second strand那列作为基因表达计数。

如果要用featureCounts重新定量：

featureCounts -T 8 -t exon -g gene_id \
  -a annotation.gtf \
  -o counts.txt \
  sample_Aligned.sortedByCoord.out.bam

最终我们得到的是一个矩阵：行为基因、列为细胞/spot，数值是UMI count。但这时候的矩阵还是原始的，不能直接用——还需要QC过滤、归一化。嗯，那些是后面几章的事。

4. 质量控制报告解读与决策

整个预处理过程，其实每个环节都在做选择。我整理了一张表格，方便你在实际项目中快速决策：

5. 空间转录组特有注意事项

普通单细胞转录组的预处理经验搬到空间分析里，有几个特别容易忽略的地方：

• 空间条形码的完整性：10x Visium的R1包含16bp条形码+12bp UMI，R2才是cDNA序列。清洗和比对时必须保证R1的条形码区域不被截断或错误过滤。
• 双端比对的策略：我建议直接用STAR的ENCODE模式（--outSAMstrandField intronMotif），可以同时捕获外显子和内含子区域，对于组织切片中未剪接的前体mRNA也有足够表达计数。
• 参考基因组选择：空间转录组通常用完整基因组（包括未注释的lncRNA和反义链），不要只用蛋白编码基因。我吃过一次亏——只用了Ensembl的“protein_coding” biotype，结果一批基因表达谱严重缺损。

预处理做到这一步，你就拿到了干净且可比较的基因表达矩阵。接下来就可以根据检测到的基因数、线粒体比例、UMI分布等指标过滤掉不合格的spot/细胞。不过那部分内容放在下一章专门讲解，这里先帮你把基础打牢。