我的snakemake管道包含31条规则,这让我发疯了。它是一个映射和snp调用管道,它使用BWA和HaplotypeCaller等。根据所使用的程序,我为每个规则创建了一个conda环境。我的代码相当长,如果需要,可以在以下地址看到:https://github.com/ltalignani/SHAVE1
具体来说,当我想构建DAG时,snakemake告诉我,haplotype_caller规则没有参考基因组作为输入。但它在档案里。以下是相关代码:
rule haplotype_caller_gvcf:
# Aim: Call germline SNPs and indels via local re-assembly of haplotypes
# Use: gatk --java-options '-Xmx{MEM_GB}g' HaplotypeCaller \
# -R Homo_sapiens_assembly38.fasta \
# -I input.bam \
# -O output.g.vcf.gz \
# -ERC GVCF # Essential to GenotypeGVCFs: produce genotype likelihoods
message:
"HaplotypeCaller calling SNVs and Indels for {wildcards.sample} sample ({wildcards.aligner}-{wildcards.mincov})"
conda:
GATK4
input:
refpath = REFPATH,
reference = REFERENCE,
bam = "results/04_Variants/{sample}_{aligner}_{mincov}X_indel-qual.bam"
output:
gvcf="results/04_Variants/haplotypecaller/{sample}_{aligner}_{mincov}X_variant-call.g.vcf"
log:
"results/11_Reports/haplotypecaller/{sample}_{aligner}_{mincov}X_variant-call.log" # optional
resources:
mem_gb= MEM_GB,
shell:
"gatk HaplotypeCaller " # --java-options '-Xmx{resources.mem_gb}g'
"-R {input.refpath}{input.reference} "
"-I {input.bam} "
"-O {output.gvcf} "
"-ERC GVCF" # Essential to GenotypeGVCFs: produce genotype likelihoods在snakefile头中定义REFPATH和引用变量如下:
REFPATH = config["consensus"]["path"] # Path to reference genome REFERENCE = config["consensus"]["reference"] # Genome reference sequence, in fasta format
.yaml中的配置文件如下所示:
consensus:
`reference: "GCA_018104305.1_AalbF3_genomic.fasta"``path: "resources/genomes/" # Path to genome reference`当我要DAG的时候:
snakemake -s workflow/rules/shave.smk --dag | dot -Tpng > test.png我知道这个错误:
`MissingInputException in line 247 of /Users/loic/snakemake/short-read-alignment-vector-pipeline/workflow/rules/shave.smk:`Missing input files for rule haplotype_caller_gvcf:
GCA_018104305.1_AalbF3_genomic.fasta
这是蛇形动物的结构:
也试图使用蛇皮棉,但输出是可以的。
发布于 2022-10-14 16:12:28
我已经查看了您的Github,并且文件夹resources/genomes/只包含一个文件GCA_018104305.1_AalbF3_genomic.fasta.fai。您是否尝试过将该文件重命名为预期的输入名GCA_018104305.1_AalbF3_genomic.fasta,例如去掉.fai扩展名?
发布于 2022-10-16 02:18:51
谢谢你的回答,fasta是1.47GB。这就是为什么它不在资源/基因组文件夹中。.fai是fasta索引,对于一些像GATK这样的程序来说是必需的。
https://stackoverflow.com/questions/74067011
复制相似问题