我最近開始使用 Snakemake。我正在構建一個管道，其目標是分離讀取的分數（基于 Kraken2 分類的 taxID）。我無法理解如何根據 taxID 串列引入新的通配符。
我有：
- 我想進一步檢查的 taxID 串列。
- 基于 taxID 和 Kraken2 分類從我的原始 .fastq 檔案中提取讀取的腳本
下面顯示的代碼是管道的一個子集，但可以解決問題。R 腳本 extract_microbial_taxids.R 生成一個檔案，其中每個 taxID 位于單獨的行中。函式 SplitGenusList() 將此檔案拆分為一個陣列，其中每個元素都是一個 taxID。
我想拆分輸入的 .fastq 檔案并將 taxID 作為新的通配符引入 - 可能基于 R 腳本中的 taxID 串列。關于我如何實作這一目標的任何線索？然后通配符將用于幾個規則，然后只有 {sample} 通配符再次相關。
rule all:
input:
expand("data/microbial_taxids/{sample}/genus_list.txt", sample = config["SAMPLE"])
rule extract_microbial_taxids:
input:
"data/kraken2_classification/{sample}/{sample}.report",
config["Rlibpath"]
output:
"data/microbial_taxids/{sample}/genus_list.txt"
script:
"code/extract_microbial_taxids.R"
def SplitGenusList(genus_list):
genus_list = open(genus_list)
contents = genus_list.read()
contents_split = contents.splitlines()
return contents_split
rule extract_microbial_reads:
input:
taxID_list = "data/microbial_taxids/{sample}/genus_list.txt",
kraken_file = "data/kraken2_classification/{sample}/{sample}.kraken",
reads = config["reads"] "{sample}.fastq"
output:
taxID_reads = expand("data/microbial_reads/{sample}/{taxID}.fastq", taxID = SplitGenusList(input.taxID_list))
shell:
"""
for taxID in genus_list
do
python code/extract_kraken_reads.py --kraken {input.kraken_file} -s {input.reads} -o "data/microbial_reads/{sample}/$taxID\_.fastq" -t $taxID --include-children
done
"""
更新
在@Oliver 和@dariober 的大量幫助下，我設法構建了一個按預期作業的作業流程。我從中學到了一些東西：
1：規則變數（如）對于 python 、 和 shell{input.taxID_list}的寫法應該不同。input.taxID_list{input.taxID_list}
2：運行命令允許以shell()一種看起來非常方便的方式結合python和shell（使用@dariober建議的功能）。
這是作業代碼：
rule all:
input:
expand("data/microbial_taxids/{sample}/genus_list.txt", sample = config["SAMPLE"]),
expand("data/microbial_reads/{sample}/extracting_microbial_reads.done", sample = config["SAMPLE"])
rule extract_microbial_taxids:
input:
"data/kraken2_classification/{sample}/{sample}.report",
config["Rlibpath"]
output:
"data/microbial_taxids/{sample}/genus_list.txt"
script:
"code/extract_microbial_taxids.R"
def SplitGenusList(genus_list):
genus_list = open(genus_list)
contents = genus_list.read()
contents_split = contents.splitlines()
return contents_split
rule extract_microbial_reads:
input:
taxID_list="data/microbial_taxids/{sample}/genus_list.txt",
kraken_file="data/kraken2_classification/{sample}/{sample}.kraken",
reads=config["reads"] "{sample}.fastq",
report_file="data/kraken2_classification/{sample}/{sample}.report"
output:
touch("data/microbial_reads/{sample}/extracting_microbial_reads.done")
run:
import os
genus_list = SplitGenusList(input.taxID_list)
for taxID in genus_list:
shell('python code/extract_kraken_reads.py -k {input.kraken_file} -s {input.reads} -o "data/microbial_reads/{wildcards.sample}/' taxID '_.fastq" -t ' taxID ' -r {input.report_file} --include-children')
uj5u.com熱心網友回復：
我不確定我是否正確理解了您，但本質上您只是在尋找一種在規則中參考變數的方法genus_list，該變數應該包含 taxID 或利息，對嗎？
好吧，我認為這里有一個主要問題，Snakemake 需要能夠匯出具體的輸出來創建 DAG，因為它決定了規則執行的順序。因此，即使您將taxID_reads = expand(...)其放入輸出中，Snakemake 仍然會嘗試找出實際結果。因此，SplitGenusList呼叫該函式，然后由于檔案不存在而產生錯誤genus_list。
因此，我建議不要在輸入或輸出中呼叫此函式，而只是在代碼塊中呼叫。為了記錄給定規則的成功執行，您可以簡單地創建一個虛擬檔案。
我的建議：
rule all:
input:
expand("data/microbial_reads/{sample}/extracting_microbial_reads.done", sample=config["SAMPLE"])
rule extract_microbial_taxids:
input:
"data/kraken2_classification/{sample}/{sample}.report",
config["Rlibpath"]
output:
"data/microbial_taxids/{sample}/genus_list.txt"
script:
"code/extract_microbial_taxids.R"
def SplitGenusList(genus_list):
genus_list = open(genus_list)
contents = genus_list.read()
contents_split = contents.splitlines()
return contents_split
rule extract_microbial_reads:
input:
taxID_list="data/microbial_taxids/{sample}/genus_list.txt",
kraken_file="data/kraken2_classification/{sample}/{sample}.kraken",
reads=config["reads"] "{sample}.fastq"
output:
touch("data/microbial_reads/{sample}/extracting_microbial_reads.done")
run:
import os
genus_list = SplitGenusList("data/microbial_taxids/" wildcards.sample "/genus_list.txt")
for taxID in genus_list:
os.system('python code/extract_kraken_reads.py --kraken ' input.kraken_file ' -s ' input.reads
' -o "data/microbial_reads/' wildcards.sample '/' taxID '_.fastq" '
'-t ' taxID ' --include-children')
轉載請註明出處，本文鏈接：https://www.uj5u.com/ruanti/509834.html
標籤：Python重击蛇制造
上一篇：如何輸出浮點結果？[復制]