nf-core 流程部署：应用场景

nf-core 流程部署：应用场景

Written By
技能练习生
技能练习生

本指南通过真实场景展示 nf-core 流程的典型应用，帮助你理解如何在不同研究项目中使用这些工具。

案例 1：RNA-seq 差异表达分析

研究背景

某生物学实验室研究药物处理对癌细胞系的影响，需要进行转录组分析以识别差异表达基因。

研究设计

  • 样本类型：人类癌细胞系
  • 实验条件：药物处理组 vs. 对照组
  • 生物学重复：每组 3 个重复（共 6 个样本）
  • 测序平台：Illumina NovaSeq 2×150 bp

执行流程

步骤 1：环境准备

# 运行环境检查
python scripts/check_environment.py

# 检查所需基因组
python scripts/manage_genomes.py check GRCh38

步骤 2：样本表准备

使用自动化脚本生成样本表：

# 从 FASTQ 文件生成样本表
python scripts/generate_samplesheet.py \
    /path/to/fastq_data \
    rnaseq \
    -o samplesheet.csv

# 验证样本表
python scripts/generate_samplesheet.py --validate samplesheet.csv rnaseq

生成的样本表格式：

sample,fastq_1,fastq_2,strandedness
control_rep1,/data/fastq/control_rep1_R1.fastq.gz,/data/fastq/control_rep1_R2.fastq.gz,auto
control_rep2,/data/fastq/control_rep2_R1.fastq.gz,/data/fastq/control_rep2_R2.fastq.gz,auto
control_rep3,/data/fastq/control_rep3_R1.fastq.gz,/data/fastq/control_rep3_R2.fastq.gz,auto
treated_rep1,/data/fastq/treated_rep1_R1.fastq.gz,/data/fastq/treated_rep1_R2.fastq.gz,auto
treated_rep2,/data/fastq/treated_rep2_R1.fastq.gz,/data/fastq/treated_rep2_R2.fastq.gz,auto
treated_rep3,/data/fastq/treated_rep3_R1.fastq.gz,/data/fastq/treated_rep3_R2.fastq.gz,auto

步骤 3：运行流程

# 运行 nf-core/rnaseq 流程
nextflow run nf-core/rnaseq \
    -r 3.22.2 \
    -profile docker \
    --input samplesheet.csv \
    --outdir results \
    --genome GRCh38 \
    --aligner star_salmon \
    --skip_pseudo_aligner false \
    -resume

步骤 4：结果解读

关键输出文件：

  1. 基因计数矩阵：results/star_salmon/salmon.merged.gene_counts.tsv

    • 用于下游差异表达分析
    • 可直接导入 DESeq2 或 edgeR

  2. TPM 表达值：results/star_salmon/salmon.merged.gene_tpm.tsv

    • 标准化的表达值
    • 便于样本间比较

  3. 质控报告：results/multiqc/multiqc_report.html

    • 整体质控评估
    • 识别潜在问题样本

  4. 比对统计：results/star_salmon/logs/alignment_stats/

    • 比对率、mapping 质量等

下游分析建议

完成流程后，使用差异表达工具分析结果：

# 使用 DESeq2 进行差异表达分析
library(DESeq2)

# 读取计数矩阵
count_data <- read.table("results/star_salmon/salmon.merged.gene_counts.tsv",
                          header=TRUE, row.names=1)

# 创建样本信息表
sample_info <- data.frame(
    condition = c(rep("control", 3), rep("treated", 3))
)

# 运行 DESeq2
dds <- DESeqDataSetFromMatrix(countData = count_data,
                              colData = sample_info,
                              design = ~ condition)
dds <- DESeq(dds)
res <- results(dds)

案例 2：全外显子组变异检测

研究背景

某遗传学研究团队分析家族性遗传病的致病突变，对家系成员进行全外显子组测序。

研究设计

  • 样本类型：人类外显子组测序数据
  • 家系结构：先证者 + 父母（trio 设计）
  • 测序深度：平均 100× 覆盖
  • 目标：识别隐性遗传致病突变

执行流程

步骤 1：环境检查

# 验证环境
python scripts/check_environment.py

# 确认基因组索引
python scripts/manage_genomes.py check GRCh38

步骤 2：样本表准备

Sarek 流程需要特殊的样本表格式，包含肿瘤/正常状态信息：

# 生成 sarek 样本表
python scripts/generate_samplesheet.py \
    /path/to/wes_data \
    sarek \
    -o samplesheet.csv

手动编辑样本表，标记样本状态：

patient,sample,lane,fastq_1,fastq_2,status
family1_proband,proband,L001,/data/wes/proband_R1.fastq.gz,/data/wes/proband_R2.fastq.gz,0
family1_father,father,L001,/data/wes/father_R1.fastq.gz,/data/wes/father_R2.fastq.gz,0
family1_mother,mother,L001,/data/wes/mother_R1.fastq.gz,/data/wes/mother_R2.fastq.gz,0

步骤 3：运行变异检测流程

# 使用 GATK HaplotypeCaller 进行生殖系变异检测
nextflow run nf-core/sarek \
    -r 3.7.1 \
    -profile docker \
    --input samplesheet.csv \
    --outdir results \
    --genome GRCh38 \
    --tools haplotypecaller,variantannotator \
    --skip_bqsr false \
    -resume

关键参数说明：

  • --tools haplotypecaller：使用 GATK 最佳实践进行变异检测
  • --skip_bqsr false：执行碱基质量分数重校正（提高准确性）
  • --variantannotator：添加功能注释

步骤 4：结果解读

关键输出：

  1. VCF 文件：results/variant_calling/haplotypecaller/

    • 包含所有检测到的变异
    • 按样本分割

  2. 联合基因分型：results/variant_calling/genotypegvcfs/

    • 家系样本的联合 VCF

  3. 质控指标：results/multiqc/multiqc_report.html

    • 覆盖度统计
    • 变异质量分布

变异筛选策略

根据家系设计进行变异筛选：

# 使用 bcftools 筛选隐性遗传突变
bcftools view -i 'GT="0/1" && INFO/AF<0.01' joint_genotyping.vcf.gz | \
bcftools filter -i 'INFO/DP>20' > candidate_variants.vcf

筛选标准：

  • 新发突变或符合隐性遗传模式
  • 人群频率 < 1%（gnomAD）
  • 在致病基因列表中
  • 预测功能影响（如错义、无义突变）

案例 3：ATAC-seq 染色质可及性分析

研究背景

研究干细胞分化过程中染色质开放性的动态变化，识别关键调控区域。

研究设计

  • 细胞类型：胚胎干细胞和分化细胞
  • 时间点：分化第 0、2、5、10 天
  • 生物学重复：每个时间点 2 个重复
  • 测序平台：Illumina HiSeq 2×50 bp

执行流程

步骤 1：环境准备

# 环境检查
python scripts/check_environment.py

# 检查基因组
python scripts/manage_genomes.py check GRCh38

步骤 2：样本表准备

# 生成 ATAC-seq 样本表
python scripts/generate_samplesheet.py \
    /path/to/atac_data \
    atacseq \
    -o samplesheet.csv

样本表示例：

sample,fastq_1,fastq_2,replicate
day0_rep1,/data/atac/day0_rep1_R1.fastq.gz,/data/atac/day0_rep1_R2.fastq.gz,1
day0_rep2,/data/atac/day0_rep2_R1.fastq.gz,/data/atac/day0_rep2_R2.fastq.gz,1
day2_rep1,/data/atac/day2_rep1_R1.fastq.gz,/data/atac/day2_rep1_R2.fastq.gz,1
day2_rep2,/data/atac/day2_rep2_R1.fastq.gz,/data/atac/day2_rep2_R2.fastq.gz,1
day5_rep1,/data/atac/day5_rep1_R1.fastq.gz,/data/atac/day5_rep1_R2.fastq.gz,1
day5_rep2,/data/atac/day5_rep2_R1.fastq.gz,/data/atac/day5_rep2_R2.fastq.gz,1
day10_rep1,/data/atac/day10_rep1_R1.fastq.gz,/data/atac/day10_rep2_R2.fastq.gz,1
day10_rep2,/data/atac/day10_rep2_R1.fastq.gz,/data/atac/day10_rep2_R2.fastq.gz,1

步骤 3：运行 ATAC-seq 流程

# 运行 nf-core/atacseq 流程
nextflow run nf-core/atacseq \
    -r 2.1.2 \
    -profile docker \
    --input samplesheet.csv \
    --outdir results \
    --genome GRCh38 \
    --read_length 50 \
    --skip_consensus_peaks false \
    -resume

关键参数：

  • --read_length 50：根据测序读长设置
  • --skip_consensus_peaks false：生成一致的峰集合

步骤 4：结果解读

关键输出：

  1. 峰文件：results/macs2/narrowPeak/

    • 每个样本的开放染色质区域
    • BED 格式，可用于下游分析

  2. 一致峰：results/consensus_peak_set/

    • 所有样本共有的峰
    • 用于定量比较

  3. 可视化文件：results/bwa/mergedLibrary/bigwig/

    • 覆盖度轨道
    • 可在基因组浏览器中查看

  4. 质控指标：

    • TSS 富集分数（评估文库质量）
    • FRiP（fraction of reads in peaks）
    • 序列片段长度分布

下游分析

使用 R/Bioconductor 进行差异可及性分析：

library(diffbind)

# 读取峰文件
samples <- read.table("samplesheet.csv", header=TRUE)
peaks <- read.table("results/consensus_peak_set/consensus_peaks.bed")

# 创建 DBA 对象
dba_obj <- dba(sampleSheet="samplesheet.csv")

# 计算差异可及性
dba_obj <- dba.count(dba_obj, peaks=peaks)
dba_obj <- dba.contrast(dba_obj, categories=DBA_CONDITION)
dba_obj <- dba.analyze(dba_obj)

# 提取显著差异的峰
diff_peaks <- dba.report(dba_obj, th=0.05)

案例 4：公共数据重分析（GEO/SRA）

研究场景

验证已发表研究的结果，或进行跨研究整合分析。

执行流程

步骤 1：获取数据

以 GEO 数据集 GSE110004 为例：

# 获取研究信息
python scripts/sra_geo_fetch.py info GSE110004

# 下载数据（交互模式）
python scripts/sra_geo_fetch.py download GSE110004 -o ./fastq -i

# 生成样本表
python scripts/sra_geo_fetch.py samplesheet GSE110004 \
    --fastq-dir ./fastq \
    -o samplesheet.csv

步骤 2：选择并运行流程

根据数据类型选择合适的流程：

# 假设是 RNA-seq 数据
nextflow run nf-core/rnaseq \
    -r 3.22.2 \
    -profile docker \
    --input samplesheet.csv \
    --outdir results \
    --genome GRCh38 \
    -resume

重分析注意事项

  1. 基因组版本一致性

    • 原研究使用的基因组版本
    • 基因注释版本对结果影响很大

  2. 流程参数对比

    • 尽可能使用与原研究相似的参数
    • 记录所有偏离参数

  3. 结果验证

    • 与原始研究的关键结果对比
    • 重点关注排序和显著性

常见问题和解决方案

问题 1：样本表格式错误

症状：流程启动后立即失败，提示 "Invalid samplesheet format"

解决方案：

# 使用验证脚本检查样本表
python scripts/generate_samplesheet.py --validate samplesheet.csv rnaseq

# 常见问题：
# - 路径不是绝对路径
# - 缺少必需列
# - 混合使用不同的行尾符（Windows/Mac/Linux）

问题 2：内存不足

症状：某些步骤被 killed，日志显示 "Out of memory"

解决方案：

# 限制资源使用
nextflow run nf-core/rnaseq \
    --input samplesheet.csv \
    --genome GRCh38 \
    --max_memory '128.GB' \
    --max_cpus 16 \
    -resume

问题 3：基因组索引不存在

症状：错误 "Genome GRCh38 not found in iGenomes cache"

解决方案：

# 下载基因组索引
python scripts/manage_genomes.py download GRCh38

# 或使用自定义基因组
nextflow run nf-core/rnaseq \
    --igenomes_base /path/to/custom/genomes \
    --genome GRCh38 \
    --input samplesheet.csv

最佳实践总结

  1. 从小规模测试开始

    • 先用测试数据验证环境
    • 确保流程完整运行后再投入真实数据

  2. 保留中间文件

    • 使用 -resume 参数可从断点恢复
    • 节省计算时间和资源

  3. 质控先行

    • 先检查 MultiQC 报告
    • 发现问题及早终止，避免浪费资源

  4. 详细记录

    • 保存所有使用的命令和参数
    • 记录环境和软件版本
    • 便于结果复现和问题排查

  5. 合理规划存储

    • 原始数据 + 中间文件 + 结果可能需要数 TB
    • 及时清理不必要的中间文件
    • 归档重要的结果文件
下一部分nf-core 流程部署：工作原理