Snakemake: Как обрабатывать динамически созданные файлы с динамически создаваемыми именами файлов? - PullRequest
Новая
       122

Snakemake: Как обрабатывать динамически созданные файлы с динамически создаваемыми именами файлов?

1 голос
/ 14 июля 2020

У меня есть файл multifasta DE.faa, содержащий неизвестное количество белковых последовательностей и заголовки fasta, например GB012883-PA, GB009065-PA, GB007275-PA. Я разделил файл multifasta (используя правило gather_proteins_with_no_BL62_hits) на отдельные файлы fasta с именами файлов: GB012883-PA.faa, GB009065-PA.faa, GB007275-PA.faa, и теперь я sh выполняю BLASTp для каждого из них удаленно.

На выходе из BLASTp я sh должен иметь файлы tsv с соответствующими именами файлов GB012883-PA.tsv, GB009065-PA.tsv, GB007275-PA.tsv и файлы стратегии поиска: GB012883-PA.asn, GB009065-PA.asn, GB007275-PA.asn.

Вот как далеко я продвинулся: 

rule gather_proteins_with_no_BL62_hits:
    input:
        all_found_prot_seqs = "DE.faa"
    output:
        directory("no_BL62_hits/"), 
        touch("fasta_expansion.done")
    script:
        "gather_proteins_with_no_BL62_hits.R"


checkpoint BLAST_missing_protein_seqs:
    input:
        flag = "fasta_expansion.done",
        seqs = "no_BL62_hits/{seq_name}.faa"
    output:
        res = "BL45_headerless/{seq_name}.tsv",
        s_t = "BL45_headerless/{seq_name}.asn"
    shell:
        r"""
        blastp -remote \
        -query {input.seqs} \
        -evalue 0.01 \
        -db nr \
        -export_search_strategy {output.s_t} \
        -word_size 2 \
        -matrix "BLOSUM45"
        -max_target_seqs 100 \
        -outfmt '6 sseqid sacc qstart' > {output.res} 2> {log}
        """


def aggregate_input(wildcards):
    with open(checkpoints.BLAST_missing_protein_seqs.get(**wildcards).output.res, 'r') as f:
        return [seq_name.rstrip() + '.tsv' for seq_name in f]


rule aggregate:
    input:
        aggregate_input
    output:
        touch("BLAST45.done")

Однако, когда я запускаю это, я получаю InputFunctionException с сообщением: WorkflowError: Missing wildcard values for seq_name.

Мой вопрос:

Как мне сообщить Snakemake, что подстановочный знак seq_name относится к идентификаторам последовательности GB012883-PA, GB009065-PA, GB007275-PA (и исправить ошибку WorkflowError)?

1 Ответ

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

Надеюсь, вы уже сами это решили. Если нет, то я надеюсь, что это поможет. Я не тестировал ничего из этого, поэтому вам может потребоваться исправить это здесь и там: 

checkpoint gather_proteins_with_no_BL62_hits:
    input:
        all_found_prot_seqs = "DE.faa"
    output:
        directory("no_BL62_hits/"), 
        touch("fasta_expansion.done")
    script:
        "gather_proteins_with_no_BL62_hits.R"


rule BLAST_missing_protein_seqs:
    input:
        flag = "fasta_expansion.done",
        seqs = "no_BL62_hits/{seq_name}.faa"
    output:
        res = "BL45_headerless/{seq_name}.tsv",
        s_t = "BL45_headerless/{seq_name}.asn"
    shell:
        r"""
        blastp -remote \
        -query {input.seqs} \
        -evalue 0.01 \
        -db nr \
        -export_search_strategy {output.s_t} \
        -word_size 2 \
        -matrix "BLOSUM45"
        -max_target_seqs 100 \
        -outfmt '6 sseqid sacc qstart' > {output.res} 2> {log}
        """


def aggregate_input(wildcards):
    checkpoint_output = checkpoints.gather_proteins_with_no_BL62_hits.get(**wildcards).output[0]
    seq_name, = glob_wildcards(checkpoint_output + "{seq_name}.faa")
    return expand("BL45_headerless/{seq_name}.{ext}", seq_name=seq_name, ext=["tsv", "asn"])


rule aggregate:
    input:
        aggregate_input
    output:
        touch("BLAST45.done")

Я изменил, какое правило является контрольной точкой. Мы можем знать / предполагать, каким должен быть наш результат после выполнения gather_proteins_with_no_BL62_hits, а не после выполнения BLAST_missing_protein_seqs, поскольку он уже предполагает, что мы знаем, какой результат мы хотим. Итак, после выполнения gather_proteins_with_no_BL62_hits мы запускаем aggregate_input. Что должно найти все seq_names, которые вам нужны.

...