生信组学分析模块 · RNA-seq 三步脚本真实运行记录

• 验真对象：content/03_生信组学分析/scripts/ 下的 run_rnaseq.sh、deseq2_de.R、enrichment_plot.R（材料自带，原文声明「仅语法检查，未端到端运行」）
• 运行环境：Ubuntu 云主机，8 核 / 15 G 内存 / 79 G 磁盘（可用 29 G），无 docker，micromamba 2.9.0，R 4.4.1 + Bioconductor 3.20
• 运行日期：2026-09-23
• 记录人：小七（子任务执行）
• 工作目录：verify/rnaseq/（真实运行在 project/；脚本原文件备份在 backup/；运行日志在 log/）

0. 结论速览（先看这段）

产物清单（均在本目录）： - run_log.md（本文件） - env_versions.txt（工具与 R 包实测版本、镜像配置、与 environment.yaml 的差异） - diff_scripts.md（三个脚本的完整 diff + 未改动文件的原因） - input_manifest.txt（数据来源、下载命令、全部输入文件 md5、样本分组映射） - backup/（改动前的原文件：3 个脚本 + environment.yaml + README.md） - log/（run_rnaseq.log、deseq2.log、enrichment.log、env_create.log、env_add_rpkg.log、fc_segfault_tmp/） - project/（完整运行目录：00_rawdata/、reference/、results/01–07、samples.tsv、dl/） - rpkg_cache/（手动补装的 3 个 Bioconductor 注释数据包） - scratch/（索引复测、featureCounts 复现性实验等临时产物）

1. 环境

1.1 环境与频道

export MAMBA_ROOT_PREFIX=/root/micromamba_root
# 环境已建好：/root/micromamba_root/envs/rnaseq（micromamba 2.9.0，频道 conda-forge + bioconda）

本次真实执行的两条安装命令（取自 $PREFIX/conda-meta/history）：

micromamba install -y -n rnaseq -c conda-forge -c bioconda \
    bioconductor-clusterprofiler bioconductor-org.hs.eg.db bioconductor-enrichplot \
    bioconductor-genomeinfodbdata bioconductor-go.db bioconductor-dose bioconductor-gosemsim
micromamba install -y -n rnaseq -c conda-forge r-ggridges   # enrichplot 画 ridgeplot 用

基础环境（fastqc / fastp / hisat2 / samtools / subread / multiqc / R 4.4.1）在本次任务之前已由同一个 rnaseq 环境装入（日志 log/env_create.log），本次没有重建；下表版本是本次实测输出。

1.2 实测版本（完整见 env_versions.txt）

FastQC v0.12.1        fastp 1.3.7           hisat2 2.2.3
samtools 1.24         featureCounts 2.1.1   multiqc 1.35
R 4.4.1 (2024-06-14)  → Bioconductor 3.20
DESeq2 1.46.0   apeglm 1.28.0   clusterProfiler 4.14.0   enrichplot 1.26.1
org.Hs.eg.db 3.20.0   GO.db 3.20.0   GenomeInfoDbData 1.2.13   fgsea 1.32.2

1.3 装 R 包时的坑（已解决，供后来者省时间）

1. bioconductor.org 直连约 5 KB/s（装 org.Hs.eg.db 的 93 MB 要跑 5 小时，上一轮任务因此被叫停）。 实测西湖镜像 https://mirrors.westlake.edu.cn/bioconductor 可达 10–24 MB/s；清华的 bioconductor 镜像返回 403，不可用。 配置写进 ~/.Rprofile 与 $PREFIX/lib/R/etc/Rprofile.site： r options(BioC_mirror = "https://mirrors.westlake.edu.cn/bioconductor") options(repos = c(CRAN = "https://mirrors.westlake.edu.cn/CRAN")) options(timeout = 3600)
2. conda 的 post-link 脚本下载注释数据包不认 BioC_mirror。它读的是 $PREFIX/share/bioconductor-data-packages/dataURLs.json，里面写死了 bioconductor.org。 改写为西湖镜像后（备份 dataURLs.json.orig_backup）依然失败：post-link 的 R CMD INSTALL 报 there is no package called 'GenomeInfoDbData'——依赖装早了，注释包根本装不进去。 最终解法是手动下载后按依赖顺序装（顺序不能颠倒）： bash R CMD INSTALL --library=$PREFIX/lib/R/library GenomeInfoDbData_1.2.13.tar.gz R CMD INSTALL --library=$PREFIX/lib/R/library GO.db_3.20.0.tar.gz R CMD INSTALL --library=$PREFIX/lib/R/library org.Hs.eg.db_3.20.0.tar.gz 装完验证：org.Hs.eg.db 有 40,839 个 ENSEMBL key，mapIds() 映射正常。
3. enrichplot 的 ridgeplot 依赖 r-ggridges，conda 的 enrichplot 没带上，需单独装。

对教程的建议：环境章节里「BiocManager::install(...)」在大陆网络下大概率超时， 建议正文直接给西湖镜像的 options(BioC_mirror=...) 两行（已实测有效），这是读者最容易卡住的一步。

2. 数据来源（真实可复现）

2.1 下载命令

cd verify/rnaseq/project
mkdir -p dl && cd dl
curl -L -C - -o griffith-data.tar.gz \
  http://data.biostarhandbook.com/rnaseq/projects/griffith/griffith-data.tar.gz
# 131,302,475 字节  md5 = 6396664c0596bd4fad89926f99ad65d5
# 实测 1200 s 超时中断于约 99 MB，用 curl -C - 续传后完成（Gitee 上的同名镜像 403，不可用）
tar xzf griffith-data.tar.gz     # 得到 reads/ 与 refs/

内容：6 个样本的双端 fastq（HBR_1..3 / UHR_1..3，102 bp）+ 22.fa / 22.gtf（人 chr22）+ ERCC92.fa / ERCC92.gtf。

2.2 参考文件拼装

cd refs
cat 22.fa  ERCC92.fa  > <project>/reference/genome.fa        # 93 条序列，50 MB
cat 22.gtf ERCC92.gtf > <project>/reference/annotation.gtf  # 56,472 行

HISAT2 索引由脚本自动构建（本次为验证该步骤，单独复测了一次）：

hisat2-build -p 8 reference/genome.fa reference/hisat2_index/genome
# real 0m49.649s → 索引 63 MB（genome.1.ht2 … genome.8.ht2）

2.3 样本表与分组映射

HBR→ctrl_*（control）、UHR→treat_*（treated）；fastq 统一改成 <sample>_1/_2.fastq.gz。

sample_id   group
ctrl_1  control      ctrl_2 control      ctrl_3 control
treat_1 treated      treat_2    treated      treat_3    treated

说明：原脚本要求样本表可带 # 注释行，所以最初的表头写成了 # sample_id group—— 这恰好踩中了 deseq2_de.R 的坑（见 §4-⑤），已一并修掉。

全部输入 md5 见 input_manifest.txt。

3. 每步真实命令 · 耗时 · 关键输出

3.1 主流程 run_rnaseq.sh

cd "/Coze/Drive/2026C-PsyBit（1）/生信教程网站/verify/rnaseq/project"
export PATH=/root/micromamba_root/envs/rnaseq/bin:$PATH
{ date '+===START %H:%M:%S==='; time bash \
  "/Coze/.../content/03_生信组学分析/scripts/run_rnaseq.sh"; echo "EXIT=$?"; } \
  > ../log/run_rnaseq.log 2>&1

结果：EXIT=0，real 3m21.448s。脚本自带时间戳，逐步真实耗时：

关键输出 1：比对率（results/03_bam/*.hisat2.log）

ctrl_1  117291 reads; 116714 (99.51%) aligned concordantly exactly 1 time   99.95% overall
ctrl_2  143086 reads; 142353 (99.49%) aligned concordantly exactly 1 time   99.95% overall
ctrl_3  128156 reads; 127516 (99.50%) aligned concordantly exactly 1 time   99.95% overall
treat_1 222259 reads; 218428 (98.28%) aligned concordantly exactly 1 time   99.94% overall
treat_2 158607 reads; 155863 (98.27%) aligned concordantly exactly 1 time   99.94% overall
treat_3 181463 reads; 178411 (98.32%) aligned concordantly exactly 1 time   99.94% overall

（比对率高是因为数据本身就是 chr22 富集的；treat 组「exactly 1 time」比例略低，符合 UHR 样本更杂的常识。）

关键输出 2：featureCounts 计数（results/04_counts/gene_counts.txt.summary）

Status                   ctrl_1  ctrl_2  ctrl_3  treat_1  treat_2  treat_3
Assigned                  94415  115858  102953   181540   120842   148710
Unassigned_NoFeatures     19145   22493   20962    31816    31714    25553
Unassigned_Ambiguity       3199    4049    3648     5535     3588     4482
Unassigned_MultiMapping    1308    1667    1439     7287     5563     6017
Unassigned_Unmapped          43      52      53      101       71       86

关键输出 3：干净计数矩阵（results/04_counts/counts_matrix.tsv，1,464 行 = 1,463 基因 + 表头）

Geneid  ctrl_1  ctrl_2  ctrl_3  treat_1 treat_2 treat_3
ENSG00000277248.1   0   0   0   0   0   0
ENSG00000274237.1   0   0   0   0   0   0
...
ERCC-00171  905 1227    1037    1842    1284    1602

（表尾出现 ERCC 基因，说明 spike-in 拼接生效；列名是样本名而不是绝对路径，这一条是修出来的，见 §4-④。）

完成目录文件数：01_fastqc/24、02_trimmed/24、03_bam/24、04_counts/3、05_qc_report/2。

3.2 deseq2_de.R

cd <project> && export PATH=/root/micromamba_root/envs/rnaseq/bin:$PATH
{ date ...; time Rscript ".../deseq2_de.R"; echo "EXIT=$?"; } > ../log/deseq2.log 2>&1

结果：EXIT=0，real 0m11.227s。关键输出（log/deseq2.log 原文）：

样本数：6；基因数：1463
分组：control vs treated
过滤后保留基因数：415
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
using 'apeglm' for LFC shrinkage...
显著差异基因：188 个（上调 89，下调 99）
全部完成，结果见：.../results/06_DE

产物 results/06_DE/（8 个文件）：DE_results_all.tsv（415 行）、DE_genes_significant.tsv（188 行，列 baseMean / log2FoldChange / lfcSE / pvalue / padj / gene）、normalized_counts.tsv、volcano.pdf+.png、 PCA.pdf、sample_distance_heatmap.pdf、top50_heatmap.pdf。p.adjust 最小的几个基因 log2FoldChange 在 +4.03 / +3.36 / +4.19 / +1.97 / −2.94 / +1.14，方向有正有负，符合两组真实样本的差异。

3.3 enrichment_plot.R

cd <project> && export PATH=/root/micromamba_root/envs/rnaseq/bin:$PATH
{ date ...; time Rscript ".../enrichment_plot.R"; echo "EXIT=$?"; } > ../log/enrichment.log 2>&1

结果：EXIT=0，real 3m44.655s（其中约 2 min 是两次 KEGG 各 60 s 超时）。关键输出：

上调基因 89，下调基因 99
'select()' returned 1:1 mapping between keys and columns
  In bitr(...) : 17.98% of input gene IDs are fail to map...
Reading KEGG annotation online: "https://rest.kegg.jp/link/hsa/pathway"...
  [跳过 KEGG] enrichKEGG 失败：cannot open the connection to 'https://rest.kegg.jp/link/hsa/pathway'
  2: In file(con, "r") : URL 'https://rest.kegg.jp/link/hsa/pathway': Timeout of 60 seconds was reached
...（下调基因同样跳过一次）
using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
preparing geneSet collections... GSEA analysis... leading edge analysis... done...
  Warning: There are ties in the preranked stats (0.89% of the list).
Picking joint bandwidth of 0.835
富集分析完成，结果见：.../results/07_enrichment

产物 results/07_enrichment/（10 个文件，条数=表行数）：

GO_BP_down.tsv 与 GO_MF_down.tsv 未生成：脚本对空结果执行 next（不落盘），属正常行为， 不是报错。GSEA 结果方向很合理——UHR（通用人类参考 RNA）相对 HBR（脑参考 RNA）在免疫相关通路上调、 神经突触相关通路下调，与两组组织的生物学背景一致。

4. 踩坑与解决（含 7 处脚本缺陷）

按「读者照着教程跑会不会踩」排序，①–⑦ 为脚本缺陷（已修），⑧ 为本机网络限制（未跑通，已在脚本内做优雅降级）。

① samtools sort 内存按「每线程」算，SORT_MEM=4G × 8 线程直接把机器打死 现象：[E::bam_sort_core] couldn't allocate memory for bam_mem，主流程中断。 原因：samtools sort -m 是单线程排序缓存，总占用 ≈ SORT_MEM × THREADS = 32 GB > 本机 15 GB。 修复：SORT_MEM="1G"（并加注释说明这个乘法关系）。16 G 内存的机器上这是必踩项。

② featureCounts --countReadPairs 多线程结果不可复现，且偶发段错误（最严重的一处） 现象：6 个 BAM 一起计数时 Segmentation fault (core dumped)，并留下 9 个 temp-core-*.tmp 残留文件 （已移到 log/fc_segfault_tmp/）。 复现实验（scratch/ 下）：同一份 BAM、同一命令行，-T 8 -p --countReadPairs 连跑 8 次得到 8 个不同的 md5；而 - -T 1 -p --countReadPairs：重复运行结果完全一致（确定性）； - -T 8 -p（去掉 --countReadPairs）：确定性； - 换 subread 2.0.6 建独立环境复测：-T 8 --countReadPairs 同样不可复现；-T 1 下 2.0.6 与 2.1.1 结果完全一致。 结论：竞态出在「按 fragment（read pair）计数」分支的多线程实现里，与版本无关；对差异表达分析来说 「同一份数据两次跑出不同计数矩阵」是致命的。 修复：新增 FC_THREADS=1 并让 featureCounts 用它。附带实测：小数据上 -T 1（16 s）比 -T 8（21 s）还快。

③ 项目路径含全角括号 → grep 把计数文件当二进制，矩阵为空 现象：主流程 EXIT=0 但 counts_matrix.tsv 是空文件，只打印 Binary file matches。 原因：featureCounts 会把输入 BAM 的绝对路径写进注释行；本项目路径含全角 （1），而本机有效 locale 是 C（en_US.UTF-8 并未真正安装），GNU grep 检测到非 ASCII 字节就判定为二进制文件， grep -v '^#' 因而什么都不输出。 修复：改用 grep -a -v '^#'。（顺带说明：同一个坑也会命中任何含中文路径的机器——国内用户几乎必踩。）

**④ 计数矩阵列名是绝对路径，deseq2_de.R 直接 stop() ** 现象：即使矩阵有内容，列名也是 /Coze/.../ctrl_1.sorted.bam，而 deseq2_de.R 用 intersect(colnames(counts), rownames(coldata)) 对齐样本，交集为 0 → 直接停止。 修复：在生成矩阵时用 awk 把表头改回样本名（sub(/.*\//,"",n); sub(/\.sorted\.bam$/,"",n)）。

⑤ 样本表表头写成 # sample_id 会被 read.delim(comment.char="#") 整行吃掉 现象：deseq2_de.R 报 replacement has 0 rows, data has 5——报错信息完全指不到真正原因。 原因：原脚本用 read.delim(coldata_file, comment.char = "#") 读样本表，而注释字符 # 恰好是表头首字符。 修复：脚本改为 readLines() + 只对首行剥 #（兼容两种写法）；同时把运行目录的 samples.tsv 表头 改成不带 # 的 sample_id。这一处值得写进教程正文：样本表自带的「# 开头为注释」约定和 R 的 comment.char 默认值撞车，属于设计层面的坑。

⑥ vst() 的 nsub=1000：小数据必然报 less than 'nsub' rows 现象：过滤后只有 415 个基因，vst(dds, blind=FALSE) 报 less than 'nsub' rows，PCA 与热图全出不来。 修复：nrow(dds) >= 1000 时用 vst()，否则回退 varianceStabilizingTransformation()（用全部基因）。 凡是用 chr22 / 子集 / 小 panel 做教学示例的，都会踩这一条。

⑦ enrichment_plot.R 里筛选了一个不存在的列（原脚本必然报错） 现象：原代码 all_res[!is.na(all_res$padj) & !is.na(all_res$entrez), c("gene","log2FoldChange")]—— DE_results_all.tsv 里没有 entrez 列，按原样运行 GSEA 段必然失败。 修复：只按 padj 过滤，Entrez 映射交给下面的 bitr() 现场做（与 deg 段保持一致的写法）。 另补：Ensembl ID 带版本号（ENSG...7）而 org.Hs.eg.db 的 key 不带版本号，脚本已加 sub("\\.[0-9]+$", "", id)；不加的话 bitr 会全量 NA、富集结果为空（修后仍未映射比例 17.98% / 12.12%， 因为 chr22 上大量 lncRNA 本就没有 Entrez ID，属正常）。

⑧ KEGG 富集：本机网络拿不到数据（未跑通，如实记录） 现象：enrichKEGG() 打印 Reading KEGG annotation online: "https://rest.kegg.jp/link/hsa/pathway"... 后 一直不回，直到超时才报 Timeout of 60 seconds was reached。（单独探针脚本 scratch/kegg_probe.R 确认能建连但无数据返回。） 处理：没有改动分析逻辑，只给 enrichKEGG 套了 tryCatch + options(timeout = 60)，失败时打印 [跳过 KEGG] ... 并继续跑完 GO 与 GSEA。这样：联网环境能出 KEGG 结果，离线环境也不会整段脚本崩掉、 更不会卡满默认超时（默认 60 s 已改；若读者的 options(timeout) 未设则会等更久）。

关于 environment.yaml 未回填版本：它的头部已声明「写的是撰写时的版本号，仅作参考，报找不到就删掉 =版本号」，属自述的软 pin；而清单里 star / nextflow / salmon / rsem 等大量包本次没装、没验证， 按实测版本回填会制造「未验证的断言」。实测版本差异（如 fastp 0.23.4 → 1.3.7、subread 2.0.6 → 2.1.1、 r-base 4.3 → 4.4.1）列在 env_versions.txt，是否回填由内容团队决定。

5. 给主 agent / 内容团队的结论与建议

已经验证的事实（可以对外说「实测跑通」的部分） 1. 三个脚本的算法与流程正确：6 个真实双端样本走完 fastp → HISAT2 → samtools → featureCounts → DESeq2（含 apeglm 收缩）→ GO/GSEA，全部 EXIT=0，比对率 99.94–99.95%，产出 188 个显著差异基因、 GO 富集表与 GSEA ridgeplot，生物学方向合理（免疫上调 / 神经突触下调）。 2. 全流程在 8 核 / 15 G 的单机上，用 chr22+ERCC92 子集的数据量（约 95 万对 reads）8 分钟内跑完 （主流程 3m21s + 差异 11s + 富集 3m45s）。教程里 README 写的「约 20–60 分钟」对全基因组数据是合理的， 本次小数据更快，不矛盾。 3. KEGG 部分只在联网时可用；本机 rest.kegg.jp 不可达，这一环节本次没有验证成功。 教程正文若给出 KEGG 结果图，需要说明「需联网、离线请只跑 GO/GSEA」，不能声称「已在离线环境验证」。

必须改的地方（不改读者会真卡住） - 三个脚本的 7 处缺陷已修好并原地更新到 content/03_生信组学分析/scripts/，原始文件在 verify/rnaseq/backup/，逐条 diff 在 verify/rnaseq/diff_scripts.md。建议正文与脚本保持同版本： 尤其 §4-②（featureCounts 多线程结果不可复现+段错误）、§4-⑤（# sample_id 表头被注释吞掉）这两条， 前者会让读者得到「跑两次结果不一样」的科学性事故，后者报错信息完全指不到原因，值得在正文里点明。 - README.md 的「五、常见问题」表已覆盖「计数矩阵列名与 sample_id 不一致」「enrichKEGG 需联网」， 建议再补 3 行（措辞可直接用）： | 现象 | 原因与处理 | | --- | --- | | samtools sort 报 couldn't allocate memory for bam_mem | -m 是单线程缓存，总内存 ≈ SORT_MEM × THREADS；16 G 机器请把 SORT_MEM 调到 1G | | featureCounts 结果每次跑都不一样 / 段错误 | --countReadPairs 在多线程下有竞态；把 FC_THREADS 设为 1（小数据上并不更慢） | | 计数矩阵是空的，日志里只有 Binary file matches | 项目路径含中文/全角字符且系统是 C locale，grep 把文件判成二进制；用 grep -a | - README 里 Rscript deseq2_de.R 的调用方式建议写成「先 cd 到项目根目录再跑」，与新的 WORKDIR="${WORKDIR:-$(pwd)}" 默认值配套；以及 §七 建议补一句「大陆网络装 Bioconductor 包请先设 options(BioC_mirror="https://mirrors.westlake.edu.cn/bioconductor")」（实测 5 KB/s → 10–24 MB/s）。

需要内容团队决定、我没有擅自改的 - environment.yaml 是否按实测版本回填 pin（我的判断：不回填更安全，因为它已经声明是参考版本， 且清单里多数包本次未验证）。 - 教程示例数据是否就用这套 griffith（chr22+ERCC92）：它的好处是在小机器上分钟级跑完、且是真实测序 数据（不是模拟），适合给读者做端到端练习；代价是基因数只有 1,463、差异基因 188，富集条目偏少。 若要更「好看」的富集图，需要换更大规模数据，否则 GO_BP_down / GO_MF_down 这类空结果会经常出现 （脚本会静默 next 不落盘，读者可能误以为脚本坏了——建议正文说明这一点）。

附：本记录对应的关键命令一览

# 0) 环境（已建好）
export MAMBA_ROOT_PREFIX=/root/micromamba_root
export PATH=/root/micromamba_root/envs/rnaseq/bin:$PATH

# 1) 数据
cd verify/rnaseq/project && mkdir -p dl && cd dl
curl -L -C - -o griffith-data.tar.gz \
  http://data.biostarhandbook.com/rnaseq/projects/griffith/griffith-data.tar.gz
tar xzf griffith-data.tar.gz
cat refs/22.fa  refs/ERCC92.fa  > ../reference/genome.fa
cat refs/22.gtf refs/ERCC92.gtf > ../reference/annotation.gtf

# 2) 主流程（cd 到项目根目录，WORKDIR 默认取当前目录）
cd verify/rnaseq/project
time bash "<repo>/content/03_生信组学分析/scripts/run_rnaseq.sh"      # EXIT=0, 3m21s

# 3) 差异表达
time Rscript "<repo>/content/03_生信组学分析/scripts/deseq2_de.R"      # EXIT=0, 11s

# 4) 富集（GO / KEGG / GSEA；KEGG 需联网）
time Rscript "<repo>/content/03_生信组学分析/scripts/enrichment_plot.R" # EXIT=0, 3m45s