Merge pull request #3 from snakemake-workflows/debug

Adapt to newest versions of scater, scran, and SingleCellExperiment.

Author: johanneskoester

Snakefile

@@ -39,9 +39,8 @@ rule all:
                covariate=cells.columns[1:]),
         expand("plots/hvg-pca.{covariate}.svg",
                covariate=cells.columns[1:]),
-        expand("plots/hvg-tsne.{covariate}.perp={perplexity}.seed={seed}.svg",
+        expand("plots/hvg-tsne.{covariate}.seed={seed}.svg",
                covariate=cells.columns[1:],
-               perplexity=[5, 10, 20],
                seed=[23213, 789789, 897354])
 
 

envs/eval.yaml

@@ -1,13 +1,14 @@
 channels:
-  - bioconda
   - conda-forge
+  - bioconda
 dependencies:
-  - bioconductor-scater =1.6
-  - bioconductor-scran =1.6
-  - bioconductor-singlecellexperiment =1.0
+  - bioconductor-scater =1.10
+  - bioconductor-scran =1.10
+  - bioconductor-singlecellexperiment =1.4
   - r-rcolorbrewer =1.1
-  - bioconductor-rbgl =1.54
-  - bioconductor-rgraphviz =2.22
+  - bioconductor-rbgl =1.58
+  - bioconductor-rgraphviz =2.26
   - r-gplots =3.0
-  - bioconductor-limma =3.34
+  - bioconductor-limma =3.38
   - r-rtsne =0.13
+  - bioconductor-biomart =2.38

report/hvg-corr-tsne.rst

@@ -3,9 +3,7 @@ where each point represents a cell.
 
 Cells were colored by {{ snakemake.wildcards.covariate }}.
 T-SNE visualizations can be misleading, due to parameter choices and the fact
-that it is non-deterministic. We therefore run t-SNE for different perplexities
-and random seeds.
-This plot shows results for perplexity={{ snakemake.wildcards.perplexity }} and
-seed={{ snakemake.wildcards.seed }}.
-T-SNE results will only be usable for you, if clustering is similar for different seeds.
+that it is non-deterministic. We therefore run t-SNE for different random seeds.
+This plot shows results for seed={{ snakemake.wildcards.seed }}.
+T-SNE results will only be usable for you if clustering is similar for different seeds.
 See `here <https://distill.pub/2016/misread-tsne/>`_ for additional details.

rules/variance.smk

@@ -80,13 +80,13 @@ rule hvg_tsne:
         sce="analysis/normalized.batch-removed.rds",
         var_cor="tables/hvg-correlations.tsv"
     output:
-        report("plots/hvg-tsne.{covariate}.perp={perplexity}.seed={seed}.svg",
+        report("plots/hvg-tsne.{covariate}.seed={seed}.svg",
                    caption="../report/hvg-corr-tsne.rst",
                    category="Dimension Reduction")
     params:
         fdr=config["model"]["fdr"]
     log:
-        "logs/hvg-tsne.{covariate}.perp={perplexity}.seed={seed}.log"
+        "logs/hvg-tsne.{covariate}.seed={seed}.log"
     conda:
         "../envs/eval.yaml"
     wildcard_constraints:

scripts/batch-effect-removal.R

@@ -16,7 +16,6 @@ model.vars <- colData(sce)
 model.vars$G1 <- cycle.assignments$scores$G1
 model.vars$G2M <- cycle.assignments$scores$G2M
 
-
 # setup design matrix
 model <- ~ 1
 for (variable in snakemake@params[["model_variables"]]) {
@@ -29,6 +28,6 @@ colnames(design) <- snakemake@params[["model_variables"]]
 saveRDS(design, file=snakemake@output[["design_matrix"]])
 
 # remove batch effects based on variables
-batch.removed <- removeBatchEffect(exprs(sce), covariates=model.vars[, snakemake@params[["model_variables"]]])
-norm_exprs(sce) <- batch.removed
+batch.removed <- removeBatchEffect(logcounts(sce), covariates=model.vars[, snakemake@params[["model_variables"]]])
+logcounts(sce) <- batch.removed
 saveRDS(sce, file=snakemake@output[["sce"]])

scripts/cell-cycle.R

@@ -26,7 +26,7 @@ if (species == "mouse") {
 pairs <- readRDS(system.file("exdata", markers, package="scran"))
 
 # obtain assignments
-assignments <- cyclone(sce, pairs, gene.names=rowData(sce)$feature_id)
+assignments <- cyclone(sce, pairs, gene.names=rowData(sce)$ensembl_gene_id)
 
 # store assignments
 saveRDS(assignments, file=snakemake@output[[1]])

scripts/filter-cells.R

@@ -13,7 +13,7 @@ sce <- readRDS(snakemake@input[[1]])
 
 # drop cells with too few counts or expressed features
 libsize.drop <- isOutlier(sce$total_counts, nmads=3, type="lower", log=TRUE)
-feature.drop <- isOutlier(sce$total_features, nmads=3, type="lower", log=TRUE)
+feature.drop <- isOutlier(sce$total_features_by_counts, nmads=3, type="lower", log=TRUE)
 
 # drop cells with too high proportion of mito genes or spike-ins expressed
 mito.drop <- isOutlier(sce$pct_counts_Mt, nmads=3, type="higher")

scripts/hvg-correlation.R

@@ -42,7 +42,7 @@ dev.off()
 # choose significantly correlated genes
 chosen <- unique(c(var.cor$gene1[sig.cor], var.cor$gene2[sig.cor]))
 # get normalized expressions
-norm.exprs <- exprs(sce)[chosen,,drop=FALSE]
+norm.exprs <- logcounts(sce)[chosen,,drop=FALSE]
 
 
 # plot heatmap

scripts/hvg-pca.R

@@ -24,6 +24,6 @@ style <- theme(
 
 # plot PCA
 svg(file=snakemake@output[[1]])
-plotPCA(sce, colour_by=covariate,
-        feature_set=chosen, ncomponents=3) + style
+plotPCA(sce[chosen, ], colour_by=covariate,
+        ncomponents=3) + style
 dev.off()

scripts/hvg-tsne.R

@@ -11,7 +11,6 @@ library(scater)
 library(scran)
 
 seed <- as.integer(snakemake@wildcards[["seed"]])
-perplexity <- as.integer(snakemake@wildcards[["perplexity"]])
 fdr <- snakemake@params[["fdr"]]
 covariate <- gsub("-", ".", snakemake@wildcards[["covariate"]])
 
@@ -29,6 +28,5 @@ style <- theme(
 # plot t-SNE
 svg(file=snakemake@output[[1]])
 set.seed(seed)
-plotTSNE(sce, perplexity=perplexity,
-         colour_by=covariate, feature_set=chosen) + style
+plotTSNE(sce[chosen, ], colour_by=covariate) + style
 dev.off()