Nextflow STAR+DESeq2流程故障：ReadsPerGene.out.tab输入异常

问题背景 #

构建了一个Nextflow流程，用于处理对照组和处理组的成对fastq文件，流程步骤包括FASTQC质控、STAR基因组索引与比对、DESeq2差异分析。但DESeq2模块执行时报错：Error: object 'Control' not found。

错误原因 #

问题出在DESEQ进程的R脚本中单引号嵌套冲突：
外层使用Rscript -e '...'包裹R代码，内部定义字符串时又使用了单引号（如c('Control', 'Treatment')），导致R解析时提前截断字符串，将Control识别为未定义的对象而非字符串常量。

修复方案 #

将R脚本内部的单引号替换为双引号，避免与外层的单引号冲突。同时补充两处流程健壮性优化：

1. 读取STAR输出的ReadsPerGene.out.tab时，跳过前4行（STAR输出的前3行是统计信息，第4行才是有效基因计数）
2. 为coldata设置行名，确保与countData的列名严格匹配，符合DESeq2的输入要求

修改后的完整代码 #

#!/usr/bin/env nextflow

nextflow.enable.dsl=2

params.reads_control_1 = "./data/control/subset_control_1.fastq"
params.reads_control_2 = "./data/control/subset_control_2.fastq"
params.reads_treatment_1 = "./data/treatment/subset_treatment_1.fastq"
params.reads_treatment_2 = "./data/treatment/subset_treatment_2.fastq"
params.genome = './data/ggal/genome.fna'
params.gtf = './data/ggal/annotations.gtf'
params.star_index = './starindex'
params.output = './results'


process FASTQC {
    container 'quay.io/biocontainers/fastqc:0.12.1--hdfd78af_0'
    input:
    path reads_control_1 
    path reads_control_2 
    path reads_treatment_1 
    path reads_treatment_2

    script:
    """
    fastqc -o . $reads_control_1 $reads_control_2 $reads_treatment_1 $reads_treatment_2
    """
}


process STARINDEX {
    container 'quay.io/biocontainers/star:2.7.11b--h43eeafb_3'
    input:
    path genome
    path gtf

    output:
    path params.star_index
        
    script:
    """
    STAR --runThreadN 8 --runMode genomeGenerate \
         --genomeDir ${params.star_index} \
         --genomeFastaFiles $genome \
         --sjdbGTFfile $gtf
    """
}
        
        
process STAR_ALIGN_CONTROL {
    container 'quay.io/biocontainers/star:2.7.11b--h43eeafb_3'
    input:
    path reads_control_1
    path reads_control_2
    path star_index

    output:
    path 'controlReadsPerGene.out.tab' 

    script:
    """
    STAR --runThreadN 8 --genomeDir $star_index \
         --readFilesIn $reads_control_1 $reads_control_2 \
         --outSAMtype BAM SortedByCoordinate \
         --quantMode GeneCounts \
         --outFileNamePrefix control 
    """
}

process STAR_ALIGN_TREATMENT {
    container 'quay.io/biocontainers/star:2.7.11b--h43eeafb_3'
    input:
    path reads_treatment_1
    path reads_treatment_2
    path star_index

    output:
    path 'treatmentReadsPerGene.out.tab' 

    script:
    """
    STAR --runThreadN 8 --genomeDir $star_index \
         --readFilesIn $reads_treatment_1 $reads_treatment_2 \
         --outSAMtype BAM SortedByCoordinate \
         --quantMode GeneCounts \
         --outFileNamePrefix treatment 
    """
}


process DESEQ {
    container 'quay.io/biocontainers/bioconductor-deseq2:1.42.0--r43hf17093f_2'
    input:
    path control_reads_ch
    path treatment_reads_ch

    output:
    path "deseq2_results.csv"

    script:
    """
    Rscript -e '
    library(DESeq2)
    # 跳过STAR输出的前4行，仅读取基因计数数据
    control_reads <- read.table("${control_reads_ch}", header=FALSE, skip=4)
    treatment_reads <- read.table("${treatment_reads_ch}", header=FALSE, skip=4)
    combined_counts <- cbind(control_reads[,4], treatment_reads[,4])
    colnames(combined_counts) <- c("Control", "Treatment")
    rownames(combined_counts) <- control_reads[,1]
    # 确保coldata行名与countData列名完全匹配
    coldata <- data.frame(condition=factor(c("Control", "Treatment")), row.names=colnames(combined_counts))
    dds <- DESeqDataSetFromMatrix(countData=combined_counts, colData=coldata, design=~condition)
    dds <- DESeq(dds)
    results <- results(dds)
    write.csv(as.data.frame(results), "deseq2_results.csv")
    '
    """
}


workflow {
    control1_ch = Channel.fromPath(params.reads_control_1)
    control2_ch = Channel.fromPath(params.reads_control_2)
    treatment1_ch = Channel.fromPath(params.reads_treatment_1)
    treatment2_ch = Channel.fromPath(params.reads_treatment_2)

    FASTQC(control1_ch, control2_ch, treatment1_ch, treatment2_ch)

    genome_ch = Channel.fromPath(params.genome)
    gtf_ch = Channel.fromPath(params.gtf)
    
    star_index = STARINDEX(genome_ch, gtf_ch)

    control_reads_ch = STAR_ALIGN_CONTROL(control1_ch, control2_ch, star_index)
    treatment_reads_ch = STAR_ALIGN_TREATMENT(treatment1_ch, treatment2_ch, star_index)
    
    DESEQ(control_reads_ch, treatment_reads_ch)
}

内容的提问来源于stack exchange，提问作者wawawawa