bioinformatics - 使用 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 脚本中。

我应该如何正确地“循环”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

是的，这个是我“自己”想出来的。神奇之处在于“统治一切”部分。

这种规则组合有效:

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 部分和顶部的定义或“读取”(样本是标准样本列表)。

我使用这个 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}
        '''

它与“全部规则”(如往常一样)具有相同的行。这行得通，感谢大家的支持，释放吧。