Commit d9ad8cd3 authored by Davis McCarthy's avatar Davis McCarthy
Browse files

Bug fix to remove-conf.Rmd

parent 2ef880b2
......@@ -3,7 +3,7 @@ output: html_document
---
```{r setup, echo=FALSE}
knitr::opts_chunk$set(out.width='90%', fig.align = 'center', eval=TRUE)
knitr::opts_chunk$set(out.width='90%', fig.align = 'center', eval=TRUE, warning=FALSE, message=FALSE)
knitr::opts_knit$set(root.dir = normalizePath(".."))
```
......@@ -459,7 +459,7 @@ Note that (due to what looks like a bug in this version of `sctransform`) we
need to convert the UMI count matrix to a sparse format to apply sctransform.
```{r sctransform-apply}
```{r sctransform-apply, warning=FALSE, message=FALSE}
umi_sparse <- as(counts(umi.qc), "dgCMatrix")
### Genes expressed in at least 5 cells will be kept
sctnorm_data <- sctransform::vst(umi = umi_sparse, min_cells = 1,
......@@ -475,7 +475,7 @@ assay(umi.qc, "sctrans_norm") <- sctnorm_data$y
Let us look at the NB GLM model parameters estimated by sctransform.
```{r sctransform-params-plot}
```{r sctransform-params-plot, warning=FALSE, message=FALSE}
#sce$log10_total_counts
## Matrix of estimated model parameters per gene (theta and regression coefficients)
sctransform::plot_model_pars(sctnorm_data)
......@@ -484,14 +484,14 @@ sctransform::plot_model_pars(sctnorm_data)
We can look at the effect of sctransform's normalization on three particular
genes, ACTB, POU5F1 (aka OCT4) and CD74.
```{r sctransform-genes-plot}
```{r sctransform-genes-plot, warning=FALSE, message=FALSE}
##c('ACTB', 'Rpl10', 'Cd74')
genes_plot <- c("ENSG00000075624", "ENSG00000204531", "ENSG00000019582")
sctransform::plot_model(sctnorm_data, umi_sparse, genes_plot,
plot_residual = TRUE, cell_attr = as.data.frame(colData(umi.qc)))
```
```{r norm-pca-sctransform, fig.cap = "PCA plot of the tung data after sctransform normalisation (Pearson residuals)."}
```{r norm-pca-sctransform, warning=FALSE, message=FALSE, fig.cap = "PCA plot of the tung data after sctransform normalisation (Pearson residuals)."}
reducedDim(umi.qc, "PCA_sctrans_norm") <- reducedDim(
runPCA(umi.qc[endog_genes, ], exprs_values = "sctrans_norm")
)
......@@ -504,7 +504,7 @@ plotReducedDim(
) + ggtitle("PCA plot: sctransform normalization")
```
```{r norm-ours-rle-sctransform, fig.cap = "Cell-wise RLE of the tung data"}
```{r norm-ours-rle-sctransform, warning=FALSE, message=FALSE, fig.cap = "Cell-wise RLE of the tung data"}
plotRLE(
umi.qc[endog_genes, ],
exprs_values = "sctrans_norm",
......
......@@ -7,7 +7,7 @@ knitr::opts_chunk$set(out.width='90%', fig.align = 'center', eval=TRUE)
knitr::opts_knit$set(root.dir = normalizePath(".."))
```
## Dealing with confounders
## Batch effects
### Introduction
......@@ -165,19 +165,19 @@ matrix in the `lm_batch_indi` slot.
umi.qc$cdr <- umi.qc$total_features_by_counts_endogenous / nrow(umi.qc)
## fit a model just accounting for batch by individual
lm_design_batch1 <- model.matrix(~batch + cdr,
data = coldata(umi.qc)[umi.qc$individual == "na19098",])
data = colData(umi.qc)[umi.qc$individual == "na19098",])
fit_indi1 <- lmfit(logcounts(umi.qc)[, umi.qc$individual == "na19098"], lm_design_batch1)
fit_indi1$coefficients[,1] <- 0 ## replace intercept with 0 to preserve reference batch
resids_lm_batch1 <- residuals(fit_indi1, logcounts(umi.qc)[, umi.qc$individual == "na19098"])
lm_design_batch2 <- model.matrix(~batch + cdr,
data = coldata(umi.qc)[umi.qc$individual == "na19101",])
data = colData(umi.qc)[umi.qc$individual == "na19101",])
fit_indi2 <- lmfit(logcounts(umi.qc)[, umi.qc$individual == "na19101"], lm_design_batch2)
fit_indi2$coefficients[,1] <- 0 ## replace intercept with 0 to preserve reference batch
resids_lm_batch2 <- residuals(fit_indi2, logcounts(umi.qc)[, umi.qc$individual == "na19101"])
lm_design_batch3 <- model.matrix(~batch + cdr,
data = coldata(umi.qc)[umi.qc$individual == "na19239",])
data = colData(umi.qc)[umi.qc$individual == "na19239",])
fit_indi3 <- lmfit(logcounts(umi.qc)[, umi.qc$individual == "na19239"], lm_design_batch3)
fit_indi3$coefficients[,1] <- 0 ## replace intercept with 0 to preserve reference batch
resids_lm_batch3 <- residuals(fit_indi3, logcounts(umi.qc)[, umi.qc$individual == "na19239"])
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment