生信组学分析

#!/usr/bin/env bash
###############################################################################
# RNA-seq 标准分析流程模板（bash）
# 用途：从双端 FASTQ 走到基因水平计数矩阵，供下游 DESeq2 使用
# 工具链：fastp -> HISAT2 -> samtools -> featureCounts -> MultiQC
# 说明：所有路径都用变量集中管理，换数据只改 CONFIG 段即可
###############################################################################
set -euo pipefail

################################  CONFIG  #####################################
# 工作目录
WORKDIR="/path/to/project"
RAW_DIR="${WORKDIR}/00_rawdata"      # 原始 fastq.gz 放在这里
REF_DIR="${WORKDIR}/reference"        # 参考基因组 fasta 与 gtf 放在这里
OUT_DIR="${WORKDIR}/results"

# 参考文件（需自行下载，例如 Ensembl 的 Homo_sapiens.GRCh38）
GENOME_FA="${REF_DIR}/genome.fa"
GTF="${REF_DIR}/annotation.gtf"

# HISAT2 索引前缀（若不存在则自动构建）
HISAT2_INDEX="${REF_DIR}/hisat2_index/genome"

# 线程数与内存
THREADS=8
SORT_MEM="4G"

# 样本表（TSV，两列：sample_id 和分组；# 开头为注释）
SAMPLESHEET="${WORKDIR}/samples.tsv"

# fastp 参数
ADAPTER_R1="AGATCGGAAGAGCACACGTCTGAACTCCAGTCA"
ADAPTER_R2="AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT"
MIN_LEN=36          # 修剪后最短保留长度
MIN_QUAL=20         # 质量阈值
################################  END CONFIG  #################################

log() { echo "[$(date '+%Y-%m-%d %H:%M:%S')] $*"; }

mkdir -p "${OUT_DIR}"/{01_fastqc,02_trimmed,03_bam,04_counts,05_qc_report}

#------------------------------------------------------------------------------
# 0. 检查输入
#------------------------------------------------------------------------------
[[ -f "${GENOME_FA}" ]] || { echo "缺少参考基因组: ${GENOME_FA}"; exit 1; }
[[ -f "${GTF}" ]]       || { echo "缺少注释文件: ${GTF}"; exit 1; }
[[ -f "${SAMPLESHEET}" ]] || { echo "缺少样本表: ${SAMPLESHEET}"; exit 1; }

log "参考基因组建立索引（若已存在会跳过）"
samtools faidx "${GENOME_FA}"

#------------------------------------------------------------------------------
# 1. 构建 HISAT2 索引（只跑一次）
#------------------------------------------------------------------------------
if [[ ! -f "${HISAT2_INDEX}.1.ht2" ]]; then
    log "构建 HISAT2 索引"
    mkdir -p "$(dirname "${HISAT2_INDEX}")"
    hisat2-build -p "${THREADS}" "${GENOME_FA}" "${HISAT2_INDEX}"
else
    log "HISAT2 索引已存在，跳过"
fi

#------------------------------------------------------------------------------
# 2. 逐样本：质控 -> 修剪 -> 比对 -> 排序
#------------------------------------------------------------------------------
while read -r sample group; do
    [[ -z "${sample}" || "${sample}" == \#* ]] && continue
    log "处理样本 ${sample}（分组 ${group}）"

    R1="${RAW_DIR}/${sample}_1.fastq.gz"
    R2="${RAW_DIR}/${sample}_2.fastq.gz"
    [[ -f "${R1}" && -f "${R2}" ]] || { echo "缺少 ${sample} 的 fastq"; exit 1; }

    # 2.1 原始数据质控（FastQC）
    fastqc -t "${THREADS}" -o "${OUT_DIR}/01_fastqc" "${R1}" "${R2}"

    # 2.2 去接头 + 质量过滤（fastp）
    fastp \
        -i "${R1}" -I "${R2}" \
        -o "${OUT_DIR}/02_trimmed/${sample}_1.trim.fq.gz" \
        -O "${OUT_DIR}/02_trimmed/${sample}_2.trim.fq.gz" \
        --adapter_sequence "${ADAPTER_R1}" \
        --adapter_sequence_r2 "${ADAPTER_R2}" \
        --qualified_quality_phred "${MIN_QUAL}" \
        --length_required "${MIN_LEN}" \
        --detect_adapter_for_pe \
        --thread "${THREADS}" \
        --json "${OUT_DIR}/02_trimmed/${sample}.fastp.json" \
        --html "${OUT_DIR}/02_trimmed/${sample}.fastp.html"

    # 2.3 比对（HISAT2），直接转 BAM
    hisat2 -p "${THREADS}" -x "${HISAT2_INDEX}" \
        -1 "${OUT_DIR}/02_trimmed/${sample}_1.trim.fq.gz" \
        -2 "${OUT_DIR}/02_trimmed/${sample}_2.trim.fq.gz" \
        2> "${OUT_DIR}/03_bam/${sample}.hisat2.log" \
      | samtools view -@ "${THREADS}" -bS - \
      | samtools sort -@ "${THREADS}" -m "${SORT_MEM}" -o "${OUT_DIR}/03_bam/${sample}.sorted.bam" -

    samtools index -@ "${THREADS}" "${OUT_DIR}/03_bam/${sample}.sorted.bam"
    # 比对统计（后续 MultiQC 会汇总）
    samtools flagstat -@ "${THREADS}" "${OUT_DIR}/03_bam/${sample}.sorted.bam" \
        > "${OUT_DIR}/03_bam/${sample}.flagstat.txt"
done < "${SAMPLESHEET}"

#------------------------------------------------------------------------------
# 3. 基因水平计数（featureCounts）
#------------------------------------------------------------------------------
log "运行 featureCounts 汇总计数矩阵"
BAM_LIST=$(awk '!/^#/ && NF>0 {print "'"${OUT_DIR}"'/03_bam/"$1".sorted.bam"}' "${SAMPLESHEET}")

featureCounts \
    -T "${THREADS}" \
    -p --countReadPairs \
    -t exon -g gene_id \
    -a "${GTF}" \
    -o "${OUT_DIR}/04_counts/gene_counts.txt" \
    ${BAM_LIST}

# 去掉行首注释，得到干净的计数矩阵（DESeq2 可直接读）
grep -v '^#' "${OUT_DIR}/04_counts/gene_counts.txt" | cut -f1,7- \
    > "${OUT_DIR}/04_counts/counts_matrix.tsv"

log "计数矩阵完成: ${OUT_DIR}/04_counts/counts_matrix.tsv"

#------------------------------------------------------------------------------
# 4. 汇总质控报告（MultiQC）
#------------------------------------------------------------------------------
log "生成 MultiQC 汇总报告"
multiqc -f -o "${OUT_DIR}/05_qc_report" "${OUT_DIR}"

log "流程结束。下一步：Rscript deseq2_de.R"
###############################################################################
# 依赖：fastqc fastp hisat2 samtools subread multiqc
# 一次性安装（conda 环境见 environment.yaml）：
#   conda create -n rnaseq -c conda-forge -c bioconda \
#       fastqc fastp hisat2 samtools subread multiqc
# 样本表 samples.tsv 示例（制表符分隔）：
#   sample_id	group
#   ctrl_1	control
#   ctrl_2	control
#   treat_1	treated
#   treat_2	treated
# 运行：
#   chmod +x run_rnaseq.sh && ./run_rnaseq.sh
###############################################################################

#!/usr/bin/env Rscript
###############################################################################
# DESeq2 差异表达分析模板
# 输入：run_rnaseq.sh 产出的 counts_matrix.tsv + samples.tsv
# 输出：标准化表达矩阵、差异结果表、火山图、PCA 图、样本距离热图
###############################################################################

suppressPackageStartupMessages({
  library(DESeq2)
  library(ggplot2)
  library(pheatmap)
})

## ------------------------------ 参数 ---------------------------------------
count_file   <- "results/04_counts/counts_matrix.tsv"   # 计数矩阵
coldata_file <- "samples.tsv"                            # 样本表（sample_id / group）
out_dir      <- "results/06_DE"
ref_level    <- "control"    # 对照组（参考水平），按你的分组改
padj_cut     <- 0.05         # 显著性阈值
lfc_cut      <- 1            # log2FC 阈值
dir.create(out_dir, showWarnings = FALSE, recursive = TRUE)

## ---------------------------- 读入数据 -------------------------------------
counts <- read.delim(count_file, row.names = 1, check.names = FALSE)
counts <- as.matrix(counts)

coldata <- read.delim(coldata_file, stringsAsFactors = FALSE,
                      comment.char = "#")
rownames(coldata) <- coldata$sample_id
coldata$group <- factor(coldata$group)
coldata$group <- relevel(coldata$group, ref = ref_level)

## 对齐样本顺序（这一步漏了会直接报错或结果错乱）
common <- intersect(colnames(counts), rownames(coldata))
if (length(common) == 0) stop("计数矩阵列名与样本表 sample_id 对不上，请检查")
counts  <- counts[, common, drop = FALSE]
coldata <- coldata[common, , drop = FALSE]

message("样本数：", ncol(counts), "；基因数：", nrow(counts))
message("分组：", paste(levels(coldata$group), collapse = " vs "))

## --------------------------- 构建 DESeq 对象 -------------------------------
dds <- DESeqDataSetFromMatrix(countData = counts,
                              colData   = coldata,
                              design    = ~ group)

## 过滤低表达基因：至少在一半样本里 count >= 10
keep <- rowSums(counts(dds) >= 10) >= floor(ncol(counts) / 2)
dds  <- dds[keep, ]
message("过滤后保留基因数：", nrow(dds))

## ---------------------- 标准化 + 差异分析 ----------------------------------
dds <- DESeq(dds)

## 提取结果（本例取 实验组 vs 对照组）
res <- results(dds, contrast = c("group",
                                 setdiff(levels(coldata$group), ref_level)[1],
                                 ref_level))
res <- res[order(res$padj), ]

## 收缩 log2FC（推荐，减少低表达基因的假阳性）
res_shrunk <- lfcShrink(dds, coef = resultsNames(dds)[2], type = "apeglm")

## ------------------------------ 导出结果 -----------------------------------
norm_counts <- counts(dds, normalized = TRUE)
write.table(norm_counts, file = file.path(out_dir, "normalized_counts.tsv"),
            sep = "\t", quote = FALSE, col.names = NA)

res_df <- as.data.frame(res_shrunk)
res_df$gene <- rownames(res_df)
write.table(res_df, file = file.path(out_dir, "DE_results_all.tsv"),
            sep = "\t", quote = FALSE, row.names = FALSE)

## 显著差异基因
sig <- subset(res_df, !is.na(padj) & padj < padj_cut & abs(log2FoldChange) > lfc_cut)
write.table(sig, file = file.path(out_dir, "DE_genes_significant.tsv"),
            sep = "\t", quote = FALSE, row.names = FALSE)
message("显著差异基因：", nrow(sig), " 个（上调 ",
        sum(sig$log2FoldChange > 0), "，下调 ", sum(sig$log2FoldChange < 0), "）")

## ------------------------------ 可视化 -------------------------------------

## 1) 火山图
volcano_df <- as.data.frame(res_shrunk)
volcano_df$significance <- "NotSig"
volcano_df$significance[volcano_df$padj < padj_cut &
                        volcano_df$log2FoldChange >  lfc_cut] <- "Up"
volcano_df$significance[volcano_df$padj < padj_cut &
                        volcano_df$log2FoldChange < -lfc_cut] <- "Down"

p_volcano <- ggplot(volcano_df,
                    aes(x = log2FoldChange, y = -log10(padj),
                        color = significance)) +
  geom_point(alpha = 0.7, size = 1.2) +
  scale_color_manual(values = c(Up = "#c0392b", Down = "#2471a3",
                                NotSig = "grey75")) +
  geom_vline(xintercept = c(-lfc_cut, lfc_cut), linetype = "dashed") +
  geom_hline(yintercept = -log10(padj_cut), linetype = "dashed") +
  labs(x = "log2 fold change", y = "-log10(padj)", color = NULL,
       title = "Differential expression") +
  theme_bw(base_size = 13)

ggsave(file.path(out_dir, "volcano.pdf"), p_volcano, width = 6, height = 5)
ggsave(file.path(out_dir, "volcano.png"), p_volcano, width = 6, height = 5,
       dpi = 300)

## 2) PCA 图（用方差稳定变换后的数据）
vsd <- vst(dds, blind = FALSE)
pca <- plotPCA(vsd, intgroup = "group", returnData = TRUE)
pct <- round(100 * attr(pca, "percentVar"))

p_pca <- ggplot(pca, aes(PC1, PC2, color = group, label = name)) +
  geom_point(size = 3) +
  geom_text(vjust = -1, size = 3, show.legend = FALSE) +
  xlab(paste0("PC1: ", pct[1], "% variance")) +
  ylab(paste0("PC2: ", pct[2], "% variance")) +
  theme_bw(base_size = 13)

ggsave(file.path(out_dir, "PCA.pdf"), p_pca, width = 6, height = 5)

## 3) 样本相关性热图
sample_dist <- dist(t(assay(vsd)))
dist_mat <- as.matrix(sample_dist)
pdf(file.path(out_dir, "sample_distance_heatmap.pdf"), width = 7, height = 6)
pheatmap(dist_mat, clustering_distance_rows = sample_dist,
         clustering_distance_cols = sample_dist,
         annotation_col = data.frame(group = coldata$group,
                                     row.names = rownames(coldata)),
         main = "Sample-to-sample distance")
dev.off()

## 4) 显著基因热图（取前 50 个最显著基因）
if (nrow(sig) > 1) {
  top_n <- head(sig$gene, 50)
  pdf(file.path(out_dir, "top50_heatmap.pdf"), width = 8, height = 10)
  pheatmap(assay(vsd)[top_n, ], scale = "row",
           annotation_col = data.frame(group = coldata$group,
                                       row.names = rownames(coldata)),
           show_rownames = TRUE, fontsize_row = 6,
           main = "Top 50 DE genes")
  dev.off()
}

message("全部完成，结果见：", normalizePath(out_dir))

###############################################################################
# 依赖安装：
#   BiocManager::install(c("DESeq2"))
#   install.packages(c("ggplot2","pheatmap"))
# 运行：
#   Rscript deseq2_de.R
# 注意：
#   1) apeglm 收缩需要 install.packages("apeglm")；没装可把 type 改成 "ashr"
#      或直接不收缩，用 res（注意此时 log2FC 更保守）
#   2) 分组水平名必须是合法 R 变量名，避免 "ctrl-x" 这类带横线的名字
###############################################################################

#!/usr/bin/env Rscript
###############################################################################
# 功能富集与可视化模板（clusterProfiler）
# 输入：deseq2_de.R 产出的 DE_genes_significant.tsv
# 输出：GO / KEGG 富集结果表 + 气泡图 + 条形图，附带 GSEA 示例
###############################################################################

suppressPackageStartupMessages({
  library(clusterProfiler)
  library(org.Hs.eg.db)      # 人类注释，换物种需换对应的 org.*.eg.db
  library(ggplot2)
  library(enrichplot)
})

## ------------------------------ 参数 ---------------------------------------
deg_file  <- "results/06_DE/DE_genes_significant.tsv"
all_file  <- "results/06_DE/DE_results_all.tsv"
out_dir   <- "results/07_enrichment"
padj_cut  <- 0.05
dir.create(out_dir, showWarnings = FALSE, recursive = TRUE)

## ---------------------------- 读入基因列表 ---------------------------------
deg <- read.delim(deg_file, stringsAsFactors = FALSE)
up_genes   <- deg$gene[deg$log2FoldChange > 0]
down_genes <- deg$gene[deg$log2FoldChange < 0]
message("上调基因 ", length(up_genes), "，下调基因 ", length(down_genes))

## 基因 ID 转换（Ensembl -> Entrez）
convert <- function(ids) {
  mapIds(org.Hs.eg.db, keys = ids, column = "ENTREZID",
         keytype = "ENSEMBL", multiVals = "first")
}

## ------------------------- 1. GO 富集（过表征分析） ------------------------
run_go <- function(ids, tag) {
  eg <- bitr(ids, fromType = "ENSEMBL", toType = "ENTREZID",
             OrgDb = org.Hs.eg.db)
  if (nrow(eg) == 0) return(NULL)

  for (ont in c("BP", "CC", "MF")) {
    res <- enrichGO(gene = eg$ENTREZID, OrgDb = org.Hs.eg.db,
                    ont = ont, keyType = "ENTREZID",
                    pAdjustMethod = "BH", pvalueCutoff = 0.05,
                    qvalueCutoff = 0.2, readable = TRUE)
    if (is.null(res) || nrow(as.data.frame(res)) == 0) next

    write.table(as.data.frame(res),
                file = file.path(out_dir, paste0("GO_", ont, "_", tag, ".tsv")),
                sep = "\t", quote = FALSE, row.names = FALSE)

    p <- dotplot(res, showCategory = 15) +
      ggtitle(paste0("GO ", ont, " enrichment (", tag, ")"))
    ggsave(file.path(out_dir, paste0("GO_", ont, "_", tag, ".pdf")),
           p, width = 8, height = 7)
  }
}

run_go(up_genes,   "up")
run_go(down_genes, "down")

## ---------------------------- 2. KEGG 富集 ---------------------------------
run_kegg <- function(ids, tag) {
  eg <- bitr(ids, fromType = "ENSEMBL", toType = "ENTREZID",
             OrgDb = org.Hs.eg.db)
  if (nrow(eg) == 0) return(NULL)

  kk <- enrichKEGG(gene = eg$ENTREZID, organism = "hsa",
                   keyType = "kegg", pvalueCutoff = 0.05,
                   qvalueCutoff = 0.2)
  if (is.null(kk) || nrow(as.data.frame(kk)) == 0) return(NULL)
  kk <- setReadable(kk, OrgDb = org.Hs.eg.db, keyType = "ENTREZID")

  write.table(as.data.frame(kk),
              file = file.path(out_dir, paste0("KEGG_", tag, ".tsv")),
              sep = "\t", quote = FALSE, row.names = FALSE)

  p <- barplot(kk, showCategory = 15, x = "GeneRatio") +
    ggtitle(paste0("KEGG enrichment (", tag, ")"))
  ggsave(file.path(out_dir, paste0("KEGG_", tag, ".pdf")),
         p, width = 8, height = 7)
}

run_kegg(up_genes,   "up")
run_kegg(down_genes, "down")

## --------------------------- 3. GSEA（进阶） -------------------------------
## GSEA 用全部基因按 log2FC 排序，不设阈值，优势是不丢信息
all_res <- read.delim(all_file, stringsAsFactors = FALSE)
all_res <- all_res[!is.na(all_res$padj) & !is.na(all_res$entrez),
                   c("gene", "log2FoldChange")]
gene_list <- all_res$log2FoldChange
names(gene_list) <- all_res$gene
gene_list <- sort(gene_list, decreasing = TRUE)

## 把 Ensembl 换成 Entrez
eg_all <- bitr(names(gene_list), fromType = "ENSEMBL", toType = "ENTREZID",
               OrgDb = org.Hs.eg.db)
gene_list <- gene_list[eg_all$ENSEMBL]
names(gene_list) <- eg_all$ENTREZID

gse <- gseGO(geneList = gene_list, OrgDb = org.Hs.eg.db, ont = "BP",
             keyType = "ENTREZID", pvalueCutoff = 0.05,
             minGSSize = 10, maxGSSize = 500)
if (!is.null(gse) && nrow(as.data.frame(gse)) > 0) {
  write.table(as.data.frame(gse),
              file = file.path(out_dir, "GSEA_GO_BP.tsv"),
              sep = "\t", quote = FALSE, row.names = FALSE)
  pdf(file.path(out_dir, "GSEA_ridgeplot.pdf"), width = 8, height = 8)
  print(ridgeplot(gse, showCategory = 15))
  dev.off()
}

message("富集分析完成，结果见：", normalizePath(out_dir))

###############################################################################
# 依赖安装：
#   BiocManager::install(c("clusterProfiler","org.Hs.eg.db","enrichplot"))
# 换物种：小鼠用 org.Mm.eg.db + organism = "mmu"；其余同理
# 运行：Rscript enrichment_plot.R
# 常见坑：
#   1) 基因 ID 类型必须说清（ENSEMBL / SYMBOL / ENTREZID），搞错会返回空结果
#   2) 网络不通时 enrichKEGG 会失败，可改 enrichGO（纯本地注释，不联网）
#   3) 富集结果没有条目数，通常是基因太少或 ID 没对上，先检查 bitr 的返回
###############################################################################

# RNA-seq 分析环境依赖清单
# 用法：mamba env create -f environment.yaml        （或 conda env create -f ...）
#       conda activate rnaseq
# 说明：下面写的是撰写时的版本号，仅作参考。若创建环境时报找不到某个版本
#       （solvability / PackagesNotFoundError），把对应包的 =版本号 删掉即可，
#       conda 会自动选当前可用版本。正式分析前请保存 conda list 的输出以便复现。
name: rnaseq
channels:
  - bioconda
  - conda-forge
  - defaults
dependencies:
  # ---- 基础工具 ----
  - python=3.11
  - pip
  # ---- 质控与修剪 ----
  - fastqc=0.12.1
  - fastp=0.23.4
  - trim-galore=0.6.10
  - multiqc=1.21
  - sortmerna=4.3.6
  # ---- 比对 ----
  - bwa=0.7.18
  - bowtie2=2.5.3
  - hisat2=2.2.1
  - star=2.7.11b
  - minimap2=2.28
  # ---- 处理 SAM/BAM/VCF ----
  - samtools=1.20
  - bcftools=1.20
  - bedtools=2.31.1
  - picard=3.1.1
  # ---- 序列工具 ----
  - seqkit=2.8.1
  - mafft=7.526
  # ---- 定量 ----
  - subread=2.0.6          # 提供 featureCounts
  - salmon=1.10.2          # 伪比对定量
  - kallisto=0.50.1
  - rsem=1.3.3
  - stringtie=2.2.3
  # ---- 表观 / 变异 ----
  - macs3=3.0.1
  - deeptools=3.5.4
  - snpeff=5.1
  # ---- 流程引擎（可选，跑 nf-core 时用）----
  - nextflow=24.10.0
  # ---- R 与常用包 ----
  - r-base=4.3
  - r-tidyverse
  - r-ggplot2
  - r-pheatmap
  - r-rcolorbrewer
  - bioconductor-deseq2
  - bioconductor-edger
  - bioconductor-tximport
  - bioconductor-clusterprofiler
  - bioconductor-org.hs.eg.db
  - bioconductor-enrichplot
  - bioconductor-complexheatmap
  - bioconductor-chipseeker