CODEWITHSUNDEEP
snakemake
Patrick Roelli
Patrick Roelli

Reputation: 46

Missing wildcards in S4 snakemake Object in R

I'm running a workflow with a main Snakefile including rules from the rules folder and calling rscripts from those included rules.

Here are a few lines and their specific files:

Snakefile:

samples = pd.read_table("samples.csv", header=0, sep=',', index_col=0)
rule extract:
    input:
        'summary/umi_expression_matrix.tsv'
include: "rules/extract_expression_single.smk"

rules/extract_expression_single.smk:

rule merge_umi:
    input:
        expand('summary/{sample}_umi_expression_matrix.tsv', sample=samples.index)
    output:
        'summary/umi_expression_matrix.tsv'
    script:
        "../scripts/merge_counts_single.R"

scripts/merge_counts_single.R:

samples = read.csv('samples.csv', header=TRUE, stringsAsFactors=FALSE)$samples
read_list = c()
for (i in 1:length(samples)){
    temp_matrix = read.table(snakemake@input[[i]][1], header=T, stringsAsFactors = F)
    cell_barcodes = colnames(temp_matrix)[-1]
    colnames(temp_matrix) = c("GENE",paste(samples[i], cell_barcodes, sep = "_"))
    read_list=c(read_list, list(temp_matrix))
}

# Little function that allows to merge unequal matrices
merge.all <- function(x, y) {
  merge(x, y, all=TRUE, by="GENE")
}

read_counts <- Reduce(merge.all, read_list)
read_counts[is.na(read_counts)] = 0
rownames(read_counts) = read_counts[,1]
read_counts = read_counts[,-1]
write.table(read_counts, file=snakemake@output[[1]], sep='\t')

The "clean" way to do it would be to call [email protected] to attribute sample names to the script. But for some reason snakemake@wildcards is an empty vector. In python:

print(type(snakemake.wildcards))
print(snakemake.wildcards)
print('done')

gives:

<class 'snakemake.io.Wildcards'>

done

which means it's also empty. So right now I have to rely on getting back to the samples.csv file and getting the sample names there. I will also have to double check matching indexes maybe using greps, don't want the samples and the files to get mixed up.

Any idea why this is happening?

Update:

I've tried adding the sample_name as params to see if this would work and it actually does.

rule merge_umi:
    input:
        expand('summary/{sample}_umi_expression_matrix.tsv', sample=samples.index)
    params:
        sample_name = lambda wildcards: samples.index
    output:
        'summary/umi_expression_matrix.tsv'
    script:
        "../scripts/merge_counts_single.R"

I'm gonna use this for now, but my guess is there is still an issue with the scope of wildcards in included rules. Or maybe I'm doing it wrong.

Upvotes: 1

Views: 827

Answers (1)

anilbey
anilbey

Reputation: 1977

The idea of using wildcards is to call a rule for each value in the wildcards.

If you use the expand function in the input of a rule, then your rule will take all of the wildcard values and create a list of strings. Which means, your rule will be invoked just for once (not for each wildcard value). Per default, expand uses the python itertools function product that yields all combinations of the provided wildcard values.

By doing so, you cannot use that wildcard inside your rule any longer. Because when that rule is invoked, it gets all of the wildcard values and convert them into a list that will be given to your R script just for once (not for each wildcard value).

In your case, using wildcards is not suitable, since your merge_count rule will be run only for once (not for each wildcard value).

Upvotes: 1

Related Questions