Snakemake Получение контрольной точки и агрегатной функции для работы - PullRequest
Новая
       16

Snakemake Получение контрольной точки и агрегатной функции для работы

0 голосов
/ 21 мая 2019

У меня проблема с тем, чтобы заставить работать мою команду агрегата 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:

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

1 Ответ

0 голосов
/ 21 мая 2019

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

...