Snakemake中用cutadapt实现Nanopore测序demultiplex的输出异常问题

你当前遇到的12个输出文件内容完全一致、规则重复运行的问题，本质是两处设计偏差：

• 原demultiplex规则仅将单个{name}.tmp文件声明为输出，Snakemake会按照通配符匹配逻辑，为12个样本名分别触发一次独立的规则运行
• 每次运行cutadapt时，仅指定了单个输出路径，没有使用cutadapt针对多重条码的拆分输出逻辑，命令会把所有匹配到任意条码的read全部写入当前对应的输出文件，不会按条码归属拆分内容，因此12次运行生成的文件内容完全相同

要实现“单组输入触发1次规则运行、一次性生成所有对应条码的独立输出文件”的目标，需要调整规则的输出声明逻辑和cutadapt的命令参数：

1. 调整规则的输出声明 #

因为规则仅需要运行1次、生成全部12个样本的输出文件，需要把所有待生成的样本文件全部列入规则的output列表，Snakemake识别到所有目标输出属于同一次运行的产物后，就不会重复触发规则。

2. 启用cutadapt的多文件输出能力 #

cutadapt支持直接读取fasta格式的条码序列文件，通过内置的占位符，自动将匹配到不同条码的read拆分到对应命名的独立文件中，无需为每个样本单独运行一次命令。
需要注意：Snakemake的shell块中，Snakemake自身的通配符是用单花括号{}标记的，传递给cutadapt的自定义占位符需要用双花括号{{}}标记，避免被Snakemake提前替换。

调整后的完整规则参考代码如下：

rule demultiplex:
    """ Demultiplex Nextera-styled barcodes using the compiled barcode fasta file """
    input:
        barcodeFile=os.path.join(outDir, "barcodes_used.fasta"),
        fastqFile=os.path.join(outDir, "reads_trim.fq")
    output:
        # 列出该规则需要生成的全部样本输出，Snakemake仅会触发1次运行
        expand(os.path.join(outDir, "{name}.tmp"), name=samples_list)
    log: "logs/cutadapt/demultiplexing.log"
    shell:
        """
        cutadapt -e 0.20 -O 10 -m 15 -M 15 \
          -g file:{input.barcodeFile} \
          -o {outDir}/{{adapter_name}}.tmp \
          {input.fastqFile}
        """

实现逻辑说明 #

• Snakemake渲染shell命令时，会把单花括号的内容替换为实际的文件路径、目录值，双花括号的{{adapter_name}}会被渲染为cutadapt可识别的{adapter_name}占位符
• cutadapt运行时会逐行读取barcodes_used.fasta中的条码序列，将序列标题作为adapter_name的值，把匹配到对应条码的read写入{adapter_name}.tmp文件，实现按样本拆分的需求
• 因为规则的output列表已经包含了所有样本的输出文件，Snakemake构建有向无环图时，会把所有目标文件归为同一次规则运行的产物，不会重复执行12次

• 请确保samples_list中的样本名，和barcodes_used.fasta文件里的序列标题完全匹配，避免cutadapt生成的文件名和Snakemake声明的输出文件名不一致，触发缺失输出的错误
• 如果需要排查未匹配read的情况，可以在cutadapt命令中增加--unmatched {outDir}/unmatched.tmp参数，将没有匹配到任何条码的read单独输出到文件
• 如果使用的cutadapt版本不支持{adapter_name}占位符，可参考对应版本的说明文档调整占位符写法，核心逻辑是保持cutadapt一次运行完成所有条码的拆分，不要为每个样本单独触发命令
• 不需要给单次运行的规则配置带通配符的日志路径，全局日志即可满足排查需求

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