在snakemake中动态设置参数

Question

我正在尝试建立一条蛇制造管道，它可以再次比对读取一个参考基因组。索引基因组的规则有一个参数，该参数需要与基因组映射的读取的最大长度。现在，我已经在config.yaml文件中手动设置了这个值，但是如果能根据管道中输入读取的内容来设置这个值就太好了。例如，通过运行一个小的awk代码片段来查找所有fastq文件中的最大长度。

管道可以接受任何给定数量的样本来与基因组进行映射，但只有样本中最长的读取长度应该用于构建索引(我宁愿只使用一次)

有什么方法可以做到这一点吗？

我当前的规则如下所示：

rule index_reference_genome:
    """
    Index reference genome.
    """
    input: 
        ref_fasta=config['ref_fasta'],
        ref_gtf=config['ref_gtf']
    output: 
        directory("{outdir}/references/star")
    threads: 40
    params: 
        len=int(config['max_read_len'])-1,
        #star=config['star']
    log: "{outdir}/logs/star/index_reference_genome.log"
    benchmark: "{outdir}/benchmarks/star/index_reference_genome.log"
    conda: config['conda_env']
    shell: 
        """
        if [ ! -d {output} ] ; then mkdir {output}; fi;
        STAR \
            --runThreadN {threads} \
            --runMode genomeGenerate \
            --genomeDir {output} \
            --genomeFastaFiles {input.ref_fasta} \
            --sjdbGTFfile {input.ref_gtf} \
            --sjdbOverhang {params.len} &> {log}; #--outTmpDir {output}/STAR_tmp
        """

Answer 1

正如@ManavalanGajapathy已经建议的那样，你可以在你的shell脚本中这样做:它已经足够复杂了，所以可以在那里添加一个小的awk脚本。

无论如何，复杂的shell脚本在Snakemake的params部分中并不是很流行，所以您仍然可以在shell部分中计算这个长度：

    params: 
        len=lambda wildcards: getMaxReadLength(wildcards)

现在您需要实现这个getMaxReadLength函数，它基于通配符(如果有)的值来检索读取的最大长度：

def getMaxReadLength(wildcards):
    return "".join(shell("Your AWK code here", iterable=True))