# 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 需要联网
```

输出目录：

```
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"` |

## 六、扩展方向

- 想省去比对，用伪比对：`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()`**，论文方法学部分要写清楚版本。