在 snakemake 运行 期间动态减少输入文件集
Reduce the set of input files dynamically during a snakemake run
这更多是关于 snakemake 功能的技术问题。我想知道是否可以在 snakemake 运行 期间动态更改输入样本集。
我想这样做的原因如下:让我们假设一组样本相关的 bam 文件。第一条规则确定每个样本的质量(基于 bam 文件),即涉及所有输入文件。
但是,根据指定的标准,只有一部分样本被认为是有效的,应该进一步处理。所以下一步(例如基因计数或其他)应该只对批准的 bam 文件进行,如下面的最小示例所示:
configfile: "config.yaml"
rule all:
input: "results/gene_count.tsv"
rule a:
input: expand( "data/{sample}.bam", sample=config['samples'])
output: "results/list_of_qual_approved_samples.out"
shell: '''command'''
rule b:
input: expand( "data/{sample}.bam", sample=config['valid_samples'])
output: "results/gene_count.tsv"
shell: '''command'''
在此示例中,规则 a 将使用有效示例名称列表扩展配置文件,即使我相信知道这是不可能的。
当然,直接的解决方案是有两个不同的输入:1.) 所有 bam 文件和 2.) 列出所有有效文件的文件。这将归结为在规则代码中进行样本选择。
rule alternative_b:
input:
expand( "data/{sample}.bam", sample=config['samples']),
"results/list_of_qual_approved_samples.out"
output: "results/gene_count.tsv"
shell: '''command'''
但是,您是否找到一种方法来设置规则以实现第一个示例的行为?
非常感谢,
拉尔夫
**这不完全是对您问题的回答,而是对实现目标的建议。 **
我认为在管道 运行 期间修改 yaml 文件是不可能的,或者至少不是微不足道的。
就个人而言,当我 运行 snakemake 工作流程时,我使用我称为 "metadata" 的外部文件。它们包括一个配置文件,还有一个包含样本列表的选项卡文件(可能还有关于所述样本的其他信息)。配置文件包含一个参数,该参数是该文件的路径。
在这样的设置中,我建议让您的 "rule a" 输出另一个包含所选示例的选项卡文件,并且该文件的路径可以包含在配置文件中(即使它没有启动工作流时存在)。规则 b 将该文件作为输入。
在你的情况下你可以:
config:
samples: "/path/to/samples.tab"
valid_samples: "/path/to/valid_samples.tab"
我不知道这是否有意义,因为它基于我自己的组织。我认为它很有用,因为它允许存储更多信息,而不仅仅是样本名称,而且如果您有 100 多个样本,管理起来会容易得多!
我想我有一个有趣的答案。
起初我认为这是不可能的。因为 Snakemake needs the final files
最后。所以不能一开始就不知道分开一组文件。
但后来我尝试使用函数 dynamic。使用 dynamic 函数,您不必知道规则将创建的文件数量。
所以我这样编码:
rule all:
input: "results/gene_count.tsv"
rule a:
input: expand( "data/{sample}.bam", sample=config['samples'])
output: dynamic("data2/{foo}.bam")
shell:
'./bloup.sh "{input}"'
rule b:
input: dynamic("data2/{foo}.bam")
output: touch("results/gene_count.tsv")
shell: '''command'''
就像你的第一个例子一样,snakefile 想要生成一个名为 results/gene_count.ts
的文件。
rule a
从配置文件中获取所有样本。此规则执行一个脚本,该脚本选择要创建的文件。我有 4 initial
个文件(geneA、geneB、geneC、geneD),它 only touches two
用于第二个目录中的输出(geneA 和 geneD 文件)。 dynamic函数没有问题
rule b
获取了rule a
创建的所有动态文件。所以你只需要生成 results/gene_count.tsv。我只是在示例中触摸它。
这里是log of Snakemake
了解更多信息:
Provided cores: 1
Rules claiming more threads will be scaled down.
Job counts:
count jobs
1 a
1 all
1 b
3
rule a:
input: data/geneA.bam, data/geneB.bam, data/geneC.bam, data/geneD.bam
output: data2/{*}.bam (dynamic)
Subsequent jobs will be added dynamically depending on the output of this rule
./bloup.sh "data/geneA.bam data/geneB.bam data/geneC.bam data/geneD.bam"
Dynamically updating jobs
Updating job b.
1 of 3 steps (33%) done
rule b:
input: data2/geneD.bam, data2/geneA.bam
output: results/gene_count.tsv
command
Touching output file results/gene_count.tsv.
2 of 3 steps (67%) done
localrule all:
input: results/gene_count.tsv
3 of 3 steps (100%) done
另一种方法,不使用 "dynamic"。
并不是您不知道要使用多少文件,而是您只使用了开始时使用的文件的一个子集。由于您能够生成所有潜在文件的 "samples.txt" 列表,我假设您有一个坚定的起点。
我做了类似的事情,我有一些初始文件,我想处理这些文件的有效性(在我的例子中,我正在提高质量~排序、索引等)。然后我想忽略除我的结果文件之外的所有内容。
为了避免创建第二个示例文件列表,我建议创建第二个数据目录 (reBamDIR)、data2 (BamDIR)。在 data2 中,您对所有有效文件进行符号链接。这样,Snake 就可以处理 data2 目录中的所有内容。使管道更容易向下移动,管道可以停止依赖示例列表,并且它可以使用通配符处理所有内容(更容易编码)。这是可能的,因为当我使用符号链接然后标准化名称时。我在输出规则中列出了符号链接文件,以便 Snakemake 了解它们,然后它可以创建 DAG。
`-- output
|-- bam
| |-- Pfeiffer2.bam -> /home/tboyarski/share/projects/tboyarski/gitRepo-LCR-BCCRC/Snakemake/buildArea/output/reBam/Pfeiffer2_realigned_sorted.bam
| `-- Pfeiffer2.bam.bai -> /home/tboyarski/share/projects/tboyarski/gitRepo-LCR- BCCRC/Snakemake/buildArea/output/reBam/Pfeiffer2_realigned_sorted.bam.bai
|-- fastq
|-- mPile
|-- reBam
| |-- Pfeiffer2_realigned_sorted.bam
| `-- Pfeiffer2_realigned_sorted.bam.bai
在这种情况下,您只需要 "validator" 中的一个 return 值,以及一个响应它的条件运算符。
我认为你已经在某个地方有了这个,因为你必须在验证步骤中使用条件。与其使用它将文件名写入 txt 文件,不如将文件符号链接到最终位置并继续。
我的原始数据在 reBamDIR 中。
我存储在 BamDIR 中的最终数据。
我只将管道中此阶段的文件符号链接到 bamDIR。
reBamDIR 中还有其他文件,但我不希望管道的其余部分看到它们,因此,我将它们过滤掉。
我不太确定如何实现 "validator" 和你的条件,因为我不知道你的情况,我也在学习。只是试图提供替代观点//方法。
from time import gmtime, strftime
rule indexBAM:
input:
expand("{outputDIR}/{reBamDIR}/{{samples}}{fileTAG}.bam", outputDIR=config["outputDIR"], reBamDIR=config["reBamDIR"], fileTAG=config["fileTAG"])
output:
expand("{outputDIR}/{reBamDIR}/{{samples}}{fileTAG}.bam.bai", outputDIR=config["outputDIR"], reBamDIR=config["reBamDIR"], fileTAG=config["fileTAG"]),
expand("{outputDIR}/{bamDIR}/{{samples}}.bam.bai", outputDIR=config["outputDIR"], bamDIR=config["bamDIR"]),
expand("{outputDIR}/{bamDIR}/{{samples}}.bam", outputDIR=config["outputDIR"], bamDIR=config["bamDIR"])
params:
bamDIR=config["bamDIR"],
outputDIR=config["outputDIR"],
logNAME="indexBAM." + strftime("%Y-%m-%d.%H-%M-%S", gmtime())
log:
"log/" + config["reBamDIR"]
shell:
"samtools index {input} {output[0]} " \
" 2> {log}/{params.logNAME}.stderr " \
"&& ln -fs $(pwd)/{output[0]} $(pwd)/{params.outputDIR}/{params.bamDIR}/{wildcards.samples}.bam.bai " \
"&& ln -fs $(pwd)/{input} $(pwd)/{params.outputDIR}/{params.bamDIR}/{wildcards.samples}.bam"
这更多是关于 snakemake 功能的技术问题。我想知道是否可以在 snakemake 运行 期间动态更改输入样本集。
我想这样做的原因如下:让我们假设一组样本相关的 bam 文件。第一条规则确定每个样本的质量(基于 bam 文件),即涉及所有输入文件。 但是,根据指定的标准,只有一部分样本被认为是有效的,应该进一步处理。所以下一步(例如基因计数或其他)应该只对批准的 bam 文件进行,如下面的最小示例所示:
configfile: "config.yaml"
rule all:
input: "results/gene_count.tsv"
rule a:
input: expand( "data/{sample}.bam", sample=config['samples'])
output: "results/list_of_qual_approved_samples.out"
shell: '''command'''
rule b:
input: expand( "data/{sample}.bam", sample=config['valid_samples'])
output: "results/gene_count.tsv"
shell: '''command'''
在此示例中,规则 a 将使用有效示例名称列表扩展配置文件,即使我相信知道这是不可能的。
当然,直接的解决方案是有两个不同的输入:1.) 所有 bam 文件和 2.) 列出所有有效文件的文件。这将归结为在规则代码中进行样本选择。
rule alternative_b:
input:
expand( "data/{sample}.bam", sample=config['samples']),
"results/list_of_qual_approved_samples.out"
output: "results/gene_count.tsv"
shell: '''command'''
但是,您是否找到一种方法来设置规则以实现第一个示例的行为?
非常感谢, 拉尔夫
**这不完全是对您问题的回答,而是对实现目标的建议。 **
我认为在管道 运行 期间修改 yaml 文件是不可能的,或者至少不是微不足道的。
就个人而言,当我 运行 snakemake 工作流程时,我使用我称为 "metadata" 的外部文件。它们包括一个配置文件,还有一个包含样本列表的选项卡文件(可能还有关于所述样本的其他信息)。配置文件包含一个参数,该参数是该文件的路径。
在这样的设置中,我建议让您的 "rule a" 输出另一个包含所选示例的选项卡文件,并且该文件的路径可以包含在配置文件中(即使它没有启动工作流时存在)。规则 b 将该文件作为输入。
在你的情况下你可以:
config:
samples: "/path/to/samples.tab"
valid_samples: "/path/to/valid_samples.tab"
我不知道这是否有意义,因为它基于我自己的组织。我认为它很有用,因为它允许存储更多信息,而不仅仅是样本名称,而且如果您有 100 多个样本,管理起来会容易得多!
我想我有一个有趣的答案。
起初我认为这是不可能的。因为 Snakemake needs the final files
最后。所以不能一开始就不知道分开一组文件。
但后来我尝试使用函数 dynamic。使用 dynamic 函数,您不必知道规则将创建的文件数量。
所以我这样编码:
rule all:
input: "results/gene_count.tsv"
rule a:
input: expand( "data/{sample}.bam", sample=config['samples'])
output: dynamic("data2/{foo}.bam")
shell:
'./bloup.sh "{input}"'
rule b:
input: dynamic("data2/{foo}.bam")
output: touch("results/gene_count.tsv")
shell: '''command'''
就像你的第一个例子一样,snakefile 想要生成一个名为 results/gene_count.ts
的文件。
rule a
从配置文件中获取所有样本。此规则执行一个脚本,该脚本选择要创建的文件。我有 4 initial
个文件(geneA、geneB、geneC、geneD),它 only touches two
用于第二个目录中的输出(geneA 和 geneD 文件)。 dynamic函数没有问题
rule b
获取了rule a
创建的所有动态文件。所以你只需要生成 results/gene_count.tsv。我只是在示例中触摸它。
这里是log of Snakemake
了解更多信息:
Provided cores: 1
Rules claiming more threads will be scaled down.
Job counts:
count jobs
1 a
1 all
1 b
3
rule a:
input: data/geneA.bam, data/geneB.bam, data/geneC.bam, data/geneD.bam
output: data2/{*}.bam (dynamic)
Subsequent jobs will be added dynamically depending on the output of this rule
./bloup.sh "data/geneA.bam data/geneB.bam data/geneC.bam data/geneD.bam"
Dynamically updating jobs
Updating job b.
1 of 3 steps (33%) done
rule b:
input: data2/geneD.bam, data2/geneA.bam
output: results/gene_count.tsv
command
Touching output file results/gene_count.tsv.
2 of 3 steps (67%) done
localrule all:
input: results/gene_count.tsv
3 of 3 steps (100%) done
另一种方法,不使用 "dynamic"。
并不是您不知道要使用多少文件,而是您只使用了开始时使用的文件的一个子集。由于您能够生成所有潜在文件的 "samples.txt" 列表,我假设您有一个坚定的起点。
我做了类似的事情,我有一些初始文件,我想处理这些文件的有效性(在我的例子中,我正在提高质量~排序、索引等)。然后我想忽略除我的结果文件之外的所有内容。
为了避免创建第二个示例文件列表,我建议创建第二个数据目录 (reBamDIR)、data2 (BamDIR)。在 data2 中,您对所有有效文件进行符号链接。这样,Snake 就可以处理 data2 目录中的所有内容。使管道更容易向下移动,管道可以停止依赖示例列表,并且它可以使用通配符处理所有内容(更容易编码)。这是可能的,因为当我使用符号链接然后标准化名称时。我在输出规则中列出了符号链接文件,以便 Snakemake 了解它们,然后它可以创建 DAG。
`-- output
|-- bam
| |-- Pfeiffer2.bam -> /home/tboyarski/share/projects/tboyarski/gitRepo-LCR-BCCRC/Snakemake/buildArea/output/reBam/Pfeiffer2_realigned_sorted.bam
| `-- Pfeiffer2.bam.bai -> /home/tboyarski/share/projects/tboyarski/gitRepo-LCR- BCCRC/Snakemake/buildArea/output/reBam/Pfeiffer2_realigned_sorted.bam.bai
|-- fastq
|-- mPile
|-- reBam
| |-- Pfeiffer2_realigned_sorted.bam
| `-- Pfeiffer2_realigned_sorted.bam.bai
在这种情况下,您只需要 "validator" 中的一个 return 值,以及一个响应它的条件运算符。
我认为你已经在某个地方有了这个,因为你必须在验证步骤中使用条件。与其使用它将文件名写入 txt 文件,不如将文件符号链接到最终位置并继续。
我的原始数据在 reBamDIR 中。 我存储在 BamDIR 中的最终数据。 我只将管道中此阶段的文件符号链接到 bamDIR。 reBamDIR 中还有其他文件,但我不希望管道的其余部分看到它们,因此,我将它们过滤掉。
我不太确定如何实现 "validator" 和你的条件,因为我不知道你的情况,我也在学习。只是试图提供替代观点//方法。
from time import gmtime, strftime
rule indexBAM:
input:
expand("{outputDIR}/{reBamDIR}/{{samples}}{fileTAG}.bam", outputDIR=config["outputDIR"], reBamDIR=config["reBamDIR"], fileTAG=config["fileTAG"])
output:
expand("{outputDIR}/{reBamDIR}/{{samples}}{fileTAG}.bam.bai", outputDIR=config["outputDIR"], reBamDIR=config["reBamDIR"], fileTAG=config["fileTAG"]),
expand("{outputDIR}/{bamDIR}/{{samples}}.bam.bai", outputDIR=config["outputDIR"], bamDIR=config["bamDIR"]),
expand("{outputDIR}/{bamDIR}/{{samples}}.bam", outputDIR=config["outputDIR"], bamDIR=config["bamDIR"])
params:
bamDIR=config["bamDIR"],
outputDIR=config["outputDIR"],
logNAME="indexBAM." + strftime("%Y-%m-%d.%H-%M-%S", gmtime())
log:
"log/" + config["reBamDIR"]
shell:
"samtools index {input} {output[0]} " \
" 2> {log}/{params.logNAME}.stderr " \
"&& ln -fs $(pwd)/{output[0]} $(pwd)/{params.outputDIR}/{params.bamDIR}/{wildcards.samples}.bam.bai " \
"&& ln -fs $(pwd)/{input} $(pwd)/{params.outputDIR}/{params.bamDIR}/{wildcards.samples}.bam"