使用 Snakemake 的 fastqc

Question

我有一个通过 Snakemake 的样本列表。当我到达我的 fastqc 步骤时，我突然发现每个样本有两个文件（R1 和 R2 文件）。考虑以下规则：

rule fastqc:
    input:
        os.path.join(fastq_dir, '{sample}_R1_001.fastq.gz'),
        os.path.join(fastq_dir, '{sample}_R2_001.fastq.gz')
    output:
        os.path.join(fastq_dir, '{sample}_R1_fastq.html'),
        os.path.join(fastq_dir, '{sample}_R2_fastq.html')
    conda:
        "../envs/fastqc.yaml"
    shell:
        '''
        #!/bin/bash
        fastqc {input} --outdir={fastqc_dir}
        '''

这不起作用。我还尝试了以下方法：

rule fastqc:
    input:
        expand([os.path.join(fastq_dir, '{sample}_R{read}_001.fastq.gz')], read=['1', '2']
    output:
        expand([os.path.join(fastq_dir, '{sample}_R{read}_fastq.html')], read=['1', '2']
    conda:
        "../envs/fastqc.yaml"
    shell:
        '''
        #!/bin/bash
        fastqc {input} --outdir={fastqc_dir}
        '''

这也不起作用，我得到：

No values given for wildcard 'sample'.

然后我尝试了：

rule fastqc:
    input:
        expand([os.path.join(fastq_dir, '{sample}_R{read}_001.fastq.gz')], read=['1', '2'], sample=samples['samples'])
    output:
        expand([os.path.join(fastqc_dir, '{sample}_R{read}_fastqc.html')], read=['1', '2'], sample=samples['samples'])
    conda:
        "../envs/fastqc.yaml"
    shell:
        '''
        #!/bin/bash
        fastqc {input} --outdir={fastqc_dir}
        '''

但这似乎将所有 fastq 文件都输入到一个 shell 脚本中。

我应该如何正确地 "loop" 处理 1 个样本的多个输入？

最崇高的问候。

编辑：

我的规则都是这样的，也许我也应该改变它，对吧（请参阅 fastqc 的最后两行）？

# Rule all is a pseudo-rule that tells snakemake what final files to generate.
rule all:
    input:
        expand([os.path.join(analyzed_dir, '{sample}.genes.results'),
                os.path.join(rseqc_dir, '{sample}.bam_stat.txt'),
                os.path.join(rseqc_dir, '{sample}.clipping_profile.xls'),
                os.path.join(rseqc_dir, '{sample}.deletion_profile.txt'),
                os.path.join(rseqc_dir, '{sample}.infer_experiment.txt'),
                os.path.join(rseqc_dir, '{sample}.geneBodyCoverage.txt'),
                os.path.join(rseqc_dir, '{sample}.inner_distance.txt'),
                os.path.join(rseqc_dir, '{sample}.insertion_profile.xls'),
                os.path.join(rseqc_dir, '{sample}.junction.xls'),
                os.path.join(rseqc_dir, '{sample}.junctionSaturation_plot.r'),
                os.path.join(rseqc_dir, '{sample}.mismatch_profile.xls'),
                os.path.join(rseqc_dir, '{sample}.read_distribution.txt'),
                os.path.join(rseqc_dir, '{sample}.pos.DupRate.xls'),
                os.path.join(rseqc_dir, '{sample}.seq.DupRate.xls'),
                os.path.join(rseqc_dir, '{sample}.GC.xls'),
                os.path.join(rseqc_dir, '{sample}.NVC.xls'),
                os.path.join(rseqc_dir, '{sample}.qual.r'),
                os.path.join(rseqc_dir, '{sample}.RNA_fragment_size.txt'),
                os.path.join(rseqc_dir, '{sample}.STAR.genome.sorted.summary.txt'),
                os.path.join(fastq_dir, '{sample}_R1_fastq.html'),
                os.path.join(fastq_dir, '{sample}_R2_fastq.html')],
                sample=samples['samples'])

Answer 1

是的，这个是我 "myself" 想出来的。魔法在 "rule all" 部分。

这种规则组合有效：

reads = ['1', '2']

# Rule all is a pseudo-rule that tells snakemake what final files to generate.
rule all:
    input:
        expand([os.path.join(analyzed_dir, '{sample}.genes.results'),
                os.path.join(rseqc_dir, '{sample}.bam_stat.txt'),
                os.path.join(rseqc_dir, '{sample}.clipping_profile.xls'),
                os.path.join(rseqc_dir, '{sample}.deletion_profile.txt'),
                os.path.join(rseqc_dir, '{sample}.infer_experiment.txt'),
                os.path.join(rseqc_dir, '{sample}.geneBodyCoverage.txt'),
                os.path.join(rseqc_dir, '{sample}.inner_distance.txt'),
                os.path.join(rseqc_dir, '{sample}.insertion_profile.xls'),
                os.path.join(rseqc_dir, '{sample}.junction.xls'),
                os.path.join(rseqc_dir, '{sample}.junctionSaturation_plot.r'),
                os.path.join(rseqc_dir, '{sample}.mismatch_profile.xls'),
                os.path.join(rseqc_dir, '{sample}.read_distribution.txt'),
                os.path.join(rseqc_dir, '{sample}.pos.DupRate.xls'),
                os.path.join(rseqc_dir, '{sample}.seq.DupRate.xls'),
                os.path.join(rseqc_dir, '{sample}.GC.xls'),
                os.path.join(rseqc_dir, '{sample}.NVC.xls'),
                os.path.join(rseqc_dir, '{sample}.qual.r'),
                os.path.join(rseqc_dir, '{sample}.RNA_fragment_size.txt'),
                os.path.join(rseqc_dir, '{sample}.STAR.genome.sorted.summary.txt'),
                os.path.join(fastqc_dir, '{sample}_R{read}_001_fastqc.html')],
                sample=samples['samples'], read=reads)

注意简单地将 {read} 添加到其他相同的 fastqc 部分和顶部的定义或 "reads"（样本是标准样本列表）。

我使用这个fastqc规则：

rule fastqc:
    input:
       os.path.join(fastq_dir, '{sample}_R{read}_001.fastq.gz')
    output:
        os.path.join(fastqc_dir, '{sample}_R{read}_001_fastqc.html')
    conda:
        "../envs/fastqc.yaml"
    shell:
        '''
        #!/bin/bash
        fastqc {input} --outdir={fastqc_dir}
        '''

它与 "rule all" 有相同的行（像往常一样）。这行得通，谢谢大家的支持，释放吧。

使用 Snakemake 的 fastqc

fastqc using Snakemake

bioinformatics

snakemake