我正在嘗試構建一個用于分析 RNAseq 資料集的蛇形管道。一些資料集可能是 (1) 單端或雙端，以及 (2) 資料集之間的擱淺程度也可能不同。我希望蛇形管道自動檢測這兩個變數并將其填充到 RSEM 輸入中。
我想出了一種方法，使用來自 RSeQC 的 infer_experiment.py 輸出檔案“strandness.txt”（注意前兩行是空的）。
這是 PairEnd Data Fraction of reads failed to determine: 0.0134 Fraction of reads 由“1 ,1--,2 -,2- ”解釋：0.0049 Fraction of reads 由“1 -,1- , 2 ,2--": 0.9817
現在我想使用 AWK 以便使用以下 AWK 代碼從檔案“strandness.txt”中提取“PairEnd”，該代碼在 bash 中運行良好：
awk 'FNR == 3 {print $3}' ./strandedness.txt
正確輸出“PairEnd”
我將代碼放在這個 AWK 腳本“pairend.awk”中：
#!/bin/awk -f
FNR == 3 {print $3}
如果我單獨呼叫它，它作業正常：
./pairend.awk ./strandedness.txt
現在我想在我的 snakemake 檔案中呼叫腳本，但我不知道如何做到這一點。google 上的大多數示例在 shell 中使用 AWK 腳本，但我實際上想使用該腳本輸出一個變數，該變數將用作我的蛇形管道的輸入。我試過這個，但不能讓它作業：
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 代碼，但沒有成功：
錯誤：snakefile 的 ... 行中的 SyntaxError：無效的語法
我點了這條線，但它指向了paired_end=config["project_path"] "workflow/pairend.awk" config["project_path"] "results/{dataset}/RSEM/{sample}_strandness.txt"
如何在snakemake的輸入部分正確運行我/一個AWK腳本（我對在shell部分運行它不感興趣。任何幫助將不勝感激！
uj5u.com熱心網友回復：
因為您正在尋找字串而不是檔案，所以應該將其放入 params 而不是輸入。您也許可以使用該函式從輸入函式呼叫 awk，但我只在指令shell中見過它。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}
"""
- 在你的 shell 指令中運行 awk
rule RSEM_calculate_expression:
input:
...
paired_end=config["project_path"] "results/{dataset}/RSEM/{sample}_strandedness.txt"
shell:
"""
paired=$(awk 'FNR == 3 {{print $3}}' {input.paired_end})
rsem-calculate-expression --star --num-threads {threads} $paired --alignments {input.bam} {input.index} {wildcards.sample} {log}
"""
我更喜歡第一個選項。如果你決定做一些更復雜的事情并且你不需要額外的 awk 知識來理解你在做什么，那么修改起來會更容易。
轉載請註明出處，本文鏈接：https://www.uj5u.com/ruanti/478104.html