diff --git a/eda_qc.md b/eda_qc.md index 5cff7c2..c184383 100644 --- a/eda_qc.md +++ b/eda_qc.md @@ -118,6 +118,8 @@ median(bcrank$total) ``` Just 2! Clearly many barcodes produce practically no output. + + ::: :::: @@ -339,21 +341,18 @@ sce$discard <- reasons$discard :::: challenge -We've removed empty cells and low-quality cells to be discarded. How many cells are we left with at this point? +Maybe our sample preparation was poor and we want the QC to be more strict. How could we change the set the QC filtering to use 4 MADs as the threshold for outlier calling? ::: solution +You set `nmads = 4` like so: -``` r -table(sce$discard) -``` - -``` output -FALSE TRUE - 2437 694 +``` r +reasons_strict <- perCellQCFilters(df, sub.fields = "subsets_Mito_percent", nmads = 4) ``` -There are 2437 cells that *haven't* been flagged to be discarded, so that's how many we have left. +You would then need to reassign the `discard` column as well, but we'll stick with the 3 MADs default for now. + ::: :::: @@ -543,12 +542,32 @@ altExpNames(0): :::: challenge -Some sophisticated experiments perform additional steps so that they can estimate size factors from so-called "spike-ins". Judging by the name, what do you think "spike-ins" are, and what additional steps are required to use them? +Fill in the blanks for normalization that uses simpler library size factors instead of deconvolution. + + +``` r +____ <- ____SizeFactors(sce) + +sizeFactors(sce) <- ____ + +sce <- ____(sce) + +sce +``` ::: solution -Spike-ins are deliberately-introduced exogeneous RNA from an exotic or synthetic source at a known concentration. This provides a known signal to normalize to. Exotic or synthetic RNA (e.g. soil bacteria RNA in a study of human cells) is used in order to avoid confusing spike-in RNA with sample RNA. This has the obvious advantage of accounting for cell-wise variation, but adds additional sample-preparation work. +``` r +lib.sf <- librarySizeFactors(sce) + +sizeFactors(sce) <- lib.sf + +sce <- logNormCounts(sce) + +sce +``` +If you run this chunk, make sure to go back and re-run the normalization with deconvolution normalization if you want your work to align with the rest of this episode. ::: :::: @@ -582,13 +601,13 @@ ggplot(mean_var_df, aes(mean, var)) + y = "Variance of log-expression") ``` - + The blue line represents the uninteresting "technical" variance for any given gene abundance. The genes with a lot of additional variance exhibit interesting "biological" variation. ### Selecting highly variable genes -The next step is to select the subset of HVGs to use in downstream analyses. A larger set will assure that we do not remove important genes, at the cost of potentially increasing noise. Typically, we restrict ourselves to the top $n$ genes, here we chose $n = 1000$, but this choice should be guided by prior biological knowledge; for instance, we may expect that only about 10% of genes to be differentially expressed across our cell populations and hence select 10% of genes as higly variable (e.g., by setting `prop = 0.1`). +The next step is to select the subset of HVGs to use in downstream analyses. A larger set will assure that we do not remove important genes, at the cost of potentially increasing noise. Typically, we restrict ourselves to the top $n$ genes, here we chose $n = 1000$, but this choice should be guided by prior biological knowledge; for instance, we may expect that only about 10% of genes to be differentially expressed across our cell populations and hence select 10% of genes as highly variable (e.g., by setting `prop = 0.1`). ``` r @@ -604,9 +623,26 @@ head(hvg.sce.var) :::: challenge -Run an internet search for some of the most highly variable genes we've identified here. See if you can identify the type of protein they produce or what sort of process they're involved in. Do they make biological sense to you? +Imagine you have data that were prepared by three people with varying level of experience, which leads to varying technical noise. How can you account for this blocking structure when selecting HVGs? -:::: +::: solution +Use the `block` argument in the call to `modelGeneVar()` like so: + + +``` r +sce$experimenter = factor(sample(c("Perry", "Merry", "Gary"), + replace = TRUE, + size = ncol(sce))) + +blocked_variance_df = modelGeneVar(sce, + block = sce$experimenter) +``` + +Blocked models are evaluated on each block separately then combined. If the experimental groups are related in some structured way, it may be preferable to use the `design` argument. See `?modelGeneVar` for more detail. + +::: + +::: ## Dimensionality Reduction @@ -658,7 +694,7 @@ ggplot(pct_var_df, labs(y = "Variance explained (%)") ``` - + You can see the first two PCs capture the largest amount of variation, but in this case you have to take the first 8 PCs before you've captured 50% of the total. @@ -669,7 +705,7 @@ And we can of course visualize the first 2-3 components, perhaps color-coding ea plotPCA(sce, colour_by = "sum") ``` - + It can be helpful to compare pairs of PCs. This can be done with the `ncomponents` argument to `plotReducedDim()`. For example if one batch or cell type splits off on a particular PC, this can help visualize the effect of that. @@ -678,7 +714,7 @@ It can be helpful to compare pairs of PCs. This can be done with the `ncomponent plotReducedDim(sce, dimred = "PCA", ncomponents = 3) ``` - + ### Non-linear methods @@ -695,7 +731,7 @@ sce <- runTSNE(sce, dimred = "PCA") plotTSNE(sce) ``` - + ``` r @@ -706,7 +742,7 @@ sce <- runUMAP(sce, dimred = "PCA") plotUMAP(sce) ``` - + It is easy to over-interpret t-SNE and UMAP plots. We note that the relative sizes and positions of the visual clusters may be misleading, as they tend to inflate dense clusters and compress sparse ones, such that we cannot use the size as a measure of subpopulation heterogeneity. @@ -793,7 +829,7 @@ sce$DoubletScore <- dbl.dens plotTSNE(sce, colour_by = "DoubletScore") ``` - + We can explicitly convert this into doublet calls by identifying large outliers for the score within each sample. Here we use the "griffiths" method to do so. @@ -816,13 +852,13 @@ sce$doublet <- dbl.calls plotColData(sce, y = "DoubletScore", colour_by = "doublet") ``` - + ``` r plotTSNE(sce, colour_by = "doublet") ``` - + One way to determine whether a cell is in a real transient state or it is a doublet is to check the number of detected genes and total UMI counts. @@ -831,13 +867,13 @@ One way to determine whether a cell is in a real transient state or it is a doub plotColData(sce, "detected", "sum", colour_by = "DoubletScore") ``` - + ``` r plotColData(sce, "detected", "sum", colour_by = "doublet") ``` - + In this case, we only have a few doublets at the periphery of some clusters. It could be fine to keep the doublets in the analysis, but it may be safer to discard them to avoid biases in downstream analysis (e.g., differential expression). @@ -873,6 +909,29 @@ The package `DropletTestFiles` includes the raw output from Cell Ranger of the p :::::::::::::::::::::::::::::::::: + +:::: challenge + +#### Extension challenge 1: Spike-ins + +Some sophisticated experiments perform additional steps so that they can estimate size factors from so-called "spike-ins". Judging by the name, what do you think "spike-ins" are, and what additional steps are required to use them? + +::: solution + +Spike-ins are deliberately-introduced exogeneous RNA from an exotic or synthetic source at a known concentration. This provides a known signal to normalize against. Exotic (e.g. soil bacteria RNA in a study of human cells) or synthetic RNA is used in order to avoid confusing spike-in RNA with sample RNA. This has the obvious advantage of accounting for cell-wise variation, but can substantially increase the amount of sample-preparation work. + +::: + +:::: + +:::: challenge + +#### Extension challenge 2: Background research + +Run an internet search for some of the most highly variable genes we identified in the feature selection section. See if you can identify the type of protein they produce or what sort of process they're involved in. Do they make biological sense to you? + +:::: + ::::::::::::::::::::::::::::::::::::: keypoints - Empty droplets, i.e. droplets that do not contain intact cells and that capture only ambient or background RNA, should be removed prior to an analysis. The `emptyDrops` function from the [DropletUtils](https://bioconductor.org/packages/DropletUtils) package can be used to identify empty droplets. diff --git a/fig/eda_qc-rendered-unnamed-chunk-21-1.png b/fig/eda_qc-rendered-unnamed-chunk-21-1.png new file mode 100644 index 0000000..181c979 Binary files /dev/null and b/fig/eda_qc-rendered-unnamed-chunk-21-1.png differ diff --git a/fig/eda_qc-rendered-unnamed-chunk-25-1.png b/fig/eda_qc-rendered-unnamed-chunk-25-1.png index e3fb09e..a8555d0 100644 Binary files a/fig/eda_qc-rendered-unnamed-chunk-25-1.png and b/fig/eda_qc-rendered-unnamed-chunk-25-1.png differ diff --git a/fig/eda_qc-rendered-unnamed-chunk-26-1.png b/fig/eda_qc-rendered-unnamed-chunk-26-1.png index ff232ae..de7cec1 100644 Binary files a/fig/eda_qc-rendered-unnamed-chunk-26-1.png and b/fig/eda_qc-rendered-unnamed-chunk-26-1.png differ diff --git a/fig/eda_qc-rendered-unnamed-chunk-27-1.png b/fig/eda_qc-rendered-unnamed-chunk-27-1.png new file mode 100644 index 0000000..c803807 Binary files /dev/null and b/fig/eda_qc-rendered-unnamed-chunk-27-1.png differ diff --git a/fig/eda_qc-rendered-unnamed-chunk-28-1.png b/fig/eda_qc-rendered-unnamed-chunk-28-1.png index 433225a..e3fb09e 100644 Binary files a/fig/eda_qc-rendered-unnamed-chunk-28-1.png and b/fig/eda_qc-rendered-unnamed-chunk-28-1.png differ diff --git a/fig/eda_qc-rendered-unnamed-chunk-29-1.png b/fig/eda_qc-rendered-unnamed-chunk-29-1.png index 8cf9120..ff232ae 100644 Binary files a/fig/eda_qc-rendered-unnamed-chunk-29-1.png and b/fig/eda_qc-rendered-unnamed-chunk-29-1.png differ diff --git a/fig/eda_qc-rendered-unnamed-chunk-31-1.png b/fig/eda_qc-rendered-unnamed-chunk-31-1.png new file mode 100644 index 0000000..433225a Binary files /dev/null and b/fig/eda_qc-rendered-unnamed-chunk-31-1.png differ diff --git a/fig/eda_qc-rendered-unnamed-chunk-32-1.png b/fig/eda_qc-rendered-unnamed-chunk-32-1.png new file mode 100644 index 0000000..8cf9120 Binary files /dev/null and b/fig/eda_qc-rendered-unnamed-chunk-32-1.png differ diff --git a/fig/eda_qc-rendered-unnamed-chunk-32-2.png b/fig/eda_qc-rendered-unnamed-chunk-32-2.png new file mode 100644 index 0000000..2ab61c9 Binary files /dev/null and b/fig/eda_qc-rendered-unnamed-chunk-32-2.png differ diff --git a/fig/eda_qc-rendered-unnamed-chunk-33-1.png b/fig/eda_qc-rendered-unnamed-chunk-33-1.png new file mode 100644 index 0000000..f9d0b4d Binary files /dev/null and b/fig/eda_qc-rendered-unnamed-chunk-33-1.png differ diff --git a/fig/eda_qc-rendered-unnamed-chunk-33-2.png b/fig/eda_qc-rendered-unnamed-chunk-33-2.png new file mode 100644 index 0000000..9cc8b17 Binary files /dev/null and b/fig/eda_qc-rendered-unnamed-chunk-33-2.png differ diff --git a/md5sum.txt b/md5sum.txt index b9507b3..5642182 100644 --- a/md5sum.txt +++ b/md5sum.txt @@ -5,7 +5,7 @@ "index.md" "495939ddd3f110be3bbcd49b60f4a7ce" "site/built/index.md" "2024-09-06" "links.md" "8184cf4149eafbf03ce8da8ff0778c14" "site/built/links.md" "2024-09-06" "episodes/intro-sce.Rmd" "e0218183fa139dbe5adbc2ed0b02b8c0" "site/built/intro-sce.md" "2024-09-16" -"episodes/eda_qc.Rmd" "1e88f395d30778f4526532deea43eb03" "site/built/eda_qc.md" "2024-09-06" +"episodes/eda_qc.Rmd" "92b605c153422bca1a40e600f0f021ff" "site/built/eda_qc.md" "2024-09-16" "episodes/cell_type_annotation.Rmd" "66af56b730aaa88e937bc1743afb471a" "site/built/cell_type_annotation.md" "2024-09-08" "episodes/multi-sample.Rmd" "2d38d9903358ea8a8067abd82a1f1f54" "site/built/multi-sample.md" "2024-09-08" "episodes/large_data.Rmd" "b9710492c6792ea435778c4e42f27e02" "site/built/large_data.md" "2024-09-09"