Snakemake "Function did not return str or list of str."-CodePudding

I'm trying to run a RNAseq snakemake pipeline. I'm getting stuck on my input function.

import pandas as pd
import os
import fnmatch
import re
 
# --- Importing Configuration Files --- #
configfile: "/DATA/config/config.yaml"
 
table_cols = ['dataset','sample','species','frr','gtf_version','fa_version']
table_samples = pd.read_table('/DATA/config/samples.tsv', header=0, sep='\t', names=table_cols)
DATASET = table_samples.dataset.values.tolist()
SAMPLE = table_samples['sample'].values.tolist()
SPECIES = table_samples.species.values.tolist()
FRR = table_samples.frr.values.tolist()
GTF_VERSION = table_samples.gtf_version.values.tolist()
FA_VERSION = table_samples.fa_version.values.tolist()
 
print(DATASET,SAMPLE,SPECIES,FRR,GTF_VERSION,FA_VERSION)
 
 
rule all:
        input:
                directory(expand(config["project_path"] "resources/starIndex_{species}_{fa_version}_{gtf_version}",zip, species=SPECIES, fa_version=FA_VERSION, gtf_version=GTF_VERSION)),
                expand(config["project_path"] "results/{dataset}/star_aligned_1pass/{sample}_{species}_Aligned.sortedByCoord.out.bam", zip, dataset=DATASET, sample=SAMPLE, species=SPECIES)
                
wildcard_constraints:
        dataset="|".join([re.escape(x) for x in DATASET]),
        sample="|".join([re.escape(x) for x in SAMPLE]),
        species="|".join([re.escape(x) for x in SPECIES]),
        gtf_version="|".join([re.escape(x) for x in GTF_VERSION]),
        fa_version="|".join([re.escape(x) for x in FA_VERSION])
 
 
## rule starIndex ##  Create star index if it does not exist yet
rule starIndex:
        priority: 1
        input:
                fasta=expand(config["project_path"] "resources/{species}.{fa_version}.dna.primary_assembly.fa",zip, species=SPECIES, fa_version=FA_VERSION),
                gtf=expand(config["project_path"] "resources/{species}.{gtf_version}.gtf",zip, species=SPECIES, gtf_version=GTF_VERSION)
        output:
                directory(config["project_path"] "resources/starIndex_{species}_{fa_version}_{gtf_version}")
        threads:
                20
        params:
                directory(config["project_path"] "resources/starIndex_{species}_{fa_version}_{gtf_version}")
        conda:
                "envs/DTPedia_bulkRNAseq.yaml"
        shell:
                """
                STAR --runThreadN {threads} --runMode genomeGenerate --genomeDir {output} --genomeFastaFiles {input.fasta} --sjdbGTFfile {input.gtf}
                """
 
rule star_1pass_alignment:
        priority: 4
        input:
                read1=config["project_path"] "resources/raw_datasets/{dataset}/{sample}_{species}_RNA-Seq_1.fastq.gz",
                read2=config["project_path"] "resources/raw_datasets/{dataset}/{sample}_{species}_RNA-Seq_2.fastq.gz",
                index=determine_species,
                prefix=config["project_path"] "results/{dataset}/star_aligned_1pass/{sample}_{species}_"
        output:
                bam=config["project_path"] "results/{dataset}/star_aligned_1pass/{sample}_{species}_Aligned.sortedByCoord.out.bam",
                log=config["project_path"] "results/{dataset}/star_aligned_1pass/{sample}_{species}_Log.final.out",
                sj_1pass=config["project_path"] "results/{dataset}/star_aligned_1pass/{sample}_{species}_SJ.out.tab"
        threads:
                12
        conda:
                "envs/DTPedia_bulkRNAseq.yaml"
        shell:
                """
                STAR --runMode alignReads --genomeDir {input.index} --genomeLoad LoadAndKeep --outSAMtype BAM SortedByCoordinate --limitBAMsortRAM 10000000000 --limitGenomeGenerateRAM 20000000000 --readFilesIn {input.read1} {input.read2} --runThreadN {threads} --readFilesCommand gunzip -c --outFileNamePrefix {input.prefix}
                """

This is the error:

['PRJNA493818_GSE120639_SRP162872', 'PRJNA493818_GSE120639_SRP162872', 'PRJNA362883_GSE93946_SRP097621', 'PRJNA362883_GSE93946_SRP097621'] ['SRR7942395_GSM3406786_sAML_Control_1', 'SRR7942395_GSM3406786_sAML_Control_1', 'SRR5195524_GSM2465521_KrasT_45649_NoDox', 'SRR5195524_GSM2465521_KrasT_45649_NoDox'] ['Homo_sapiens', 'Homo_sapiens', 'Mus_musculus', 'Mus_musculus'] [1, 2, 1, 2] ['GRCh38.106', 'GRCh38.106', 'GRCm39.107', 'GRCm39.107'] ['GRCh38', 'GRCh38', 'GRCm39', 'GRCm39']

Building DAG of jobs...
WorkflowError in line 113 of /DATA/workflow/snakefileV21:
Function did not return str or list of str.

I for example tried modifying the line after return without success and the same error output:

# function determine_species_fasta # function for determining fasta input of correct species to rule starIndex
def determine_species(wildcards):
        read1 = config["project_path"] "resources/raw_datasets/{wildcards.dataset}/{wildcards.sample}_{wildcards.species}_RNA-Seq_1.fastq.gz"
        if fnmatch.fnmatch(read1, '*Homo_sapiens*'):
                return "/DATA/resources/starIndex_Homo_sapiens_GRCh38_GRCh38.106"
        elif fnmatch.fnmatch(read1, '*Mus_musculus*'):
                return "/DATA/resources/starIndex_Mus_musculus_GRCm39_GRCm39.107"

Perhaps the wildcards in read1 = config["project_path"] "resources/raw_datasets/{wildcards.dataset}/{wildcards.sample}_{wildcards.species}_RNA-Seq_1.fastq.gz" are not properly filled in? I also tried unpack() without succes https://snakemake.readthedocs.io/en/v6.0.0/snakefiles/rules.html#input-functions-and-unpack.

I hope you could help (:

EDIT 1

I've changed the code to this after suggestions from @SultanOrazbayev. This snakemake pipeline analyzes RNAseq data from mice and humans. This python input function determines which species' starIndex to use and pastes the relevant directory (not file) that was outputted in the rule starIndex.:

# function determine_species_fasta # function for determining fasta input of correct species to rule starIndex
def determine_species(wildcards):
        read1 = config["project_path"] "resources/raw_datasets/{wildcards.dataset}/{wildcards.sample}_{wildcards.species}_RNA-Seq_1.fastq.gz"
        if fnmatch.fnmatch(read1, '*Homo_sapiens*'):
                return "/DATA/resources/starIndex_Homo_sapiens_GRCh38_GRCh38.106"
        elif fnmatch.fnmatch(read1, '*Mus_musculus*'):
                return "/DATA/resources/starIndex_Mus_musculus_GRCm39_GRCm39.107"
        else:
                raise ValueError(f"Wildcards do not match the desired pattern: {wildcards}")

now give this error:

(base) @darwin:/DATA/workflow$ snakemake -s snakefileV21 --use-conda 
['PRJNA493818_GSE120639_SRP162872', 'PRJNA493818_GSE120639_SRP162872', 'PRJNA362883_GSE93946_SRP097621', 'PRJNA362883_GSE93946_SRP097621'] ['SRR7942395_GSM3406786_sAML_Control_1', 'SRR7942395_GSM3406786_sAML_Control_1', 'SRR5195524_GSM2465521_KrasT_45649_NoDox', 'SRR5195524_GSM2465521_KrasT_45649_NoDox'] ['Homo_sapiens', 'Homo_sapiens', 'Mus_musculus', 'Mus_musculus'] [1, 2, 1, 2] ['GRCh38.106', 'GRCh38.106', 'GRCm39.107', 'GRCm39.107'] ['GRCh38', 'GRCh38', 'GRCm39', 'GRCm39']

Building DAG of jobs...
InputFunctionException in line 115 of /DATA/workflow/snakefileV21:
ValueError: Wildcards do not match the desired pattern: PRJNA493818_GSE120639_SRP162872 SRR7942395_GSM3406786_sAML_Control_1 Homo_sapiens
Wildcards:
dataset=PRJNA493818_GSE120639_SRP162872
sample=SRR7942395_GSM3406786_sAML_Control_1
species=Homo_sapiens

So I suspect this line is not properly executed read1 = config["project_path"] "resources/raw_datasets/{wildcards.dataset}/{wildcards.sample}_{wildcards.species}_RNA-Seq_1.fastq.gz" Any suggestions to correct this?

CodePudding user response：

You are not formatting the string you concatenate to the config value

# function determine_species_fasta # function for determining fasta input of correct species to rule starIndex
def determine_species(wildcards):
        read1 = config["project_path"] f"resources/raw_datasets/{wildcards.dataset}/{wildcards.sample}_{wildcards.species}_RNA-Seq_1.fastq.gz"
                                     # ^ was missing this f!
        if fnmatch.fnmatch(read1, '*Homo_sapiens*'):
                return "/DATA/resources/starIndex_Homo_sapiens_GRCh38_GRCh38.106"
        elif fnmatch.fnmatch(read1, '*Mus_musculus*'):
                return "/DATA/resources/starIndex_Mus_musculus_GRCm39_GRCm39.107"
        else:
                raise ValueError(f"Wildcards do not match the desired pattern: {wildcards}")

However, this function could be simplified by directly considering the wildcards.species value:

# function determine_species_fasta # function for determining fasta input of correct species to rule starIndex
def determine_species(wildcards):
    if wildcards.species == 'Homo_spaiens':
        return "/DATA/resources/starIndex_Homo_sapiens_GRCh38_GRCh38.106"
    if wildcards.species == 'Mus_musculus':
        return "/DATA/resources/starIndex_Mus_musculus_GRCm39_GRCm39.107"
    raise ValueError("")

If you have problems with your functions, the easiest way to debug is add in some print statements to ensure variables look how you expect.

CodePudding user response：

A mechanical problem with the current function is that the conditions do not capture all cases:

def test(x):
   if x>10:
      ...
   elif x>5:
      ...
   # will return None for x<=5

I don't know what the function determine_species should return in your case, but perhaps this will work:

def determine_species(wildcards):
   read1 = ...
   if fnmatch.fnmatch(read1, '*Homo_sapiens*'):
      return ...
    elif fnmatch.fnmatch(read1, '*Mus_musculus*'):
      return ...
    else:
      raise ValueError(f"Wildcards do not match the desired pattern: {wildcards}")

This will err on cases that do not conform to your fnmatch condition.