Динамически устанавливаемые параметры в snakemake - PullRequest
0 голосов
/ 13 июля 2020

Я пытаюсь настроить конвейер snakemake, который снова выравнивает и читает эталонный геном. Правило, которое индексирует геном, имеет параметр, которому требуется максимальная длина считываний, которые будут сопоставлены с геномом. Сейчас я вручную установил это в файле 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
        """

1 Ответ

0 голосов
/ 14 июля 2020

Как уже советовал @ManavalanGajapathy, вы можете сделать это в своем сценарии оболочки: он уже достаточно сложен, поэтому туда можно добавить небольшой сценарий awk.

В любом случае, сложные сценарии оболочки не так популярны в Раздел Snakemake shell, поэтому вы все равно можете рассчитать эту длину в разделе params:

    params: 
        len=lambda wildcards: getMaxReadLength(wildcards)

Теперь вам нужно реализовать эту функцию getMaxReadLength, которая основана на значениях подстановочных знаков (если есть ) получает максимальную длину чтения:

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