# scripts/ 使用说明

本目录是「生信组学分析」的 RNA-seq 流程脚本模板，bash + R 标准工具链，可直接改路径跑。

## 文件清单

| 文件 | 作用 |
|---|---|
| `run_rnaseq.sh` | 主流程：fastp 修剪 → HISAT2 比对 → samtools 排序 → featureCounts 计数 → MultiQC 汇总 |
| `deseq2_de.R` | 差异表达：DESeq2 建模、标准化、lfcShrink、火山图 / PCA / 热图 |
| `enrichment_plot.R` | 功能富集：clusterProfiler 做 GO、KEGG、GSEA，出气泡图与条形图 |
| `environment.yaml` | 全部依赖的 conda 环境清单（bioconda 渠道） |

## 一、装环境

```bash
# 推荐用 mamba，比 conda 快很多
mamba env create -f environment.yaml
conda activate rnaseq
```

如果只想跑 bash 主流程，最小依赖是：

```bash
conda create -n rnaseq -c conda-forge -c bioconda \
    fastqc fastp hisat2 samtools subread multiqc
```

R 若不用 conda 装，就走 BiocManager：

```r
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install(c("DESeq2","tximport","clusterProfiler",
                       "org.Hs.eg.db","enrichplot","ComplexHeatmap"))
install.packages(c("ggplot2","pheatmap","apeglm","tidyverse"))
```

## 二、准备数据

目录结构（照 `run_rnaseq.sh` 里的变量改即可）：

```
project/
├── 00_rawdata/          # 双端 fastq.gz，命名 <sample>_1.fastq.gz / <sample>_2.fastq.gz
├── reference/
│   ├── genome.fa        # 参考基因组（如 Ensembl GRCh38 primary assembly）
│   ├── annotation.gtf   # 同版本注释（版本必须和基因组一致！）
│   └── hisat2_index/    # 脚本自动生成，无需手动建
├── samples.tsv          # 样本表
└── results/             # 全部输出
```

样本表 `samples.tsv`（制表符分隔，两列）：

```
sample_id	group
ctrl_1	control
ctrl_2	control
treat_1	treated
treat_2	treated
```

> 参考数据下载：Ensembl 的基因组与 GTF 必须来自**同一个 release**，否则 featureCounts 会报基因对不上。

## 三、跑流程

```bash
chmod +x run_rnaseq.sh
./run_rnaseq.sh                      # 约 20–60 分钟，取决于数据量和机器
Rscript deseq2_de.R                  # 差异表达，几分钟
Rscript enrichment_plot.R            # 富集分析，enrichKEGG 需要联网
```

> 三条命令都在项目根目录下执行。R 脚本按相对路径读 `results/` 和 `samples.tsv`，换个目录跑会直接报找不到输入。

耗时参考：本项目实测用 chr22 + ERCC92 子集（6 个双端样本），主流程 3 分 21 秒、差异表达 11 秒、富集 3 分 45 秒（8 核 / 15 GB 单机）。换人类全基因组会明显更久，主要卡在比对与计数。

输出目录：

```
results/
├── 01_fastqc/          # 原始数据质控报告
├── 02_trimmed/         # 修剪后 fastq + fastp 报告
├── 03_bam/             # 排序后 BAM + flagstat
├── 04_counts/          # gene_counts.txt 与干净矩阵 counts_matrix.tsv
├── 05_qc_report/       # MultiQC 汇总 HTML（先看这个）
├── 06_DE/              # 标准化矩阵、差异表、火山图、PCA、热图
└── 07_enrichment/      # GO/KEGG/GSEA 结果表与图
```

先打开 `results/05_qc_report/multiqc_report.html` 和 `results/01_fastqc/*.html`，质控不过关就别急着看差异结果。

## 四、必改的三个地方

1. `run_rnaseq.sh` 顶部的 `WORKDIR` / `REF_DIR` / `SAMPLESHEET` 改成自己的绝对路径。
2. `run_rnaseq.sh` 里的接头序列 `ADAPTER_R1/R2`：用自己文库的接头，或直接靠 `--detect_adapter_for_pe` 自动检测。
3. `deseq2_de.R` 里的 `ref_level` 改成你的对照组名（要和 `samples.tsv` 里完全一致）。

## 五、常见问题

| 现象 | 原因与处理 |
|---|---|
| `featureCounts` 结果全是 0 | 基因组与 GTF 版本不匹配，或建库类型不对（链特异性文库需加 `-s 1/2`） |
| HISAT2 报内存不足 | 人类基因组索引约需 8–16GB，换 STAR 更吃内存，小机器改用 Salmon |
| DESeq2 报 `counts` 列名对不上 | 计数矩阵列名（BAM 文件名）与 `samples.tsv` 的 sample_id 不一致 |
| 富集结果为空 | 基因 ID 类型搞错，或基因数太少（<10 个几乎出不来结果） |
| `enrichKEGG` 报错 | 需要联网；离线改跑 `enrichGO` |
| apeglm 报错 | 没装 `apeglm`，`install.packages("apeglm")`，或把 `lfcShrink` 的 `type` 改成 `"ashr"` |
| `samtools sort` 报 `couldn't allocate memory for bam_mem` | `-m` 是**单个线程**的缓存，总内存 ≈ `SORT_MEM × THREADS`。16 GB 机器把 `SORT_MEM` 调到 `1G` |
| featureCounts 每次跑的结果都不一样 / 偶发段错误 | `--countReadPairs` 在多线程下有竞态，同一份数据能跑出不同计数。把 `FC_THREADS` 设为 1（小数据上并不更慢） |
| 计数矩阵是空的，日志里只有 `Binary file matches` | 项目路径含中文或全角字符、系统又是 C locale 时，`grep` 会把文本文件判成二进制，改用 `grep -a` |

> 表格最后三条是这套脚本在真实数据上跑的时候真踩到的坑（完整记录见站点「代码验证」页），不是推测。

## 六、扩展方向

- 想省去比对，用伪比对：`salmon index` + `salmon quant`，再用 `deseq2_de.R` 前加一段 `tximport` 汇总。
- 想要工业级整合流程：直接跑 nf-core/rnaseq（https://nf-co.re/rnaseq），`nextflow run nf-core/rnaseq -profile test,docker` 可先验环境。
- 想看剪接事件：`StringTie` 组装 + `Ballgown`，或直接用 `rMATS`。

## 七、依赖版本

`environment.yaml` 里给的是撰写时的版本，仅作参考。若创建环境时报找不到某版本，把对应的 `=版本号` 删掉即可，conda 会自动选当前可用版本。**正式分析前请保存 `conda list > env_locked.txt` 和 R 的 `sessionInfo()`**，论文方法学部分要写清楚版本。

实测时同一份 `environment.yaml` 解析到的版本普遍比清单里新（fastp 1.3.7、hisat2 2.2.3、samtools 1.24、subread 2.1.1、multiqc 1.35、R 4.4.1 + Bioconductor 3.20）。照清单建环境没问题，版本号以自己机器解析出来的为准。

### 大陆网络下装 R 包

Bioconductor 官方源国内实测只有 5 KB/s 左右，一个 `org.Hs.eg.db`（93 MB）能拖几个小时。换西湖大学镜像后实测 10–24 MB/s：

```r
options(BioC_mirror = "https://mirrors.westlake.edu.cn/bioconductor")
options(repos = c(CRAN = "https://mirrors.westlake.edu.cn/CRAN"))
options(timeout = 3600)
```

清华的 bioconductor 镜像实测返回 403，用不了。把上面的 `options()` 写进 `~/.Rprofile` 可以长期生效。

另外 `enrichKEGG` 需要联网（`rest.kegg.jp`），离线或网络受限时这一步会卡住超时；脚本已做超时跳过处理，GO 与 GSEA 不受影响。