Я пишу конвейер snakemake, чтобы взять общедоступные файлы sra, преобразовать их в файлы fastq, а затем запустить их через выравнивание, пиковый вызов и регрессию оценки LD.
У меня проблема в правиле под названием SRA2fastq
ниже, в котором я использую parallel-fastq-dump
для преобразования файлов SRA в файлы с парным концом fastq. Это правило генерирует два вывода для каждого файла SRA: SRRXXXXXXX_1
и SRRXXXXXXX_2
.
Вот мой файл конфигурации:
samples:
fullard2018_NpfcATAC_1: SRR5367824
fullard2018_NpfcATAC_2: SRR5367798
fullard2018_NpfcATAC_3: SRR5367778
fullard2018_NpfcATAC_4: SRR5367754
fullard2018_NpfcATAC_5: SRR5367729
И вот несколько первых правил моего Snakefile:
# read config info into this namespace
configfile: "config.yaml"
print (config['samples'])
rule all:
input:
expand("fastq_files/{SRA}_{num}.fastq.gz", SRA=[config['samples'][x] for x in config['samples']], num=[1,2]),
expand("FastQC/{SRA}_{num}_fastqc.html", SRA=[config['samples'][x] for x in config['samples']], num=[1,2]),
"FastQC/fastq_multiqc.html",
expand("peak_files/{sample}_peaks.blrm.narrowPeak", sample=config['samples']),
"peak_files/Fullard2018_peaks.mrgd.blrm.narrowPeak",
expand("LD_annotation_files/Fullard_2018.{chr}.l2.ldscore.gz", chr=range(1,23))
rule SRA_prefetch:
params:
SRA="{SRA}"
output:
"/home/c1477909/ncbi/public/sra/{SRA}.sra"
log:
"logs/prefetch/{SRA}.log"
shell:
"prefetch {params.SRA}"
rule SRA2fastq:
input:
"/home/c1477909/ncbi/public/sra/{SRA}.sra"
output:
"fastq_files/{SRA}_1.fastq.gz",
"fastq_files/{SRA}_2.fastq.gz"
log:
"logs/SRA2fastq/{SRA}.log"
shell:
"""
parallel-fastq-dump --sra-id {input} --threads 8 \
--outdir fastq_files --split-files --gzip
"""
rule fastqc:
input:
rules.SRA2fastq.output
output:
# Output needs to end in '_fastqc.html' for multiqc to work
html="FastQC/{SRA}_{num}_fastqc.html"
log:
"logs/FASTQC/{SRA}_{num}.log"
wrapper:
"0.27.1/bio/fastqc"
rule multiqc_fastq:
input:
lambda wildcards: expand("FastQC/{SRA}_{num}_fastqc.html", SRA=[config['samples'][x] for x in config['samples']], num=[1,2])
output:
"FastQC/fastq_multiqc.html"
wrapper:
"0.27.1/bio/multiqc"
rule bowtie2:
input:
sample=lambda wildcards: expand("fastq_files/{SRA}_{num}.fastq.gz", SRA=config['samples'][wildcards.sample], num=[1,2])
output:
"bam_files/{sample}.bam"
log:
"logs/bowtie2/{sample}.txt"
params:
index=config["index"], # prefix of reference genome index (built with bowtie2-build),
extra=""
threads: 8
wrapper:
"0.27.1/bio/bowtie2/align"
Однако, когда я запускаю Snakefile, я получаю следующую ошибку:
Error in job SRA2fastq while creating output files fastq_files/SRR5367754_1.fastq.gz, fastq_files/SRR5367754_2.fastq.gz
Я видел эту ошибку много раз раньше, и обычно она возникает, когда имя выходного файла, сгенерированного программой, не совсем соответствует имени выходного файла, которое вы указываете в соответствующем правиле snakemake. Однако здесь дело обстоит не так, как если бы я запускаю команду snakemake, генерируемую для этого конкретного правила отдельно, файлы создаются, как ожидалось, и имена файлов совпадают. Вот пример одного экземпляра правила, взятого после запуска snakemake -np
:
rule SRA2fastq:
input: /home/c1477909/ncbi/public/sra/SRR5367779.sra
output: fastq_files/SRR5367779_1.fastq.gz, fastq_files/SRR5367779_2.fastq.gz
log: logs/SRA2fastq/SRR5367779.log
jobid: 18
wildcards: SRA=SRR5367779
parallel-fastq-dump --sra-id /home/c1477909/ncbi/public/sra/SRR5367779.sra --threads 8 --outdir fastq_files --split-files --gzip
Обратите внимание, что выходные файлы, сгенерированные командой parallel-fastq-dump
, запущенной отдельно (т.е. без использования snakemake), имеют имена, указанные в правиле SRA2fastq
:
ls fastq_files
SRR5367729_1.fastq.gz SRR5367729_2.fastq.gz
Я немного озадачен этим, поскольку эту ошибку обычно легко исправить, но я не могу понять, в чем проблема. Я попытался изменить раздел вывода SRA2fastq
на:
output:
file1="fastq_files/{SRA}_1.fastq.gz",
file2="fastq_files/{SRA}_2.fastq.gz"
Однако это вызывает ту же ошибку. Я также пробовал просто указать один выходной файл, но позже это повлияет на правило bowtie2
, так как я получаю input files missing
ошибку.
Есть идеи, что здесь происходит? Что-то мне не хватает при поиске нескольких выходных файлов в одном правиле?
Огромное спасибо
parallel-fastq-dump
, а затемecho $?
. Возможно, программное обеспечение возвращает код ошибки, даже если оно работает. - person The Unfun Cat   schedule 21.08.2018