Snakemake заставляет контрольную точку и агрегатную функцию работать

У меня проблема с работой моей совокупной команды snakemake. Я надеюсь взять данный файл GTF, найти отдельные регионы внутри GTF и, если они найдутся, записать эти регионы в отдельный файл. Таким образом, я не уверен, сколько выходных файлов GTF будет создавать каждый входной файл GTF. Чтобы решить эту проблему, я пытаюсь использовать контрольную точку змейки.

Для этого я написал краткий сценарий под названием collapse_gtf_file.py, который просто принимает файл GTF и генерирует N файлов, соответствующих количеству найденных отдельных регионов. Таким образом, если дан файл test_regions.gtf с тремя областями, он сгенерирует test_regions_1.gtf, test_regions_2.gtf test_regions_3.gtf соответственно.

После указанного разделения все файлы GTF должны быть преобразованы в файлы fasta и объединены.

Однако мне не удалось заставить мою команду контрольной точки работать. Я могу заставить примеры работать, но когда я пытаюсь построить более крупный конвейер вокруг этой контрольной точки, он ломается.

До сих пор я пробовал следовать руководству по контрольным точкам, которое можно найти здесь https://snakemake.readthedocs.io/en/stable/snakefile/rules.html#dynamic-files

sample=config["samples"]
reference=config["reference"]

rule all:
    input:
    ¦   expand("string_tie_assembly/{sample}.gtf", sample=sample),
    ¦                                                                                                 expand("string_tie_assembly_merged/merged_{sample}.gtf",sample=sample),
    ¦   #expand("split_gtf_file/{sample}", sample=sample),
    ¦   #expand("lncRNA_fasta/{sample}.fa", sample=sample),
    ¦   "combined_fasta/all_fastas_combined.fa",
    ¦   "run_CPC2/cpc_calc_output.txt"

rule samtools_sort:
    input:
    ¦   "mapped_reads/{sample}.bam"
    output:
    ¦   "sorted_reads/{sample}.sorted.bam"
    shell:
    ¦   "samtools sort -T sorted_reads/{wildcards.sample} {input} > {output}"

rule samtools_index:
    input:
        "sorted_reads/{sample}.sorted.bam"
    output:
        "sorted_reads/{sample}.sorted.bam.bai"
    shell:
        "samtools index {input}"

rule generate_fastq:
    input:
        "sorted_reads/{sample}.sorted.bam"
    output:
        "fastq_reads/{sample}.fastq"
    shell:
        "samtools fastq {input} > {output}"

rule string_tie_assembly:
    input:
        "sorted_reads/{sample}.sorted.bam"
    output:
        "string_tie_assembly/{sample}.gtf"
    shell:
        "stringtie {input} -f 0.0 -a 0 -m 50 -c 3.0 -f 0.0 -o {output}"

rule merge_gtf_file_features:
    input:
    ¦   "string_tie_assembly/{sample}.gtf"
    output:
    ¦   "string_tie_assembly_merged/merged_{sample}.gtf"
    shell:
    ¦   #prevents errors when there's no sequence
    ¦   """
    ¦   set +e
    ¦   stringtie --merge -o {output} -m 25 -c 3.0 {input}
    ¦   exitcode=$?
    ¦   if [ $exitcode -eq 1 ]
    ¦   then
    ¦   ¦   exit 0
    ¦   else
    ¦   ¦   exit 0
    ¦   fi
    ¦   """

#This is where the issue seems to arise from. Modeled after https://snakemake.readthedocs.io/en/stable/snakefile/rules.html#dynamic-files 

checkpoint clustering:
    input:
    ¦   "string_tie_assembly_merged/merged_{sample}.gtf"
    output:
    ¦   clusters = directory("split_gtf_file/{sample}")
    shell:
    ¦   """
    ¦   mkdir -p split_gtf_file/{wildcards.sample} ;

    ¦   python collapse_gtf_file.py -gtf {input} -o split_gtf_file/{wildcards.sample}/{wildcards.sample}
    ¦   """

rule gtf_to_fasta:
    input:
    ¦   "split_gtf_file/{sample}/{sample}_{i}.gtf"
    output:
    ¦   "lncRNA_fasta/{sample}/canidate_{sample}_{i}.fa"
    wildcard_constraints:
    ¦   i="\d+"
    shell:
    ¦   "gffread -w {output} -g {reference} {input}"


def aggregate_input(wildcards):
    checkpoint_output = checkpoints.clustering.get(**wildcards).output[0]
    x = expand("lncRNA_fasta/{sample}/canidate_{sample}_{i}.fa",
    ¦   sample=wildcards.sample,
    ¦   i=glob_wildcards(os.path.join(checkpoint_output, "{i}.fa")).i)
    return x

rule combine_fasta_file:
    input:
    ¦   aggregate_input
    output:
    ¦   "combined_fasta/all_fastas_combined.fa"
    shell:
    ¦   "cat {input} > {output}"


Error Message:

InputFunctionException in combine_fasta_file:
WorkflowError: Missing wildcard values for sample
Wildcards:

Мне кажется, что это указывает на то, что что-то не так с тем, как я назвал подстановочные знаки выше в совокупной команде, но я не могу понять, что именно. Любые указатели будут очень признательны.


person pabster212    schedule 21.05.2019    source источник


Ответы (2)


Функция aggregate_input ожидает переменную wildcards.sample, но у вас нет подстановочных знаков, указанных в rule combine_fasta_file. Вы можете либо указать подстановочный знак в этом правиле, либо реорганизовать функцию для использования глобальной переменной, если это применимо.

person Manavalan Gajapathy    schedule 21.05.2019
comment
Спасибо за ваш вклад, но я не думаю, что это причина. На основе найденного кода snakemake.readthedocs.io/en/stable / snakefiles / - похоже, что команда, объединяющая все файлы (здесь combine_fasta_file), не должна нуждаться в подстановочных знаках. Или я неправильно понимаю вашу точку зрения? - person pabster212; 21.05.2019
comment
Ваша ссылка примерно dynamic(), но она не используется в вашем коде. Вы имеете в виду checkpoint вместо? - person Manavalan Gajapathy; 21.05.2019
comment
Да! Моя ошибка, вот правильная ссылка. snakemake.readthedocs.io/en/stable/ snakefiles /. Я также обновил вопрос примером, который мне пришлось запускать изолированно. - person pabster212; 21.05.2019
comment
Удалил тестовый пример. Это был всего лишь пример того, что мне довелось поработать, и это было очень похоже. По-прежнему возникает та же проблема / ошибка - person pabster212; 21.05.2019
comment
Вы уверены, что отсутствие подстановочных знаков в output не проблема? Код в вашей ссылке о контрольной точке зависит от подстановочного знака sample в rule aggregate. - person Manavalan Gajapathy; 21.05.2019
comment
Да!! Это было. Итак, я думаю, что мою проблему можно резюмировать в том, что я пытался использовать свой rule combine_fasta_file для объединения нескольких шагов. Вместо этого мне нужно было сначала собрать все .fa файлы, а ЗАТЕМ собрать все .fa файлы из разных образцов. Огромное спасибо!! - person pabster212; 21.05.2019

Правило, которое собирает выходные данные контрольной точки (combine_fasta_file, в вашем случае), должно содержать подстановочные знаки, присутствующие в выходных данных правила контрольной точки (clustering, в вашем случае). Это связано с тем, что Snakemake смотрит на желаемый результат, а затем пытается вернуться через ваши правила к доступным входным файлам. Если подстановочные знаки в правиле сбора не включают подстановочные знаки правила контрольной точки, точная ссылка на контрольную точку может быть неоднозначной.

Чтобы уточнить комментарии к контрольной точке Snakemake docs пример Манавалана Гаджапати и многочисленные агрегации, необходимые для pabster212, я создал минимальный пример для выделения контрольных точек и множественного сбора.

Настройка входа:

mkdir -p input/sampleA_run{1,2,3}
mkdir -p input/sampleB_run{4,5}
touch input/sampleA_run1/{apple,banana,cucumber,beet,basil}.txt
touch input/sampleA_run2/{alpha,beta,gamma}.txt
touch input/sampleA_run3/{red,blue,green,brown}.txt
touch input/sampleB_run4/{alice,betty,beth,barbara,christy}.txt
touch input/sampleB_run5/{alex,brad,bart,bennett,clive}.txt

Snakefile, который проходит через контрольные точки и мультиагрегацию:

rule all:
    input:
        # Final
        expand("output/gathered_samples/{sample}.txt", sample = ["sampleA","sampleB"])

# My real-world analysis at this step (Nanopore guppy basecaller) 
# Takes as input a directory 
# With a inconsistent numbers of files and names for files (.fast5 in real analysis)
# Due to filtering, there are more inputs than outputs (.fastq in real analysis)
# And the output filenames are not predictably named or numbered
# Therefore I use a checkpoint to evaluate which output files have been created
# In this toy example, only files beginning with 'b' are maintained
checkpoint unpredictable_file_generator:
    input:
        "input/{sample}_{run_id}"
    output:
        directory("output/unpredictable_files/{sample}_{run_id}"),
    shell:
        """
        mkdir -p {output}

        for some_file in {input}/*
        do
            outfile=$(basename "$some_file")
            if [[ "$outfile" =~ ^b.* ]]
            then
                touch "{output}/$outfile"
            fi
        done
        """

# For my real-world analysis, this intermediate step 
# Represents generation of bams from fastqs
# No intermediate step is necessarily required following the checkpoint
# Just modify the input-gathering function below (gather_from_checkpoint) 
# To expand 'output/unpredictable_files' instead of 'output/intermediate_files'
rule intermediate_step:
    input:
        "output/unpredictable_files/{sample}_{run_id}/{file_id}.txt"
    output:
        "output/intermediate_files/{sample}_{run_id}/{file_id}.txt"
    shell:
        """
        cp {input} {output}
        """

# Function to gather all outputs from checkpoint rule
# The output of the rule referring back to the checkpoint (gather_one, below)
# Must contain the all wildcards in the checkpoint rule
# To gather further (e.g. gather by sample rather than by run_id)
# An additional gather rule is required
def gather_from_checkpoint(wildcards):
    checkpoint_output = checkpoints.unpredictable_file_generator.get(**wildcards).output[0]
    return expand("output/intermediate_files/{sample}_{run_id}/{file_id}.txt",
           sample=wildcards.sample,
           run_id=wildcards.run_id,
           file_id=glob_wildcards(os.path.join(checkpoint_output, "{file_id}.txt")).file_id)

# Gathers all files by run_id
# But each sample is still divided into runs
# For my real-world analysis, this could represent a 'samtools merge bam'
rule gather_one:
    input:
        gather_from_checkpoint
    output:
        "output/gathered_runs/{sample}_{run_id}.txt"
    shell:
        """
        # Generate a file for each run 
        # Containing the names of all files output by 'unknown_file_generator' 
        echo {input} | tr ' ' $'\\n' > {output}
        """

# Gathers all previously-run_id-gathered files by sample
# For my real-world analysis, this could represent a second 'samtools merge bam'
# After which point the first merge could potentially be removed if it were marked 'temporary'
############
# Note wildcard restriction was required 
# To prevent the {run_id} match from descending into diretories
# And generating 'ChildIOException' errors
rule gather_two:
    input:
        lambda wildcards: expand("output/gathered_runs/{sample}_{run_id}.txt",
            sample = wildcards.sample,
            run_id = glob_wildcards("input/" + wildcards.sample + "_{run_id,[^/]+}/").run_id)
    output:
        "output/gathered_samples/{sample}.txt"
    shell:
        """
        # This output should contain 6 filenames for each sample
        # Each beginning with 'b'
        cat {input} > {output}
        """
person GobiJerboa    schedule 17.04.2020