Snakemake 无法将功能符号链接识别为输出并删除部分但不是全部输出

如何解决Snakemake 无法将功能符号链接识别为输出并删除部分但不是全部输出

我和我的同事使用 Snakemake 构建了一个读取映射和变体调用管道。我们使用了两组不同的短读全基因组样本。一组样本由新生成的 WGS 数据组成，另一组由来自 NCBI SRA 存储库的样本组成。所有这些都被映射到相同的参考基因组。为清楚起见，样本名称如下：

sra_ids = [
    "SRR1657028","SRR1657029","SRR1575526","SRR1575545","SRR1575527","SRR1575528","SRR1575543","SRR1575544","SRR1575541","SRR1575542","SRR1575539","SRR1575540","SRR1575532","SRR1575531","SRR1575534","SRR1575533","SRR1575538","SRR1575537","SRR1575536","SRR1575535","SRR1575530","SRR1575529"]

    new_samples = [
    "RANO330-OMH_Pedw1","RANO332-OMH_Pedw2","RANO54-OMH_Pedw3","TSINJ32-OMH_Pdia1","TSINJ38-OMH_Pdia2","TSINJ47-OMH_Pdia3","JEJ01-OMH_Pcan1","JEJ3-11-OMH_Pcan2","MERY3-OMH_Pcan3","ANAL10-OMH_Pper1","TOBI5-1-OMH_Pper2","TOBI5-3-OMH_Pper3","DAR4-11-OMH_Ptat1","DAR4-39-OMH_Ptat2","DAR4-5-OMH_Ptat3","JAM4-16-OMH_Pcor1","JAM4-20-OMH_Pcor2","JAM4-7-OMH_Pcor3","KIBO15-OMH_Pdec1","KIBO36-OMH_Pdec2","KIBO44-OMH_Pdec3","KMTEA7-10-OMH_Pver1","KMTEA7-2-OMH_Pver2","KMTEA7-4-OMH_Pver3","DASI5-08-OMH_Alang1","DASI5-16-OMH_Alang2","DASI5-21-OMH_Alang3"]

在我们最初的规则中，我们下载和处理 SRA 数据。我们现在想要通过符号链接，将所有 fastqs（来自两个集合：sra_ids 和 new_samples）拉到一个目录中，并使用相同的文件名格式，以允许在所有下游规则中以相同方式处理所有样本。以下规则旨在创建所述符号链接（注意 sra1 和 sra2 输入函数指向相对路径，而 new1 和 new2 指向绝对路径）：

rule consolidate_fastqs:
    input:
        sra1 = expand(
            "renamed_fastqs/{sample}_fixed_1.fastq.gz",sample=sra_ids),sra2 = expand(
            "renamed_fastqs/{sample}_fixed_2.fastq.gz",new1 = expand(
            os.path.join(fastq_directory,"{sample}_read1.fastq.gz"),sample=new_samples),new2 = expand(
            os.path.join(fastq_directory,"{sample}_read2.fastq.gz"),sample=new_samples)
    output:
        expand(
            "fastqs_consolidated/{sample}_{read}.fastq.gz",sample=initial_sample_list,read=["read1","read2"])
    params:
        threads = 1,mem = 4,t = very_short
    run:
        for i in input.sra1:
            original = i
            basename = i.split("/")[-1].split("_")[0]
            new_name = "fastqs_consolidated/{}_read1.fastq.gz".format(basename)
            shell(
                "ln -srf {original} {new_name} && touch -h {new_name}")
        for i in input.sra2:
            original = i
            basename = i.split("/")[-1].split("_")[0]
            new_name = "fastqs_consolidated/{}_read2.fastq.gz".format(basename)
            shell(
                "ln -srf {original} {new_name} && touch -h {new_name}")
        for i in input.new1:
            original = i
            basename = i.split("/")[-1].split("_")[0]
            new_name = "fastqs_consolidated/{}_read1.fastq.gz".format(basename)
            shell(
                "ln -sf {original} {new_name} && touch -h {new_name}")
        for i in input.new2:
            original = i
            basename = i.split("/")[-1].split("_")[0]
            new_name = "fastqs_consolidated/{}_read2.fastq.gz".format(basename)
            shell(
                "ln -sf {original} {new_name} && touch -h {new_name}")

我们正在使用 CentOS 3.9 和 Slurm 进行作业管理的集群上运行此管道。然而，管道在上述规则下失败，并出现以下错误：

MissingOutputException in line 241 of /scratch/general/lustre/u6035429/DissAssembly/Snakefile:
Job Missing files after 5 seconds:
fastqs_consolidated/RANO330-OMH_Pedw1_read1.fastq.gz
fastqs_consolidated/RANO330-OMH_Pedw1_read2.fastq.gz
fastqs_consolidated/RANO332-OMH_Pedw2_read1.fastq.gz
fastqs_consolidated/RANO332-OMH_Pedw2_read2.fastq.gz
fastqs_consolidated/RANO54-OMH_Pedw3_read1.fastq.gz
fastqs_consolidated/RANO54-OMH_Pedw3_read2.fastq.gz
fastqs_consolidated/TSINJ32-OMH_Pdia1_read1.fastq.gz
fastqs_consolidated/TSINJ32-OMH_Pdia1_read2.fastq.gz
fastqs_consolidated/TSINJ38-OMH_Pdia2_read1.fastq.gz
fastqs_consolidated/TSINJ38-OMH_Pdia2_read2.fastq.gz
fastqs_consolidated/TSINJ47-OMH_Pdia3_read1.fastq.gz
fastqs_consolidated/TSINJ47-OMH_Pdia3_read2.fastq.gz
fastqs_consolidated/JEJ01-OMH_Pcan1_read1.fastq.gz
fastqs_consolidated/JEJ01-OMH_Pcan1_read2.fastq.gz
fastqs_consolidated/JEJ3-11-OMH_Pcan2_read1.fastq.gz
fastqs_consolidated/JEJ3-11-OMH_Pcan2_read2.fastq.gz
fastqs_consolidated/MERY3-OMH_Pcan3_read1.fastq.gz
fastqs_consolidated/MERY3-OMH_Pcan3_read2.fastq.gz
fastqs_consolidated/ANAL10-OMH_Pper1_read1.fastq.gz
fastqs_consolidated/ANAL10-OMH_Pper1_read2.fastq.gz
fastqs_consolidated/TOBI5-1-OMH_Pper2_read1.fastq.gz
fastqs_consolidated/TOBI5-1-OMH_Pper2_read2.fastq.gz
fastqs_consolidated/TOBI5-3-OMH_Pper3_read1.fastq.gz
fastqs_consolidated/TOBI5-3-OMH_Pper3_read2.fastq.gz
fastqs_consolidated/DAR4-11-OMH_Ptat1_read1.fastq.gz
fastqs_consolidated/DAR4-11-OMH_Ptat1_read2.fastq.gz
fastqs_consolidated/DAR4-39-OMH_Ptat2_read1.fastq.gz
fastqs_consolidated/DAR4-39-OMH_Ptat2_read2.fastq.gz
fastqs_consolidated/DAR4-5-OMH_Ptat3_read1.fastq.gz
fastqs_consolidated/DAR4-5-OMH_Ptat3_read2.fastq.gz
fastqs_consolidated/JAM4-16-OMH_Pcor1_read1.fastq.gz
fastqs_consolidated/JAM4-16-OMH_Pcor1_read2.fastq.gz
fastqs_consolidated/JAM4-20-OMH_Pcor2_read1.fastq.gz
fastqs_consolidated/JAM4-20-OMH_Pcor2_read2.fastq.gz
fastqs_consolidated/JAM4-7-OMH_Pcor3_read1.fastq.gz
fastqs_consolidated/JAM4-7-OMH_Pcor3_read2.fastq.gz
fastqs_consolidated/KIBO15-OMH_Pdec1_read1.fastq.gz
fastqs_consolidated/KIBO15-OMH_Pdec1_read2.fastq.gz
fastqs_consolidated/KIBO36-OMH_Pdec2_read1.fastq.gz
fastqs_consolidated/KIBO36-OMH_Pdec2_read2.fastq.gz
fastqs_consolidated/KIBO44-OMH_Pdec3_read1.fastq.gz
fastqs_consolidated/KIBO44-OMH_Pdec3_read2.fastq.gz
fastqs_consolidated/KMTEA7-10-OMH_Pver1_read1.fastq.gz
fastqs_consolidated/KMTEA7-10-OMH_Pver1_read2.fastq.gz
fastqs_consolidated/KMTEA7-2-OMH_Pver2_read1.fastq.gz
fastqs_consolidated/KMTEA7-2-OMH_Pver2_read2.fastq.gz
fastqs_consolidated/KMTEA7-4-OMH_Pver3_read1.fastq.gz
fastqs_consolidated/KMTEA7-4-OMH_Pver3_read2.fastq.gz
fastqs_consolidated/DASI5-08-OMH_Alang1_read1.fastq.gz
fastqs_consolidated/DASI5-08-OMH_Alang1_read2.fastq.gz
fastqs_consolidated/DASI5-16-OMH_Alang2_read1.fastq.gz
fastqs_consolidated/DASI5-16-OMH_Alang2_read2.fastq.gz
fastqs_consolidated/DASI5-21-OMH_Alang3_read1.fastq.gz
fastqs_consolidated/DASI5-21-OMH_Alang3_read2.fastq.gz
This might be due to filesystem latency. If that is the case,consider to increase the wait time with --latency-wait.

Removing output files of failed job consolidate_fastqs since they might be corrupted:
fastqs_consolidated/SRR1657028_read1.fastq.gz,fastqs_consolidated/SRR1657028_read2.fastq.gz,fastqs_consolidated/SRR1657029_read1.fastq.gz,fastqs_consolidated/SRR1657029_read2.fastq.gz,fastqs_consolidated/SRR1575526_read1.fastq.gz,

其他 SRA 样本依此类推。

我们一直在尝试增加延迟时间，但无济于事。此外，我们实时监控了 fastqs_consolidated 目录，并且规则正在为所有 fastq 文件生成有效的符号链接。更奇怪的是，无论出于何种原因，它都无法识别新样本的输出，而只是删除了 SRA 样本的符号链接。即使我们手动创建符号链接，它也会这样做。尽管假设缺少包含新样本的输出，但这些 fastq 文件的功能符号链接仍保留在输出目录 (fastqs_consolidated) 中。

我们尝试了符号链接命令的各种迭代，包括使用和不使用“-r”标志以及使用和不使用从符号链接到原始文件的相对路径（即，../{original} 和 {original }).

对于可能出现的问题，您有什么建议吗？

解决方法

在我看来，依赖命名约定是不可靠且难以维护的。这就是我要做的...

准备一个将样本 ID 链接到 fastq 文件的样本表（表格）。在此表中，您还可以添加一些有关处理可能需要的样品的信息。使用 pandas（更好）或使用您自己的 csv/tsv 解析器读取此文件。样本表中的信息将指导snakemake如何进行。请注意，示例名称可能与 fastq 文件名无关，并且本地 fastq 文件可以位于任何目录中。我认为这是一个干净而灵活的设置。例如，这是sample_sheet.tsv：

sample_id  sra_id      fastq_r1                                 fastq_r2
sample1    SRR1657028  sra/SRR1657028.fastq.R1.gz               sra/SRR1657028.fastq.R2.gz
sample2    NA          fastq_dir/RANO330-OMH_Pedw1.R1.fastq.gz  fastq_dir/RANO330-OMH_Pedw1.R2.fastq.gz

蛇文件：

import pandas

ss = pandas.read_csv('sample_sheet.tsv',sep= '\t')

rule all:
    input:
        expand('vcf/{sample_id}.vcf',sample_id= ss.sample_id),rule sra_download:
    output:
        'sra/{sra_id}.fastq.R1.gz','sra/{sra_id}.fastq.R2.gz',shell:
        r"""
        fastq-dump --sra-id {wildcards.sra_id} ...
        # or better,query ENA archive
        """

rule align:
    input:
        # Use the sample sheet to link sample_id to fastq files
        fq1= lambda wc: ss[ss.sample_id == wc.sample_id].fastq_r1.iloc[0],fq2= lambda wc: ss[ss.sample_id == wc.sample_id].fastq_r2.iloc[0],ref= 'genome.fa'
    output:
        bam= 'bam/{sample_id}.bam',shell:
        r"""
        bwa mem or whatever {input.ref} {input.fq1} {input.fq2} > {output.bam}
        """

rule call_variants:
    input:
        bam= 'bam/{sample_id}.bam',output:
        vcf= 'vcf/{sample_id}.vcf',shell:
        r"""
        call variants {input.bam} > {output.vcf}
        """

无需重命名文件，也无需依赖任何命名方案。如果文件系统上存在fastq文件，snakemake会进行对齐。否则，它将使用 rule sra_download 来获取文件。我还没有测试过这个，可能会有错误，但希望你能明白。

Snakemake 无法将功能符号链接识别为输出并删除部分但不是全部输出

如何解决Snakemake 无法将功能符号链接识别为输出并删除部分但不是全部输出

解决方法

相关推荐