如何运行 snakemake 输入部分中的 AWK 脚本（以检索 snakemake 输入变量）？

Question

我正在尝试构建用于分析 RNAseq 数据集的 snakemake 管道。一些数据集可能是 (1) 单端或双端，以及 (2) 数据集之间的搁浅也可能不同。我希望 snakemake 管道自动检测这 2 个变量并将其填充到 RSEM 输入中。

我想出了一个方法，使用来自 RSeQC 的 infer_experiment.py 输出文件“strandedness.txt”（注意前两行是空的）。

This is PairEnd Data Fraction of reads failed to determine: 0.0134 Fraction of reads explained by "1++,1--,2+-,2-+": 0.0049 Fraction of reads explained by "1+-,1-+,2++,2--": 0.9817

现在我想使用 AWK，以便使用以下 AWK 代码从文件“strandedness.txt”中提取例如“PairEnd”，它在 bash 中运行良好：

awk 'FNR == 3 {print }' ./strandedness.txt

正确输出“PairEnd”

我把代码放在这个 AWK 脚本中 "pairend.awk":

#!/bin/awk -f
FNR == 3 {print }

如果我单独调用它，它工作正常：

./pairend.awk ./strandedness.txt

现在我想调用我的 snakemake 文件中的脚本，但我不知道该怎么做。 google 上的大多数示例在 shell 中使用 AWK 脚本，但我实际上想使用该脚本输出一个变量，该变量将用作我的 snakemake 管道的输入。我试过了，但无法正常工作：

DATASET,SAMPLE,FRR, =glob_wildcards(config["project_path"]+"resources/raw_datasets/{dataset}/{sample}_Homo_sapiens_RNA-Seq_{frr}.fastq.gz")
print(DATASET,SAMPLE,FRR)

rule all:
        input:
                expand(config["project_path"]+"results/{dataset}/RSEM/{sample}_strandedness.txt",dataset=DATASET, sample=SAMPLE)

## Transcript quantification
rule RSEM_calculate_expression:
        input:
                fasta=config["project_path"]+"resources/Homo_sapiens.GRCh38.dna.primary_assembly.fa",
                gtf=config["project_path"]+"resources/Homo_sapiens.GRCh38.106.gtf",
                bam=config["project_path"]+"results/{dataset}/star_aligned_2pass/{sample}_Aligned.sortedByCoord.out.bam",
                paired_end=config["project_path"]+"workflow/pairend.awk" config["project_path"]+"results/{dataset}/RSEM/{sample}_strandedness.txt"
        output:
                gene_results=config["project_path"]+"results/{dataset}/RSEM/{sample}.genes.results",
                isoform_results=config["project_path"]+"results/{dataset}/RSEM/{sample}.isoforms.results"
        params:
                index=config["project_path"]+"resources/starIndex/"
        threads:
                12
        log:
                config["project_path"]+"results/{dataset}/RSEM/{sample}.log"
        shell:
                """
                rsem-calculate-expression --star --num-threads {threads} {paired_end} --alignments {input.bam} {input.index} {wildcards.sample} {log}
                """

我也试过从 AWK 脚本调用 AWK 代码但没有成功：

Error: SyntaxError in line ... of snakefile: invalid syntax

我划线，但它指向 paired_end=config["project_path"]+"workflow/pairend.awk" config["project_path "]+"results/{dataset}/RSEM/{sample}_strandedness.txt"

如何在 snakemake 的输入部分正确地运行 my/an AWK 脚本（我对运行在 shell 部分中使用它不感兴趣。任何帮助将不胜感激！

Answer 1

因为您要查找的是字符串而不是文件，所以应将其放入参数而不是输入。您也许可以使用 shell 函数从输入函数调用 awk，但我只在 run 指令中看到过它。

我认为最简单的选择是

将您的 awk 命令转换为 python 中的输入函数。

def get_pairend(wildcards, input):
    # since your file is small, just read the whole thing.  If it's larger don't read everything.
    # get third line, third token
    return open(input.paired_end).readlines()[2].split()[2]

rule RSEM_calculate_expression:
    input:
        ...
        paired_end=config["project_path"]+"results/{dataset}/RSEM/{sample}_strandedness.txt"
        ...
    params:
        index=config["project_path"]+"resources/starIndex/",
        paired_end=get_pairend,
    shell:
        """
        rsem-calculate-expression --star --num-threads {threads} {params.paired_end} --alignments {input.bam} {input.index} {wildcards.sample} {log}
        """

运行 awk 在你的 shell 指令中

rule RSEM_calculate_expression:
    input:
        ...
        paired_end=config["project_path"]+"results/{dataset}/RSEM/{sample}_strandedness.txt"

    shell:
        """
        paired=$(awk 'FNR == 3 {{print }}' {input.paired_end})
        rsem-calculate-expression --star --num-threads {threads} $paired --alignments {input.bam} {input.index} {wildcards.sample} {log}
        """

我更喜欢第一个选项。如果您决定做一些更复杂的事情并且您不需要额外的 awk 知识来理解您在做什么，那么修改起来会更容易。

如何 运行 snakemake 输入部分中的 AWK 脚本（以检索 snakemake 输入变量）？

How to run an AWK script within the input section of snakemake (in order to retrieve a snakemake input variable)?

bash

awk

snakemake

如何运行 snakemake 输入部分中的 AWK 脚本（以检索 snakemake 输入变量）？