samples = ["DM8909"]

rule all:
    input:
        "/home/huangqiang/whole_genome_pipline/temp/seqkit/second_sequence_statistics.txt"
      
rule fastp_map:
    input:
        r1 = "seq/{sample}_1.fastq.gz",
        r2 = "seq/{sample}_2.fastq.gz"
    output:
        r1_clean = "temp/fastp_output/{sample}_R1_clean.fastq.gz",
        r2_clean = "temp/fastp_output/{sample}_R2_clean.fastq.gz",
        html_report = "temp/fastp_output/{sample}_fastp.html",
        json_stats = "temp/fastp_output/{sample}_fastp.json"
    shell:
        "fastp \
        -i {input.r1} \
        -I {input.r2} \
        -o {output.r1_clean} \
        -O {output.r2_clean} \
        -h {output.html_report} \
        -j {output.json_stats} \
        -q 15 -u 40 -l 50 --dedup"

rule mkdir_seqkit:
    output:
        "temp/seqkit"
    shell:
        "mkdir -p {output}"

rule second_sequence_statistics:
    input:
        expand("/home/huangqiang/whole_genome_pipline/temp/fastp_output/{sample}_{pair}_clean.fastq.gz", sample=samples, pair=["R2","R1"])
    output:
        "/home/huangqiang/whole_genome_pipline/temp/seqkit/second_sequence_statistics.txt"
    shell:
        "seqkit \
        stats -aT -i {input} > {output}"

the result showed R1 has a result but R2 missed

file    format  type    num_seqs    sum_len min_len avg_len max_len Q1  Q2  Q3  sum_gap N50 N50_num Q20(%)  Q30(%)  AvgQual GC(%)   sum_n
/home/huangqiang/whole_genome_pipline/temp/fastp_output/DM8909_R1_clean.fastq.gz    FASTQ   DNA 4121308 618107404   50  150.0   150 150.0   150.0   150.0   0   150 1   98.24   94.91   27.70   49.67   9786

samples = ["DM8909"]

rule all:
    input:
        "/home/huangqiang/whole_genome_pipline/temp/seqkit/second_sequence_statistics.txt"
      
rule fastp_map:
    input:
        r1 = "seq/{sample}_1.fastq.gz",
        r2 = "seq/{sample}_2.fastq.gz"
    output:
        r1_clean = "temp/fastp_output/{sample}_R1_clean.fastq.gz",
        r2_clean = "temp/fastp_output/{sample}_R2_clean.fastq.gz",
        html_report = "temp/fastp_output/{sample}_fastp.html",
        json_stats = "temp/fastp_output/{sample}_fastp.json"
    shell:
        "fastp \
        -i {input.r1} \
        -I {input.r2} \
        -o {output.r1_clean} \
        -O {output.r2_clean} \
        -h {output.html_report} \
        -j {output.json_stats} \
        -q 15 -u 40 -l 50 --dedup"

rule mkdir_seqkit:
    output:
        "temp/seqkit"
    shell:
        "mkdir -p {output}"

rule second_sequence_statistics:
    input:
        expand("/home/huangqiang/whole_genome_pipline/temp/fastp_output/{sample}_{pair}_clean.fastq.gz", sample=samples, pair=["R2","R1"])
    output:
        "/home/huangqiang/whole_genome_pipline/temp/seqkit/second_sequence_statistics.txt"
    shell:
        "seqkit \
        stats -aT -i {input} > {output}"

the result showed R1 has a result but R2 missed

file    format  type    num_seqs    sum_len min_len avg_len max_len Q1  Q2  Q3  sum_gap N50 N50_num Q20(%)  Q30(%)  AvgQual GC(%)   sum_n
/home/huangqiang/whole_genome_pipline/temp/fastp_output/DM8909_R1_clean.fastq.gz    FASTQ   DNA 4121308 618107404   50  150.0   150 150.0   150.0   150.0   0   150 1   98.24   94.91   27.70   49.67   9786

• snakemake

1 Answer 1

I am not sure this solves your problem, but three things come to mind while looking at your code. A) You should stay consistent with absolute and relative paths (I would say, only use absolute paths if really necessary). B) You do not need to create folders, snakemake does this for you. C) Use triple quotes for multi line code blocks

This will results in:

samples = ["DM8909"]

rule all:
    input:
        "temp/seqkit/second_sequence_statistics.txt"
      
rule fastp_map:
    input:
        r1 = "seq/{sample}_1.fastq.gz",
        r2 = "seq/{sample}_2.fastq.gz"
    output:
        r1_clean = "temp/fastp_output/{sample}_R1_clean.fastq.gz",
        r2_clean = "temp/fastp_output/{sample}_R2_clean.fastq.gz",
        html_report = "temp/fastp_output/{sample}_fastp.html",
        json_stats = "temp/fastp_output/{sample}_fastp.json"
    shell:
        """
        fastp \
        -i {input.r1} \
        -I {input.r2} \
        -o {output.r1_clean} \
        -O {output.r2_clean} \
        -h {output.html_report} \
        -j {output.json_stats} \
        -q 15 -u 40 -l 50 --dedup
        """


rule second_sequence_statistics:
    input:
        expand("temp/fastp_output/{sample}_{pair}_clean.fastq.gz", sample=samples, pair=["R2","R1"])
    output:
        "temp/seqkit/second_sequence_statistics.txt"
    shell:
        """
        seqkit \
        stats -aT -i {input} > {output}
        """

I don't know the syntax for seqkit and fastp, but this will run the following commands:

fastp -i seq/DM8909_1.fastq.gz -I seq/DM8909_2.fastq.gz -o temp/fastp_output/DM8909_R1_clean.fastq.gz -O temp/fastp_output/DM8909_R2_clean.fastq.gz -h temp/fastp_output/DM8909_fastp.html -j temp/fastp_output/DM8909_fastp.json -q 15 -u 40 -l 50 --dedup

seqkit stats -aT -i temp/fastp_output/DM8909_R2_clean.fastq.gz temp/fastp_output/DM8909_R1_clean.fastq.gz > temp/seqkit/second_sequence_statistics.txt