Page not found
+The page you requested cannot be found (perhaps it was moved or renamed).
+You may want to try searching to find the page's new location, or use +the table of contents to find the page you are looking for.
+diff --git a/docs/01_sesion2_diapositivas.html b/docs/01_sesion2_diapositivas.html new file mode 100644 index 0000000..a765ba0 --- /dev/null +++ b/docs/01_sesion2_diapositivas.html @@ -0,0 +1,604 @@ + + +
+The page you requested cannot be found (perhaps it was moved or renamed).
+You may want to try searching to find the page's new location, or use +the table of contents to find the page you are looking for.
+Instructora: Laura Gómez-Romero
+Este contenido está basado en las diapositivas de Peter Hickey. Ve las diapositivas aquí. Y en el curso de OSCA, lee el material aquí
+library(BiocFileCache)
+bfc <- BiocFileCache()
+raw.path <- bfcrpath(bfc, file.path(
+ "http://cf.10xgenomics.com/samples",
+ "cell-exp/2.1.0/pbmc4k/pbmc4k_raw_gene_bc_matrices.tar.gz"
+))
+untar(raw.path, exdir = file.path(tempdir(), "pbmc4k"))
+
+library(DropletUtils)
+library(Matrix)
+fname <- file.path(tempdir(), "pbmc4k/raw_gene_bc_matrices/GRCh38")
+sce.pbmc <- read10xCounts(fname, col.names = TRUE)
Dataset de células mononucleares de sangre periférica humana (PBMC) de 10X Genomics
+Descripción aquí
+Zheng, G. X. Y. et al. Massively parallel digital transcriptional profiling of single cells. Nat. Commun. 8, 14049 (2017)
+# gene-annotation
+library(scater)
+rownames(sce.pbmc) <- uniquifyFeatureNames(
+ rowData(sce.pbmc)$ID, rowData(sce.pbmc)$Symbol
+)
+library(EnsDb.Hsapiens.v86)
+location <- mapIds(EnsDb.Hsapiens.v86,
+ keys = rowData(sce.pbmc)$ID,
+ column = "SEQNAME", keytype = "GENEID"
+)
+
+# cell-detection
+set.seed(100)
+e.out <- emptyDrops(counts(sce.pbmc))
+sce.pbmc <- sce.pbmc[, which(e.out$FDR <= 0.001)]
# quality-control
+stats <- perCellQCMetrics(sce.pbmc,
+ subsets = list(Mito = which(location == "MT"))
+)
+high.mito <- isOutlier(stats$subsets_Mito_percent,
+ type = "higher"
+)
+sce.pbmc <- sce.pbmc[, !high.mito]
+
+# normalization
+library(scran)
+set.seed(1000)
+clusters <- quickCluster(sce.pbmc)
+sce.pbmc <- computeSumFactors(sce.pbmc, cluster = clusters)
+sce.pbmc <- logNormCounts(sce.pbmc)
# variance modelling
+set.seed(1001)
+dec.pbmc <- modelGeneVarByPoisson(sce.pbmc)
+top.pbmc <- getTopHVGs(dec.pbmc, prop = 0.1)
# dimensionality-reduction
+set.seed(10000)
+sce.pbmc <- denoisePCA(sce.pbmc,
+ subset.row = top.pbmc,
+ technical = dec.pbmc
+)
+
+set.seed(100000)
+sce.pbmc <- runTSNE(sce.pbmc, dimred = "PCA")
+
+set.seed(1000000)
+sce.pbmc <- runUMAP(sce.pbmc, dimred = "PCA")
¿Para qué sirve el método que estamos usando para reducir la dimensionalidad de los datos?
¿Los HGVs están almacenados en nuestro objeto sce.pbmc?
Clustering es un procedimiento no supervisado par definir grupos de células con perfiles de expresión similares
+Su propósito principal es resumir los datos en un formato digerido susceptible a interpretación humana
+Nos permite asignar etiquetas (por ejemplo, tipos celulares) a las células
+Las técnicas de t-SNE/UMAP han comprimido datos altamente multi-dimensionales en dos dimensiones
+Esta compresión inevitablemente ha provocado la perdida de información
+Por lo tanto, agrupamos sobre los PCs y después visualizamos las identidades de los clusters en la gráfica t-SNE/UMAP
++ +
+Un cluster no implica un tipo celular
+Nosotros podemos definir tantos clusters como queramos y podemos utilizar el algoritmo que más nos acomode
+El clustering, como un microscopio, simplemente es una herramienta para explorar los datos
+Preguntar por el mejor clustering es similar a preguntar cuál es la mejor magnificación en un microscopio sin contenido
+El clustering basado en grafos fue popularizado (más no inventado) por su uso en Seurat
+Objetivo: Construir un grafo en el que cada nodo es una célula que está conectada a sus vecinos más cercanos (otras células con perfiles de expresión similares) en el espacio multidimensional
+Ilustremos como funciona para 20 células
++ +
++ +
++ +
++ +
+De una gráfica KNN se puede construir una gráfica SNN. En este tipo de grafo, dos células estarán conectadas por una arista si comparten alguno de sus vecinos más próximos.
++ +
++ +
+Podemos asignar pesos a cada arista del grafo, basándonos en la similaridad de las células involucradas, dándole pesos más altos a células que están más cercanamente relacionadas
+Para ver los distintos esquemas de pesado puedes consultar la documentación de la función makeSNNGraph del paquete bluster. Algunos ejemplos son:
+A partir de una gráfica SNN pesada podemos aplicar algoritmos para identificar comunidades de células
+Una comunidad es un grupo de células que están más conectadas a células en el mismo grupo que lo que están a células de un grupo diferente
+Cada comunidad representa un cluster
++ +
+Después de la construcción del grafo, no se almacena información adicional más alla de las células vecinas. Esto puede producir subclusters artificiales en regiones con muchas células
+library(scran)
+# Construir grafo usando k= 10 vecinos más cercanos en el espacio definido por el PCA
+g <- buildSNNGraph(sce.pbmc, k = 10, use.dimred = "PCA")
+# Identificar comunidades utilizando el método Walktrap
+clust <- igraph::cluster_walktrap(g)$membership
# Visualizar clusters en una gráfica t-SNE
+library(scater)
+sce.pbmc$cluster <- factor(clust)
+plotReducedDim(sce.pbmc, "TSNE", colour_by = "cluster")
¿Qué pasa si utilizas una k más grande o más pequeña?
+library(scran)
+#Construir grafo usando k= 50 vecinos más cercanos en el espacio definido por el PCA
+g50 <- buildSNNGraph(sce.pbmc, k = 50, use.dimred = "PCA")
+# Identificar comunidades utilizando el método Walktrap
+clust50 <- igraph::cluster_walktrap(g50)$membership
# Visualizar clusters en una gráfica t-SNE
+library(scater)
+sce.pbmc$cluster50 <- factor(clust50)
+plotReducedDim(sce.pbmc, "TSNE", colour_by = "cluster50")
En esta implementación:
+El valor de k puede ser toscamente interpretado como el tamaño anticipado de la subpoblación más pequeña
Si una subpoblación tiene menos que (k+1) células entonces el método será forzado a construir aristas entre células de esa subpoblación y células de otras subpoblaciones, incrementando el riesgo de que la subpoblación en cuestión no forme su propio cluster
# Pesos definidos estilo Jaccard seguidos por clustering de Louvain
+# aka 'clustering estilo Seurat'
+g2 <- buildSNNGraph(sce.pbmc, k = 10, use.dimred = "PCA", type = "jaccard")
+clust2 <- igraph::cluster_louvain(g2)$membership
+sce.pbmc$cluster2 <- factor(clust2)
+plotReducedDim(sce.pbmc, "TSNE", colour_by = "cluster2")
Pipelines basados en Seurat:
+Pipelines basados en Scran:
+Para detalles sobre la seleccion de parámetros y comparaciones: visitar esta página.
+library("patchwork")
+
+## Estilo scran vs estilo Seurat
+plotReducedDim(sce.pbmc, "TSNE", colour_by = "cluster") +
+plotReducedDim(sce.pbmc, "TSNE", colour_by = "cluster2")
Distintas métricas de distancia
+g.num <- buildSNNGraph(sce.pbmc, use.dimred = "PCA", type = "number")
+g.jaccard <- buildSNNGraph(sce.pbmc, use.dimred = "PCA", type = "jaccard")
+g.none <- buildKNNGraph(sce.pbmc, use.dimred = "PCA")
Distintos métodos de clustering
+ +Modularidad es una métrica natural para evaluar la separación entre comunidades/clusters
+La modularidad se define como la diferencia (escalada) entre el peso total observado de las aristas entre los nodos en el mismo cluster y el peso total esperado si los pesos fueran distribuidos aleatoriamente entre todos los pares de nodos
+Nosotros calcularemos un score de modularidad para cada cluster usando las tasas en vez de las diferencias, debido a que las tasas no se ven tan fuertemente influenciadas por el tamaño de los clusters
+library(bluster)
+
+# obteniendo la métrica de modularidad
+ratio <- pairwiseModularity(g, clust, as.ratio = TRUE)
+dim(ratio)
## [1] 19 19
+library(pheatmap)
+pheatmap(log2(ratio + 1),
+ cluster_rows = FALSE,
+ cluster_cols = FALSE,
+ color = colorRampPalette(c("white", "blue"))(100)
+)
Un dataset que contiene clusters bien separados debería contener la mayoría del peso total observado en las entradas diagonales, i.e la mayoría de las aristas ocurren entre células del mismo cluster
+Para más detalles sobre evaluación de la separación entre clusters visite esta página
+Clustering por k-means
+Clustering jerárquico
+Una propiedad deseable de cualquier cluster es que éste sea estable a las perturbaciones en los datos de entrada, de esta manera:
+Uno puede hacer un proceso de bootstrap para evaluar la estabilidad de un algoritmo de clustering en un dataset dado y calcular la coasignación.
+La coasignación es la probabilidad de que células elegidas al azar del cluster X y Y sean asignadas al mismo cluster en la réplica del proceso de bootstrap
+myClusterFUN <- function(x) {
+ g <- buildSNNGraph(x, use.dimred = "PCA", type = "jaccard")
+ igraph::cluster_louvain(g)$membership
+}
+
+originals <- myClusterFUN(sce.pbmc)
+set.seed(0010010100)
+coassign <- bootstrapStability(sce.pbmc,
+ FUN = myClusterFUN,
+ clusters = originals
+)
pheatmap(coassign,
+ cluster_row = FALSE, cluster_col = FALSE,
+ color = rev(viridis::magma(100))
+)
Probabilidad alta de coasignación indica que X es estable con respecto a su separación de Y.
+Queremos altas probabilidades de coasignación en la diagonal
+Debes considerar que el bootstraping solo considera el efecto del ruido de muestreo e ignora otros factores que pueden afectar la reproducibilidad (como efectos de batch)
+Además, una pobre separación puede ser altamente estable
+Mejora la resolución al repetir el proceso de feature selection y clustering dentro de un único cluster
Se enfoca en los HGVs y PCs que son los más relevantes para un cluster específico
Veamos como se comporta la expresión de ciertos genes en nuestros clusters
+CD3E, CCR7, CD69, y CD44 son marcadores de células T de memoria.
+g.full <- buildSNNGraph(sce.pbmc, use.dimred = "PCA")
+clust.full <- igraph::cluster_walktrap(g.full)$membership
+sce.pbmc$clust.full <- factor(clust.full)
+plotExpression(sce.pbmc,
+ features = c("CD3E", "CCR7", "CD69", "CD44"),
+ x = "clust.full", colour_by = "clust.full"
+)
De esta gráfica deducimos que las células T de memoria se encuentran en el cluster 10
+Dentro de las células T de memoria, ¿dónde están las subpoblaciones CD4+ y CD8+?
+# Repetimos TODO el proceso de clustering en el subconjunto de células que hemos identificado como células T de memoria (cluster 10).
+memory <- 10
+sce.memory <- sce.pbmc[, clust.full == memory]
+dec.memory <- modelGeneVar(sce.memory)
+sce.memory <- denoisePCA(sce.memory,
+ technical = dec.memory,
+ subset.row = getTopHVGs(dec.memory, prop = 0.1)
+)
+
+g.memory <- buildSNNGraph(sce.memory, use.dimred = "PCA")
+clust.memory <- igraph::cluster_walktrap(g.memory)$membership
+sce.memory$clust.memory <- factor(clust.memory)
Expresión de CD4 es bajo, por lo tanto, su cambio es modesto, pero la interpretación es clara
+Si tipos celulares o estados celulares se extienden sobre las fronteras de los clusters, entonces un subcluster podría representar contaminación de un tipo celular en un cluster separado
+Un cluster no implica un tipo celular
Nosotros podemos definir tantos clusters como queramos y podemos utilizar el algoritmo que más nos acomode
El clustering, como un microscopio, simplemente es una herramienta para explorar los datos
Preguntar por el mejor clustering es similar a preguntar cuál es la mejor magnificación en un microscopio sin contenido
Clustering basado en grafos es rápido y evita tener que hacer suposiciones fuertes sobre la forma de los clusters o la distribución de las células dentro de cada cluster:
+scran::buildSNNGraph()
igraph::cluster_walktrap() o igraph::cluster_louvain()
Modularidad y estabilidad de los clusters son diagnósticos útiles
El proceso de subclustering podría mejorar la resolución dentro de clusters grandes
## [1] "2023-07-31 15:47:25 EDT"
+
+## user system elapsed
+## 620.242 22.904 1043.971
+
+## ─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────
+## setting value
+## version R version 4.3.1 (2023-06-16)
+## os macOS Ventura 13.4.1
+## system aarch64, darwin20
+## ui RStudio
+## language (EN)
+## collate en_US.UTF-8
+## ctype en_US.UTF-8
+## tz America/New_York
+## date 2023-07-31
+## rstudio 2023.06.0+421 Mountain Hydrangea (desktop)
+## pandoc 3.1.1 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/ (via rmarkdown)
+##
+## ─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────
+## package * version date (UTC) lib source
+## AnnotationDbi * 1.62.2 2023-07-02 [1] Bioconductor
+## AnnotationFilter * 1.24.0 2023-05-08 [1] Bioconductor
+## AnnotationHub * 3.8.0 2023-05-08 [1] Bioconductor
+## beachmat 2.16.0 2023-05-08 [1] Bioconductor
+## beeswarm 0.4.0 2021-06-01 [1] CRAN (R 4.3.0)
+## Biobase * 2.60.0 2023-05-08 [1] Bioconductor
+## BiocFileCache * 2.8.0 2023-05-08 [1] Bioconductor
+## BiocGenerics * 0.46.0 2023-06-04 [1] Bioconductor
+## BiocIO 1.10.0 2023-05-08 [1] Bioconductor
+## BiocManager 1.30.21.1 2023-07-18 [1] CRAN (R 4.3.0)
+## BiocNeighbors 1.18.0 2023-05-08 [1] Bioconductor
+## BiocParallel 1.34.2 2023-05-28 [1] Bioconductor
+## BiocSingular 1.16.0 2023-05-08 [1] Bioconductor
+## BiocVersion 3.17.1 2022-12-20 [1] Bioconductor
+## biomaRt 2.56.1 2023-06-11 [1] Bioconductor
+## Biostrings 2.68.1 2023-05-21 [1] Bioconductor
+## bit 4.0.5 2022-11-15 [1] CRAN (R 4.3.0)
+## bit64 4.0.5 2020-08-30 [1] CRAN (R 4.3.0)
+## bitops 1.0-7 2021-04-24 [1] CRAN (R 4.3.0)
+## blob 1.2.4 2023-03-17 [1] CRAN (R 4.3.0)
+## bluster * 1.10.0 2023-05-08 [1] Bioconductor
+## bookdown 0.34 2023-05-09 [1] CRAN (R 4.3.0)
+## bslib 0.5.0 2023-06-09 [1] CRAN (R 4.3.0)
+## cachem 1.0.8 2023-05-01 [1] CRAN (R 4.3.0)
+## cli 3.6.1 2023-03-23 [1] CRAN (R 4.3.0)
+## cluster 2.1.4 2022-08-22 [1] CRAN (R 4.3.1)
+## codetools 0.2-19 2023-02-01 [1] CRAN (R 4.3.1)
+## colorspace 2.1-0 2023-01-23 [1] CRAN (R 4.3.0)
+## cowplot 1.1.1 2020-12-30 [1] CRAN (R 4.3.0)
+## crayon 1.5.2 2022-09-29 [1] CRAN (R 4.3.0)
+## curl 5.0.1 2023-06-07 [1] CRAN (R 4.3.0)
+## DBI 1.1.3 2022-06-18 [1] CRAN (R 4.3.0)
+## dbplyr * 2.3.3 2023-07-07 [1] CRAN (R 4.3.0)
+## DelayedArray 0.26.6 2023-07-02 [1] Bioconductor
+## DelayedMatrixStats 1.22.0 2023-05-08 [1] Bioconductor
+## digest 0.6.33 2023-07-07 [1] CRAN (R 4.3.0)
+## dplyr * 1.1.2 2023-04-20 [1] CRAN (R 4.3.0)
+## dqrng 0.3.0 2021-05-01 [1] CRAN (R 4.3.0)
+## DropletUtils * 1.20.0 2023-05-08 [1] Bioconductor
+## edgeR 3.42.4 2023-06-04 [1] Bioconductor
+## ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.3.0)
+## EnsDb.Hsapiens.v86 * 2.99.0 2023-07-29 [1] Bioconductor
+## ensembldb * 2.24.0 2023-05-08 [1] Bioconductor
+## evaluate 0.21 2023-05-05 [1] CRAN (R 4.3.0)
+## ExperimentHub 2.8.1 2023-07-16 [1] Bioconductor
+## fansi 1.0.4 2023-01-22 [1] CRAN (R 4.3.0)
+## farver 2.1.1 2022-07-06 [1] CRAN (R 4.3.0)
+## fastmap 1.1.1 2023-02-24 [1] CRAN (R 4.3.0)
+## filelock 1.0.2 2018-10-05 [1] CRAN (R 4.3.0)
+## FNN 1.1.3.2 2023-03-20 [1] CRAN (R 4.3.0)
+## generics 0.1.3 2022-07-05 [1] CRAN (R 4.3.0)
+## GenomeInfoDb * 1.36.1 2023-07-02 [1] Bioconductor
+## GenomeInfoDbData 1.2.10 2023-06-08 [1] Bioconductor
+## GenomicAlignments 1.36.0 2023-05-08 [1] Bioconductor
+## GenomicFeatures * 1.52.1 2023-07-02 [1] Bioconductor
+## GenomicRanges * 1.52.0 2023-05-08 [1] Bioconductor
+## ggbeeswarm 0.7.2 2023-04-29 [1] CRAN (R 4.3.0)
+## ggplot2 * 3.4.2 2023-04-03 [1] CRAN (R 4.3.0)
+## ggrepel * 0.9.3 2023-02-03 [1] CRAN (R 4.3.0)
+## glue 1.6.2 2022-02-24 [1] CRAN (R 4.3.0)
+## gridExtra 2.3 2017-09-09 [1] CRAN (R 4.3.0)
+## gtable 0.3.3 2023-03-21 [1] CRAN (R 4.3.0)
+## HDF5Array 1.28.1 2023-05-08 [1] Bioconductor
+## here 1.0.1 2020-12-13 [1] CRAN (R 4.3.0)
+## highr 0.10 2022-12-22 [1] CRAN (R 4.3.0)
+## hms 1.1.3 2023-03-21 [1] CRAN (R 4.3.0)
+## htmltools 0.5.5 2023-03-23 [1] CRAN (R 4.3.0)
+## httpuv 1.6.11 2023-05-11 [1] CRAN (R 4.3.0)
+## httr 1.4.6 2023-05-08 [1] CRAN (R 4.3.0)
+## igraph 1.5.0.1 2023-07-23 [1] CRAN (R 4.3.0)
+## interactiveDisplayBase 1.38.0 2023-05-08 [1] Bioconductor
+## IRanges * 2.34.1 2023-07-02 [1] Bioconductor
+## irlba 2.3.5.1 2022-10-03 [1] CRAN (R 4.3.0)
+## jquerylib 0.1.4 2021-04-26 [1] CRAN (R 4.3.0)
+## jsonlite 1.8.7 2023-06-29 [1] CRAN (R 4.3.0)
+## KEGGREST 1.40.0 2023-05-08 [1] Bioconductor
+## knitr 1.43 2023-05-25 [1] CRAN (R 4.3.0)
+## labeling 0.4.2 2020-10-20 [1] CRAN (R 4.3.0)
+## later 1.3.1 2023-05-02 [1] CRAN (R 4.3.0)
+## lattice 0.21-8 2023-04-05 [1] CRAN (R 4.3.1)
+## lazyeval 0.2.2 2019-03-15 [1] CRAN (R 4.3.0)
+## lifecycle 1.0.3 2022-10-07 [1] CRAN (R 4.3.0)
+## limma 3.56.2 2023-06-04 [1] Bioconductor
+## locfit 1.5-9.8 2023-06-11 [1] CRAN (R 4.3.0)
+## magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.3.0)
+## Matrix * 1.6-0 2023-07-08 [1] CRAN (R 4.3.0)
+## MatrixGenerics * 1.12.2 2023-07-02 [1] Bioconductor
+## matrixStats * 1.0.0 2023-06-02 [1] CRAN (R 4.3.0)
+## memoise 2.0.1 2021-11-26 [1] CRAN (R 4.3.0)
+## metapod 1.8.0 2023-04-25 [1] Bioconductor
+## mime 0.12 2021-09-28 [1] CRAN (R 4.3.0)
+## munsell 0.5.0 2018-06-12 [1] CRAN (R 4.3.0)
+## patchwork * 1.1.2 2022-08-19 [1] CRAN (R 4.3.0)
+## PCAtools * 2.12.0 2023-05-08 [1] Bioconductor
+## pheatmap * 1.0.12 2019-01-04 [1] CRAN (R 4.3.0)
+## pillar 1.9.0 2023-03-22 [1] CRAN (R 4.3.0)
+## pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.3.0)
+## plyr 1.8.8 2022-11-11 [1] CRAN (R 4.3.0)
+## png 0.1-8 2022-11-29 [1] CRAN (R 4.3.0)
+## prettyunits 1.1.1 2020-01-24 [1] CRAN (R 4.3.0)
+## progress 1.2.2 2019-05-16 [1] CRAN (R 4.3.0)
+## promises 1.2.0.1 2021-02-11 [1] CRAN (R 4.3.0)
+## ProtGenerics 1.32.0 2023-05-08 [1] Bioconductor
+## purrr 1.0.1 2023-01-10 [1] CRAN (R 4.3.0)
+## R.methodsS3 1.8.2 2022-06-13 [1] CRAN (R 4.3.0)
+## R.oo 1.25.0 2022-06-12 [1] CRAN (R 4.3.0)
+## R.utils 2.12.2 2022-11-11 [1] CRAN (R 4.3.0)
+## R6 2.5.1 2021-08-19 [1] CRAN (R 4.3.0)
+## rappdirs 0.3.3 2021-01-31 [1] CRAN (R 4.3.0)
+## RColorBrewer 1.1-3 2022-04-03 [1] CRAN (R 4.3.0)
+## Rcpp 1.0.11 2023-07-06 [1] CRAN (R 4.3.0)
+## RCurl 1.98-1.12 2023-03-27 [1] CRAN (R 4.3.0)
+## reshape2 1.4.4 2020-04-09 [1] CRAN (R 4.3.0)
+## restfulr 0.0.15 2022-06-16 [1] CRAN (R 4.3.0)
+## rhdf5 2.44.0 2023-05-08 [1] Bioconductor
+## rhdf5filters 1.12.1 2023-05-08 [1] Bioconductor
+## Rhdf5lib 1.22.0 2023-05-08 [1] Bioconductor
+## rjson 0.2.21 2022-01-09 [1] CRAN (R 4.3.0)
+## rlang 1.1.1 2023-04-28 [1] CRAN (R 4.3.0)
+## rmarkdown 2.23 2023-07-01 [1] CRAN (R 4.3.0)
+## rprojroot 2.0.3 2022-04-02 [1] CRAN (R 4.3.0)
+## Rsamtools 2.16.0 2023-06-04 [1] Bioconductor
+## RSQLite 2.3.1 2023-04-03 [1] CRAN (R 4.3.0)
+## rstudioapi 0.15.0 2023-07-07 [1] CRAN (R 4.3.0)
+## rsvd 1.0.5 2021-04-16 [1] CRAN (R 4.3.0)
+## rtracklayer 1.60.0 2023-05-08 [1] Bioconductor
+## Rtsne 0.16 2022-04-17 [1] CRAN (R 4.3.0)
+## S4Arrays 1.0.4 2023-05-14 [1] Bioconductor
+## S4Vectors * 0.38.1 2023-05-08 [1] Bioconductor
+## sass 0.4.7 2023-07-15 [1] CRAN (R 4.3.0)
+## ScaledMatrix 1.8.1 2023-05-08 [1] Bioconductor
+## scales 1.2.1 2022-08-20 [1] CRAN (R 4.3.0)
+## scater * 1.28.0 2023-04-25 [1] Bioconductor
+## scran * 1.28.2 2023-07-23 [1] Bioconductor
+## scRNAseq * 2.14.0 2023-04-27 [1] Bioconductor
+## scuttle * 1.9.4 2023-01-23 [1] Bioconductor
+## sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.3.0)
+## shiny 1.7.4.1 2023-07-06 [1] CRAN (R 4.3.0)
+## SingleCellExperiment * 1.22.0 2023-05-08 [1] Bioconductor
+## sparseMatrixStats 1.12.2 2023-07-02 [1] Bioconductor
+## statmod 1.5.0 2023-01-06 [1] CRAN (R 4.3.0)
+## stringi 1.7.12 2023-01-11 [1] CRAN (R 4.3.0)
+## stringr 1.5.0 2022-12-02 [1] CRAN (R 4.3.0)
+## SummarizedExperiment * 1.30.2 2023-06-06 [1] Bioconductor
+## tibble 3.2.1 2023-03-20 [1] CRAN (R 4.3.0)
+## tidyselect 1.2.0 2022-10-10 [1] CRAN (R 4.3.0)
+## utf8 1.2.3 2023-01-31 [1] CRAN (R 4.3.0)
+## uwot 0.1.16 2023-06-29 [1] CRAN (R 4.3.0)
+## vctrs 0.6.3 2023-06-14 [1] CRAN (R 4.3.0)
+## vipor 0.4.5 2017-03-22 [1] CRAN (R 4.3.0)
+## viridis 0.6.4 2023-07-22 [1] CRAN (R 4.3.0)
+## viridisLite 0.4.2 2023-05-02 [1] CRAN (R 4.3.0)
+## withr 2.5.0 2022-03-03 [1] CRAN (R 4.3.0)
+## xfun 0.39 2023-04-20 [1] CRAN (R 4.3.0)
+## XML 3.99-0.14 2023-03-19 [1] CRAN (R 4.3.0)
+## xml2 1.3.5 2023-07-06 [1] CRAN (R 4.3.0)
+## xtable 1.8-4 2019-04-21 [1] CRAN (R 4.3.0)
+## XVector 0.40.0 2023-05-08 [1] Bioconductor
+## yaml 2.3.7 2023-01-23 [1] CRAN (R 4.3.0)
+## zlibbioc 1.46.0 2023-05-08 [1] Bioconductor
+##
+## [1] /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library
+##
+## ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
+
+Dra. Evelia Coss
+07 de agosto de 2023
+Diapositivas de Peter Hickey: Ve las diapositivas aqui.
+Curso 2021 impartido por Leonardo Collado Torres: [Ver informacion] aqui(https://comunidadbioinfo.github.io/cdsb2021_scRNAseq/control-de-calidad.html)
+ +Detectar células verdaderas y de alta calidad ✅🍪, de modo que cuando agrupemos nuestras células sea más fácil identificar poblaciones de tipos celulares distintos.
+Identificar las muestras fallidas e intentar salvar los datos o eliminarlas del análisis, además de intentar comprender por qué falló la muestra.
+Se responden durante el análisis de SingleCell RNA-Seq.
+Cada droplets debe contener una sola célula 🍪.
++Figura obtenida de Single-Cell best practices.
+Una o varias células podemos encontrar:
+Mas informacion sobre ERCC spike-in control transcripts.
+addPerCellQC
La función addPerCellQC
provenie del paquete scater
. Agrega las siguientes metricas en cada célula y por cada gen dentro del mismo archivo.
sum
: Número de cuentas (lecturas) totales de cada célula.detected
: Genes expresados con al menos una cuenta.altexps_ERCC_percent
: Porcentaje de cuentas mapeadas de las secuencias control (ERCC spike-in control transcripts).subsets_Mito_percent
: Porcentaje de cuentas mapeadas provenientes de la mitocondria.Realizaremos las siguientes actividades con este ejemplo:
++Línea celular de células mieloides progenitoras inmortalizadas de ratón usando SmartSeq2.
+Los paquetes que vamos a emplear para esta sección son:
+library(scRNAseq) ## para descargar datos de ejemplo
+library(DropletUtils) ## para detectar droplets
+library(Matrix) ## para leer datos en formatos comprimidos
+library(AnnotationHub) ## para obtener información de genes
+library(scater) ## para gráficas y control de calidad
+library(BiocFileCache) ## para descargar datos
+library(EnsDb.Hsapiens.v86) ## Archivo de anotacion en humanos en Ensembl
+library(dplyr) ## Modificacion de archivos dataframe
Cargar los datos empleando el paquete scRNAseq
## snapshotDate(): 2023-04-24
+## see ?scRNAseq and browseVignettes('scRNAseq') for documentation
+## loading from cache
+## see ?scRNAseq and browseVignettes('scRNAseq') for documentation
+## loading from cache
+## see ?scRNAseq and browseVignettes('scRNAseq') for documentation
+## loading from cache
+## snapshotDate(): 2023-04-24
+## see ?scRNAseq and browseVignettes('scRNAseq') for documentation
+## loading from cache
+## snapshotDate(): 2023-04-24
+## loading from cache
+
+El paquete scRNAseq
contiene multiples data sets compilados en funciones, para saber mas da click aqui.
## snapshotDate(): 2023-04-24
+
+## AnnotationHub with 1 record
+## # snapshotDate(): 2023-04-24
+## # names(): AH73905
+## # $dataprovider: Ensembl
+## # $species: Mus musculus
+## # $rdataclass: EnsDb
+## # $rdatadateadded: 2019-05-02
+## # $title: Ensembl 97 EnsDb for Mus musculus
+## # $description: Gene and protein annotations for Mus musculus based on Ensembl version 97.
+## # $taxonomyid: 10090
+## # $genome: GRCm38
+## # $sourcetype: ensembl
+## # $sourceurl: http://www.ensembl.org
+## # $sourcesize: NA
+## # $tags: c("97", "AHEnsDbs", "Annotation", "EnsDb", "Ensembl", "Gene", "Protein", "Transcript")
+## # retrieve record with 'object[["AH73905"]]'
+# Anotacion de genes con la localizacion de cada cromosoma
+ens.mm.v97 <- ah[["AH73905"]] # solo un cromosoma
## loading from cache
+
+## Warning: Unable to map 563 of 46604 requested IDs.
+
+addPerCellQC
Al imprimir la variable sce.416b
podemos observar que es de tipo SingleCellExperiment.
# Agregar la informacion de la calidad por cada gen y celula en un mismo archivo
+# Identificar genes mitocondriales
+sce.416b <- addPerCellQC(sce.416b,
+ subsets = list(Mito = is.mito)
+)
+sce.416b
## class: SingleCellExperiment
+## dim: 46604 192
+## metadata(0):
+## assays(1): counts
+## rownames(46604): ENSMUSG00000102693 ENSMUSG00000064842 ... ENSMUSG00000095742 CBFB-MYH11-mcherry
+## rowData names(1): Length
+## colnames(192): SLX-9555.N701_S502.C89V9ANXX.s_1.r_1 SLX-9555.N701_S503.C89V9ANXX.s_1.r_1 ...
+## SLX-11312.N712_S508.H5H5YBBXX.s_8.r_1 SLX-11312.N712_S517.H5H5YBBXX.s_8.r_1
+## colData names(21): Source Name cell line ... altexps_SIRV_percent total
+## reducedDimNames(0):
+## mainExpName: endogenous
+## altExpNames(2): ERCC SIRV
+Podemos observar la informacion de este data frame usando colData
:
## DataFrame with 192 rows and 21 columns
+## Source Name cell line cell type single cell well quality
+## <character> <character> <character> <character>
+## SLX-9555.N701_S502.C89V9ANXX.s_1.r_1 SLX-9555.N701_S502.C.. 416B embryonic stem cell OK
+## SLX-9555.N701_S503.C89V9ANXX.s_1.r_1 SLX-9555.N701_S503.C.. 416B embryonic stem cell OK
+## SLX-9555.N701_S504.C89V9ANXX.s_1.r_1 SLX-9555.N701_S504.C.. 416B embryonic stem cell OK
+## SLX-9555.N701_S505.C89V9ANXX.s_1.r_1 SLX-9555.N701_S505.C.. 416B embryonic stem cell OK
+## SLX-9555.N701_S506.C89V9ANXX.s_1.r_1 SLX-9555.N701_S506.C.. 416B embryonic stem cell OK
+## ... ... ... ... ...
+## SLX-11312.N712_S505.H5H5YBBXX.s_8.r_1 SLX-11312.N712_S505... 416B embryonic stem cell OK
+## SLX-11312.N712_S506.H5H5YBBXX.s_8.r_1 SLX-11312.N712_S506... 416B embryonic stem cell OK
+## SLX-11312.N712_S507.H5H5YBBXX.s_8.r_1 SLX-11312.N712_S507... 416B embryonic stem cell OK
+## SLX-11312.N712_S508.H5H5YBBXX.s_8.r_1 SLX-11312.N712_S508... 416B embryonic stem cell OK
+## SLX-11312.N712_S517.H5H5YBBXX.s_8.r_1 SLX-11312.N712_S517... 416B embryonic stem cell OK
+## genotype phenotype strain spike-in addition
+## <character> <character> <character> <character>
+## SLX-9555.N701_S502.C89V9ANXX.s_1.r_1 Doxycycline-inducibl.. wild type phenotype B6D2F1-J ERCC+SIRV
+## SLX-9555.N701_S503.C89V9ANXX.s_1.r_1 Doxycycline-inducibl.. wild type phenotype B6D2F1-J ERCC+SIRV
+## SLX-9555.N701_S504.C89V9ANXX.s_1.r_1 Doxycycline-inducibl.. wild type phenotype B6D2F1-J ERCC+SIRV
+## SLX-9555.N701_S505.C89V9ANXX.s_1.r_1 Doxycycline-inducibl.. induced CBFB-MYH11 o.. B6D2F1-J ERCC+SIRV
+## SLX-9555.N701_S506.C89V9ANXX.s_1.r_1 Doxycycline-inducibl.. induced CBFB-MYH11 o.. B6D2F1-J ERCC+SIRV
+## ... ... ... ... ...
+## SLX-11312.N712_S505.H5H5YBBXX.s_8.r_1 Doxycycline-inducibl.. induced CBFB-MYH11 o.. B6D2F1-J Premixed
+## SLX-11312.N712_S506.H5H5YBBXX.s_8.r_1 Doxycycline-inducibl.. induced CBFB-MYH11 o.. B6D2F1-J Premixed
+## SLX-11312.N712_S507.H5H5YBBXX.s_8.r_1 Doxycycline-inducibl.. induced CBFB-MYH11 o.. B6D2F1-J Premixed
+## SLX-11312.N712_S508.H5H5YBBXX.s_8.r_1 Doxycycline-inducibl.. induced CBFB-MYH11 o.. B6D2F1-J Premixed
+## SLX-11312.N712_S517.H5H5YBBXX.s_8.r_1 Doxycycline-inducibl.. wild type phenotype B6D2F1-J Premixed
+## block sum detected subsets_Mito_sum subsets_Mito_detected
+## <factor> <numeric> <numeric> <numeric> <numeric>
+## SLX-9555.N701_S502.C89V9ANXX.s_1.r_1 20160113 865936 7618 78790 20
+## SLX-9555.N701_S503.C89V9ANXX.s_1.r_1 20160113 1076277 7521 98613 20
+## SLX-9555.N701_S504.C89V9ANXX.s_1.r_1 20160113 1180138 8306 100341 19
+## SLX-9555.N701_S505.C89V9ANXX.s_1.r_1 20160113 1342593 8143 104882 20
+## SLX-9555.N701_S506.C89V9ANXX.s_1.r_1 20160113 1668311 7154 129559 22
+## ... ... ... ... ... ...
+## SLX-11312.N712_S505.H5H5YBBXX.s_8.r_1 20160325 776622 8174 48126 20
+## SLX-11312.N712_S506.H5H5YBBXX.s_8.r_1 20160325 1299950 8956 112225 25
+## SLX-11312.N712_S507.H5H5YBBXX.s_8.r_1 20160325 1800696 9530 135693 23
+## SLX-11312.N712_S508.H5H5YBBXX.s_8.r_1 20160325 46731 6649 3505 16
+## SLX-11312.N712_S517.H5H5YBBXX.s_8.r_1 20160325 1866692 10964 150375 29
+## subsets_Mito_percent altexps_ERCC_sum altexps_ERCC_detected altexps_ERCC_percent
+## <numeric> <numeric> <numeric> <numeric>
+## SLX-9555.N701_S502.C89V9ANXX.s_1.r_1 9.09882 65278 39 6.80658
+## SLX-9555.N701_S503.C89V9ANXX.s_1.r_1 9.16242 74748 40 6.28030
+## SLX-9555.N701_S504.C89V9ANXX.s_1.r_1 8.50248 60878 42 4.78949
+## SLX-9555.N701_S505.C89V9ANXX.s_1.r_1 7.81190 60073 42 4.18567
+## SLX-9555.N701_S506.C89V9ANXX.s_1.r_1 7.76588 136810 44 7.28887
+## ... ... ... ... ...
+## SLX-11312.N712_S505.H5H5YBBXX.s_8.r_1 6.19684 61575 39 7.17620
+## SLX-11312.N712_S506.H5H5YBBXX.s_8.r_1 8.63302 94982 41 6.65764
+## SLX-11312.N712_S507.H5H5YBBXX.s_8.r_1 7.53559 113707 40 5.81467
+## SLX-11312.N712_S508.H5H5YBBXX.s_8.r_1 7.50037 7580 44 13.48898
+## SLX-11312.N712_S517.H5H5YBBXX.s_8.r_1 8.05569 48664 39 2.51930
+## altexps_SIRV_sum altexps_SIRV_detected altexps_SIRV_percent total
+## <numeric> <numeric> <numeric> <numeric>
+## SLX-9555.N701_S502.C89V9ANXX.s_1.r_1 27828 7 2.90165 959042
+## SLX-9555.N701_S503.C89V9ANXX.s_1.r_1 39173 7 3.29130 1190198
+## SLX-9555.N701_S504.C89V9ANXX.s_1.r_1 30058 7 2.36477 1271074
+## SLX-9555.N701_S505.C89V9ANXX.s_1.r_1 32542 7 2.26741 1435208
+## SLX-9555.N701_S506.C89V9ANXX.s_1.r_1 71850 7 3.82798 1876971
+## ... ... ... ... ...
+## SLX-11312.N712_S505.H5H5YBBXX.s_8.r_1 19848 7 2.313165 858045
+## SLX-11312.N712_S506.H5H5YBBXX.s_8.r_1 31729 7 2.224004 1426661
+## SLX-11312.N712_S507.H5H5YBBXX.s_8.r_1 41116 7 2.102562 1955519
+## SLX-11312.N712_S508.H5H5YBBXX.s_8.r_1 1883 7 3.350892 56194
+## SLX-11312.N712_S517.H5H5YBBXX.s_8.r_1 16289 7 0.843271 1931645
+
+## [1] "Source Name" "cell line" "cell type" "single cell well quality"
+## [5] "genotype" "phenotype" "strain" "spike-in addition"
+## [9] "block" "sum" "detected" "subsets_Mito_sum"
+## [13] "subsets_Mito_detected" "subsets_Mito_percent" "altexps_ERCC_sum" "altexps_ERCC_detected"
+## [17] "altexps_ERCC_percent" "altexps_SIRV_sum" "altexps_SIRV_detected" "altexps_SIRV_percent"
+## [21] "total"
+
+Genes detectados en cada uno de los bloques de secuenciación.
+ + +Genes detectados por cada tratamiento en cada bloque de secuenciación.
+plotColData(sce.416b,
+ x = "block",
+ y = "detected",
+ other_fields = "phenotype"
+) +
+ scale_y_log10() + # Cambiar la escala en Y a logaritmo
+ facet_wrap(~phenotype) + # Dividir tratamiento
+ labs(x = "Bloques de secuenciación", y="Genes detectados \n(log)", title = "Genes detectados") # Etiquetas en X, Y y titulo
Para el filtrado y eliminación de las células de baja calidad emplearemos las variables sum
, detected
, altexps_ERCC_percent
y
+subsets_Mito_percent
.
# --- Fixed thresholds ---
+# Valores de corte fijos y estrictos
+qc.lib <- sce.416b$sum < 100000 # menos de 100 mil cuentas (lecturas)
+qc.nexprs <- sce.416b$detected < 5000 # menos de 5 mil genes
+qc.spike <- sce.416b$altexps_ERCC_percent > 10 # 10 % de las cuentas alineando a ERCC
+qc.mito <- sce.416b$subsets_Mito_percent > 10 # 10 % alineando al genoma mitocondrial
+discard <- qc.lib | qc.nexprs | qc.spike | qc.mito # si falla alguno de estos eliminarlo
+
+# secuencias control = (ERCC spike-in control transcripts)
+
+# Resumen del número de células
+# Número de células eliminadas por cada filtro
+DataFrame(
+ LibSize = sum(qc.lib),
+ NExprs = sum(qc.nexprs),
+ SpikeProp = sum(qc.spike),
+ MitoProp = sum(qc.mito),
+ Total = sum(discard)
+)
## DataFrame with 1 row and 5 columns
+## LibSize NExprs SpikeProp MitoProp Total
+## <integer> <integer> <integer> <integer> <integer>
+## 1 3 0 19 14 33
+Valores de cortes adaptados al comportamiento de nuestros datos.
+# --- Adaptative thresholds ---
+## Usando isOutlier() para determinar los valores de corte
+qc.lib2 <- isOutlier(sce.416b$sum, log = TRUE, type = "lower")
+qc.nexprs2 <- isOutlier(sce.416b$detected,
+ log = TRUE,
+ type = "lower"
+)
+qc.spike2 <- isOutlier(sce.416b$altexps_ERCC_percent,
+ type = "higher"
+)
+qc.mito2 <- isOutlier(sce.416b$subsets_Mito_percent,
+ type = "higher"
+)
+discard2 <- qc.lib2 | qc.nexprs2 | qc.spike2 | qc.mito2
+
+# Extraemos los límites de valores (thresholds)
+attr(qc.lib2, "thresholds")
## lower higher
+## 434082.9 Inf
+
+## lower higher
+## 5231.468 Inf
+## lower higher
+## 5231.468 Inf
+# Obtenemos un resumen del número de células
+# eliminadas por cada filtro
+DataFrame(
+ LibSize = sum(qc.lib2),
+ NExprs = sum(qc.nexprs2),
+ SpikeProp = sum(qc.spike2),
+ MitoProp = sum(qc.mito2),
+ Total = sum(discard2)
+)
## DataFrame with 1 row and 5 columns
+## LibSize NExprs SpikeProp MitoProp Total
+## <integer> <integer> <integer> <integer> <integer>
+## 1 4 0 1 2 6
+## discard2
+## discard FALSE TRUE Sum
+## FALSE 159 0 159
+## TRUE 27 6 33
+## Sum 186 6 192
+isOutlier()
?Supongamos que la mayor parte del conjunto de datos está formado por células de alta calidad 1
+QC = Quality Control +1 But we’ll relax that assumption in a few moments +2 MAD is similar to standard deviation +3 Loosely, 3 MADs will retain 99% of non-outliers values that follow a Normal distribution.
+Figura explicativa, Diapositiva 24.
+Batch
## Determino el bloque (batch) de muestras
+batch <- paste0(sce.416b$phenotype, "-", sce.416b$block)
+
+## Versión de isOutlier() que toma en cuenta los bloques de muestras
+qc.lib3 <- isOutlier(sce.416b$sum,
+ log = TRUE,
+ type = "lower",
+ batch = batch
+)
+qc.nexprs3 <- isOutlier(sce.416b$detected,
+ log = TRUE,
+ type = "lower",
+ batch = batch
+)
+
+qc.spike3 <- isOutlier(sce.416b$altexps_ERCC_percent,
+ type = "higher",
+ batch = batch
+)
+qc.mito3 <- isOutlier(sce.416b$subsets_Mito_percent,
+ type = "higher",
+ batch = batch
+)
+discard3 <- qc.lib3 | qc.nexprs3 | qc.spike3 | qc.mito3
+
+# Extraemos los límites de valores (thresholds)
+attr(qc.lib3, "thresholds")
## induced CBFB-MYH11 oncogene expression-20160113 induced CBFB-MYH11 oncogene expression-20160325
+## lower 461073.1 399133.7
+## higher Inf Inf
+## wild type phenotype-20160113 wild type phenotype-20160325
+## lower 599794.9 370316.5
+## higher Inf Inf
+# Obtenemos un resumen del número de células
+# eliminadas por cada filtro
+DataFrame(
+ LibSize = sum(qc.lib3),
+ NExprs = sum(qc.nexprs3),
+ SpikeProp = sum(qc.spike3),
+ MitoProp = sum(qc.mito3),
+ Total = sum(discard3)
+)
## DataFrame with 1 row and 5 columns
+## LibSize NExprs SpikeProp MitoProp Total
+## <integer> <integer> <integer> <integer> <integer>
+## 1 5 4 6 2 9
+Usaremos el filtro estricto y visualizaremos la distribucion de los datos.
+# Reducir el nombre en los fenotipos / tratamientos
+sce.416b_edited <- sce.416b
+
+sce.416b_edited$phenotype <- ifelse(grepl("induced", sce.416b_edited$phenotype),
+ "induced", "wild type")
+# Verificar
+ unique(sce.416b_edited$phenotype)
## [1] "wild type" "induced"
+ # Agregar columna de genes descartados
+sce.416b_edited$discard <- discard
+
+# Visualizacion grafica
+plotColData(sce.416b_edited, x="block", y="sum", colour_by="discard",
+ other_fields="phenotype") + facet_wrap(~phenotype) +
+ labs(x = "Bloques de secuenciacion", y="Cuentas totales", title = "Numero de cuentas") # Etiquetas en X, Y y titulo
+Visualizacion de todas las variables en una sola grafica.
+# Visualizacion grafica
+gridExtra::grid.arrange(
+ plotColData(sce.416b_edited, x="block", y="sum", colour_by="discard",
+ other_fields="phenotype") + facet_wrap(~phenotype) +
+ scale_y_log10() + labs(x = "Bloques de secuenciacion", y="Cuentas totales", title = "Numero de cuentas"), # Etiquetas en X, Y y titulo
+ plotColData(sce.416b_edited, x="block", y="detected", colour_by="discard",
+ other_fields="phenotype") + facet_wrap(~phenotype) +
+ scale_y_log10() + labs(x = "Bloques de secuenciacion", y="Genes expresados", title = "Numero de genes /n(log)"), # Etiquetas en X, Y y titulo
+ plotColData(sce.416b_edited, x="block", y="subsets_Mito_percent",
+ colour_by="discard", other_fields="phenotype") +
+ facet_wrap(~phenotype) + labs(x = "Bloques de secuenciacion", y="Porcentaje de genes mitocondriales /n(%)", title = "Contenido mitocondrial"), # Etiquetas en X, Y y titulo
+ plotColData(sce.416b_edited, x="block", y="altexps_ERCC_percent",
+ colour_by="discard", other_fields="phenotype") +
+ facet_wrap(~phenotype) + labs(x = "Bloques de secuenciacion", y="Porcentaje de ERCC /n(%)", title = "Contenido de ERCC"), # Etiquetas en X, Y y titulo,
+ ncol=1
+)
Otra grafica de visualizacion grafica es:
+ + ++Descripción gráfica la tecnología Next GEM de 10x Genomics. Fuente: 10x Genomics.
+ +Opciones algorítmicas para detecar los droplets vacíos. Fuente: Lun et al, Genome Biology, 2019.
+El dataset provenie de celulas humanas, senalando un total de 4,340 celulas con buena calidad. Para mas informacion consulta la pagina web en 10Xgenomics.
+A continuacion te presento el reporte de este dataset.
+ +Captura del reporte de secuenciacion del dataset.
+10Xgenomics realiza un alineamiento genomico con STAR y emplea CellRanger para realizar la limpieza y cuantificacion de los datos. Encontrado la informacion dividia en dos carpetas importantes:
+Dejo a tu consideracion la posibilidad de emplear alguno de estos archivos. Ambas carpetas contienen 3 archivos importantes para el analisis de los datos:
+barcode.tsv
= Secuencia que define a cada celula por identificadores de cada celula.feature.tsv
= genes los identfiicadores, se les mio expresionmatrix.mtx
= Union entre las secuencias y genes.+## 2.13.3 Importar datos en R
+Almacenaremos la informacion de 10Xgenomics en la memoria Cache de la computadora y posteriormente lo almacenaremos en una variable.
+# Datos crudos, sin procesar
+bfc <- BiocFileCache()
+raw.path <- bfcrpath(bfc, file.path(
+ "http://cf.10xgenomics.com/samples",
+ "cell-exp/2.1.0/pbmc4k/pbmc4k_raw_gene_bc_matrices.tar.gz"
+
+))
+untar(raw.path, exdir = file.path(tempdir(), "pbmc4k"))
# Archivos contenidos en Temporales
+test_dir <- tempdir()
+test_dir <- paste0(test_dir, "/pbmc4k/raw_gene_bc_matrices/GRCh38")
+list.files(test_dir) # Encontraremos 3 tipos de archivos: barcodes.tsv, genes.tsv y matrix.mtx
## [1] "barcodes.tsv" "genes.tsv" "matrix.mtx"
+# Cargar datos
+fname <- file.path(tempdir(), "pbmc4k/raw_gene_bc_matrices/GRCh38")
+sce.pbmc <- read10xCounts(fname, col.names = TRUE) # Cargar datos de 10Xgenomics de CellRanger
+sce.pbmc
## class: SingleCellExperiment
+## dim: 33694 737280
+## metadata(1): Samples
+## assays(1): counts
+## rownames(33694): ENSG00000243485 ENSG00000237613 ... ENSG00000277475 ENSG00000268674
+## rowData names(2): ID Symbol
+## colnames(737280): AAACCTGAGAAACCAT-1 AAACCTGAGAAACCGC-1 ... TTTGTCATCTTTAGTC-1 TTTGTCATCTTTCCTC-1
+## colData names(2): Sample Barcode
+## reducedDimNames(0):
+## mainExpName: NULL
+## altExpNames(0):
+Contiene 737280 droplets, aun sin filtrar y Evaluando un total de 33,694 genes en este dataset.
+Si quisieras descargar estos mismos datos pero filtrados, te dejo la ruta:
+# Datos filtrados, sin celulas muertas, sin un alto contenido mitocondrial, sin droplets vacios o con multiples celulas
+bfc <- BiocFileCache()
+raw.path <- bfcrpath(bfc, file.path(
+ "http://cf.10xgenomics.com/samples",
+ "cell-exp/2.1.0/pbmc4k/pbmc4k_filtered_gene_bc_matrices.tar.gz"
+
+))
+untar(raw.path, exdir = file.path(tempdir(), "pbmc4k"))
## DataFrame with 737280 rows and 3 columns
+## rank total fitted
+## <numeric> <integer> <numeric>
+## AAACCTGAGAAACCAT-1 219086 1 NA
+## AAACCTGAGAAACCGC-1 504862 0 NA
+## AAACCTGAGAAACCTA-1 219086 1 NA
+## AAACCTGAGAAACGAG-1 504862 0 NA
+## AAACCTGAGAAACGCC-1 219086 1 NA
+## ... ... ... ...
+## TTTGTCATCTTTACAC-1 150232 2 NA
+## TTTGTCATCTTTACGT-1 26415 33 NA
+## TTTGTCATCTTTAGGG-1 504862 0 NA
+## TTTGTCATCTTTAGTC-1 504862 0 NA
+## TTTGTCATCTTTCCTC-1 219086 1 NA
+# Mostremos solo los puntos únicos para acelerar
+# el proceso de hacer esta gráfica
+uniq <- !duplicated(bcrank$rank)
+plot(
+ bcrank$rank[uniq],
+ bcrank$total[uniq],
+ log = "xy",
+ xlab = "Rank",
+ ylab = "Total UMI count",
+ cex.lab = 1.2
+)
## Warning in xy.coords(x, y, xlabel, ylabel, log): 1 y value <= 0 omitted from logarithmic plot
+# Agregarle los puntos de corte
+abline(
+ h = metadata(bcrank)$inflection,
+ col = "darkgreen",
+ lty = 2
+)
+abline(
+ h = metadata(bcrank)$knee,
+ col = "dodgerblue",
+ lty = 2
+)
+legend(
+ "bottomleft",
+ legend = c("Inflection", "Knee"),
+ col = c("darkgreen", "dodgerblue"),
+ lty = 2,
+ cex = 1.2
+)
+UMI = Secuencia unica de cada union en la bead (unique molecular identifiers).
+Encontremos los droplets vacíos usando emptyDrops()
. Los siguientes pasos demoran un tiempo.
## Usemos DropletUtils para encontrar los droplets
+set.seed(100)
+e.out <- emptyDrops(counts(sce.pbmc))
+
+# Revisa ?emptyDrops para una explicación de porque hay valores NA
+summary(e.out$FDR <= 0.001)
## Mode FALSE TRUE NA's
+## logical 1006 4283 731991
+set.seed(100)
+limit <- 100
+all.out <-
+ emptyDrops(counts(sce.pbmc), lower = limit, test.ambient = TRUE)
+
+# Idealmente, este histograma debería verse uniforme.
+# Picos grandes cerca de cero indican que los _barcodes_
+# con un número total de cuentas menor a "lower" no son
+# de origen ambiental.
+hist(all.out$PValue[all.out$Total <= limit &
+ all.out$Total > 0],
+xlab = "P-value",
+main = "",
+col = "grey80"
+)
# Anotación de los genes
+rownames(sce.pbmc) <- uniquifyFeatureNames(
+ rowData(sce.pbmc)$ID, rowData(sce.pbmc)$Symbol
+)
+
+location <- mapIds(EnsDb.Hsapiens.v86,
+ keys = rowData(sce.pbmc)$ID,
+ column = "SEQNAME", keytype = "GENEID"
+)
## Warning: Unable to map 144 of 33694 requested IDs.
+
+# Obtener las estadisticas en un archivo aparte
+stats <- perCellQCMetrics(sce.pbmc,
+ subsets = list(Mito = which(location == "MT"))
+)
+high.mito <- isOutlier(stats$subsets_Mito_percent,
+ type = "higher"
+)
+
+# Eliminar secuencias con alto contenido de genes mitocondriales
+sce.pbmc <- sce.pbmc[, !high.mito]
## [1] "2023-07-31 15:42:42 EDT"
+
+## user system elapsed
+## 351.016 14.223 761.248
+
+## ─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────
+## setting value
+## version R version 4.3.1 (2023-06-16)
+## os macOS Ventura 13.4.1
+## system aarch64, darwin20
+## ui RStudio
+## language (EN)
+## collate en_US.UTF-8
+## ctype en_US.UTF-8
+## tz America/New_York
+## date 2023-07-31
+## rstudio 2023.06.0+421 Mountain Hydrangea (desktop)
+## pandoc 3.1.1 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/ (via rmarkdown)
+##
+## ─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────
+## package * version date (UTC) lib source
+## AnnotationDbi * 1.62.2 2023-07-02 [1] Bioconductor
+## AnnotationFilter * 1.24.0 2023-05-08 [1] Bioconductor
+## AnnotationHub * 3.8.0 2023-05-08 [1] Bioconductor
+## beachmat 2.16.0 2023-05-08 [1] Bioconductor
+## beeswarm 0.4.0 2021-06-01 [1] CRAN (R 4.3.0)
+## Biobase * 2.60.0 2023-05-08 [1] Bioconductor
+## BiocFileCache * 2.8.0 2023-05-08 [1] Bioconductor
+## BiocGenerics * 0.46.0 2023-06-04 [1] Bioconductor
+## BiocIO 1.10.0 2023-05-08 [1] Bioconductor
+## BiocManager 1.30.21.1 2023-07-18 [1] CRAN (R 4.3.0)
+## BiocNeighbors 1.18.0 2023-05-08 [1] Bioconductor
+## BiocParallel 1.34.2 2023-05-28 [1] Bioconductor
+## BiocSingular 1.16.0 2023-05-08 [1] Bioconductor
+## BiocVersion 3.17.1 2022-12-20 [1] Bioconductor
+## biomaRt 2.56.1 2023-06-11 [1] Bioconductor
+## Biostrings 2.68.1 2023-05-21 [1] Bioconductor
+## bit 4.0.5 2022-11-15 [1] CRAN (R 4.3.0)
+## bit64 4.0.5 2020-08-30 [1] CRAN (R 4.3.0)
+## bitops 1.0-7 2021-04-24 [1] CRAN (R 4.3.0)
+## blob 1.2.4 2023-03-17 [1] CRAN (R 4.3.0)
+## bluster 1.10.0 2023-05-08 [1] Bioconductor
+## bookdown 0.34 2023-05-09 [1] CRAN (R 4.3.0)
+## bslib 0.5.0 2023-06-09 [1] CRAN (R 4.3.0)
+## cachem 1.0.8 2023-05-01 [1] CRAN (R 4.3.0)
+## cli 3.6.1 2023-03-23 [1] CRAN (R 4.3.0)
+## cluster 2.1.4 2022-08-22 [1] CRAN (R 4.3.1)
+## codetools 0.2-19 2023-02-01 [1] CRAN (R 4.3.1)
+## colorspace 2.1-0 2023-01-23 [1] CRAN (R 4.3.0)
+## cowplot 1.1.1 2020-12-30 [1] CRAN (R 4.3.0)
+## crayon 1.5.2 2022-09-29 [1] CRAN (R 4.3.0)
+## curl 5.0.1 2023-06-07 [1] CRAN (R 4.3.0)
+## DBI 1.1.3 2022-06-18 [1] CRAN (R 4.3.0)
+## dbplyr * 2.3.3 2023-07-07 [1] CRAN (R 4.3.0)
+## DelayedArray 0.26.6 2023-07-02 [1] Bioconductor
+## DelayedMatrixStats 1.22.0 2023-05-08 [1] Bioconductor
+## digest 0.6.33 2023-07-07 [1] CRAN (R 4.3.0)
+## dplyr * 1.1.2 2023-04-20 [1] CRAN (R 4.3.0)
+## dqrng 0.3.0 2021-05-01 [1] CRAN (R 4.3.0)
+## DropletUtils * 1.20.0 2023-05-08 [1] Bioconductor
+## edgeR 3.42.4 2023-06-04 [1] Bioconductor
+## ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.3.0)
+## EnsDb.Hsapiens.v86 * 2.99.0 2023-07-29 [1] Bioconductor
+## ensembldb * 2.24.0 2023-05-08 [1] Bioconductor
+## evaluate 0.21 2023-05-05 [1] CRAN (R 4.3.0)
+## ExperimentHub 2.8.1 2023-07-16 [1] Bioconductor
+## fansi 1.0.4 2023-01-22 [1] CRAN (R 4.3.0)
+## farver 2.1.1 2022-07-06 [1] CRAN (R 4.3.0)
+## fastmap 1.1.1 2023-02-24 [1] CRAN (R 4.3.0)
+## filelock 1.0.2 2018-10-05 [1] CRAN (R 4.3.0)
+## generics 0.1.3 2022-07-05 [1] CRAN (R 4.3.0)
+## GenomeInfoDb * 1.36.1 2023-07-02 [1] Bioconductor
+## GenomeInfoDbData 1.2.10 2023-06-08 [1] Bioconductor
+## GenomicAlignments 1.36.0 2023-05-08 [1] Bioconductor
+## GenomicFeatures * 1.52.1 2023-07-02 [1] Bioconductor
+## GenomicRanges * 1.52.0 2023-05-08 [1] Bioconductor
+## ggbeeswarm 0.7.2 2023-04-29 [1] CRAN (R 4.3.0)
+## ggplot2 * 3.4.2 2023-04-03 [1] CRAN (R 4.3.0)
+## ggrepel 0.9.3 2023-02-03 [1] CRAN (R 4.3.0)
+## glue 1.6.2 2022-02-24 [1] CRAN (R 4.3.0)
+## gridExtra 2.3 2017-09-09 [1] CRAN (R 4.3.0)
+## gtable 0.3.3 2023-03-21 [1] CRAN (R 4.3.0)
+## HDF5Array 1.28.1 2023-05-08 [1] Bioconductor
+## here 1.0.1 2020-12-13 [1] CRAN (R 4.3.0)
+## highr 0.10 2022-12-22 [1] CRAN (R 4.3.0)
+## hms 1.1.3 2023-03-21 [1] CRAN (R 4.3.0)
+## htmltools 0.5.5 2023-03-23 [1] CRAN (R 4.3.0)
+## httpuv 1.6.11 2023-05-11 [1] CRAN (R 4.3.0)
+## httr 1.4.6 2023-05-08 [1] CRAN (R 4.3.0)
+## igraph 1.5.0.1 2023-07-23 [1] CRAN (R 4.3.0)
+## interactiveDisplayBase 1.38.0 2023-05-08 [1] Bioconductor
+## IRanges * 2.34.1 2023-07-02 [1] Bioconductor
+## irlba 2.3.5.1 2022-10-03 [1] CRAN (R 4.3.0)
+## jquerylib 0.1.4 2021-04-26 [1] CRAN (R 4.3.0)
+## jsonlite 1.8.7 2023-06-29 [1] CRAN (R 4.3.0)
+## KEGGREST 1.40.0 2023-05-08 [1] Bioconductor
+## knitr 1.43 2023-05-25 [1] CRAN (R 4.3.0)
+## labeling 0.4.2 2020-10-20 [1] CRAN (R 4.3.0)
+## later 1.3.1 2023-05-02 [1] CRAN (R 4.3.0)
+## lattice 0.21-8 2023-04-05 [1] CRAN (R 4.3.1)
+## lazyeval 0.2.2 2019-03-15 [1] CRAN (R 4.3.0)
+## lifecycle 1.0.3 2022-10-07 [1] CRAN (R 4.3.0)
+## limma 3.56.2 2023-06-04 [1] Bioconductor
+## locfit 1.5-9.8 2023-06-11 [1] CRAN (R 4.3.0)
+## magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.3.0)
+## Matrix * 1.6-0 2023-07-08 [1] CRAN (R 4.3.0)
+## MatrixGenerics * 1.12.2 2023-07-02 [1] Bioconductor
+## matrixStats * 1.0.0 2023-06-02 [1] CRAN (R 4.3.0)
+## memoise 2.0.1 2021-11-26 [1] CRAN (R 4.3.0)
+## metapod 1.8.0 2023-04-25 [1] Bioconductor
+## mime 0.12 2021-09-28 [1] CRAN (R 4.3.0)
+## munsell 0.5.0 2018-06-12 [1] CRAN (R 4.3.0)
+## pillar 1.9.0 2023-03-22 [1] CRAN (R 4.3.0)
+## pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.3.0)
+## png 0.1-8 2022-11-29 [1] CRAN (R 4.3.0)
+## prettyunits 1.1.1 2020-01-24 [1] CRAN (R 4.3.0)
+## progress 1.2.2 2019-05-16 [1] CRAN (R 4.3.0)
+## promises 1.2.0.1 2021-02-11 [1] CRAN (R 4.3.0)
+## ProtGenerics 1.32.0 2023-05-08 [1] Bioconductor
+## purrr 1.0.1 2023-01-10 [1] CRAN (R 4.3.0)
+## R.methodsS3 1.8.2 2022-06-13 [1] CRAN (R 4.3.0)
+## R.oo 1.25.0 2022-06-12 [1] CRAN (R 4.3.0)
+## R.utils 2.12.2 2022-11-11 [1] CRAN (R 4.3.0)
+## R6 2.5.1 2021-08-19 [1] CRAN (R 4.3.0)
+## rappdirs 0.3.3 2021-01-31 [1] CRAN (R 4.3.0)
+## Rcpp 1.0.11 2023-07-06 [1] CRAN (R 4.3.0)
+## RCurl 1.98-1.12 2023-03-27 [1] CRAN (R 4.3.0)
+## restfulr 0.0.15 2022-06-16 [1] CRAN (R 4.3.0)
+## rhdf5 2.44.0 2023-05-08 [1] Bioconductor
+## rhdf5filters 1.12.1 2023-05-08 [1] Bioconductor
+## Rhdf5lib 1.22.0 2023-05-08 [1] Bioconductor
+## rjson 0.2.21 2022-01-09 [1] CRAN (R 4.3.0)
+## rlang 1.1.1 2023-04-28 [1] CRAN (R 4.3.0)
+## rmarkdown 2.23 2023-07-01 [1] CRAN (R 4.3.0)
+## rprojroot 2.0.3 2022-04-02 [1] CRAN (R 4.3.0)
+## Rsamtools 2.16.0 2023-06-04 [1] Bioconductor
+## RSQLite 2.3.1 2023-04-03 [1] CRAN (R 4.3.0)
+## rstudioapi 0.15.0 2023-07-07 [1] CRAN (R 4.3.0)
+## rsvd 1.0.5 2021-04-16 [1] CRAN (R 4.3.0)
+## rtracklayer 1.60.0 2023-05-08 [1] Bioconductor
+## S4Arrays 1.0.4 2023-05-14 [1] Bioconductor
+## S4Vectors * 0.38.1 2023-05-08 [1] Bioconductor
+## sass 0.4.7 2023-07-15 [1] CRAN (R 4.3.0)
+## ScaledMatrix 1.8.1 2023-05-08 [1] Bioconductor
+## scales 1.2.1 2022-08-20 [1] CRAN (R 4.3.0)
+## scater * 1.28.0 2023-04-25 [1] Bioconductor
+## scran * 1.28.2 2023-07-23 [1] Bioconductor
+## scRNAseq * 2.14.0 2023-04-27 [1] Bioconductor
+## scuttle * 1.9.4 2023-01-23 [1] Bioconductor
+## sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.3.0)
+## shiny 1.7.4.1 2023-07-06 [1] CRAN (R 4.3.0)
+## SingleCellExperiment * 1.22.0 2023-05-08 [1] Bioconductor
+## sparseMatrixStats 1.12.2 2023-07-02 [1] Bioconductor
+## statmod 1.5.0 2023-01-06 [1] CRAN (R 4.3.0)
+## stringi 1.7.12 2023-01-11 [1] CRAN (R 4.3.0)
+## stringr 1.5.0 2022-12-02 [1] CRAN (R 4.3.0)
+## SummarizedExperiment * 1.30.2 2023-06-06 [1] Bioconductor
+## tibble 3.2.1 2023-03-20 [1] CRAN (R 4.3.0)
+## tidyselect 1.2.0 2022-10-10 [1] CRAN (R 4.3.0)
+## utf8 1.2.3 2023-01-31 [1] CRAN (R 4.3.0)
+## vctrs 0.6.3 2023-06-14 [1] CRAN (R 4.3.0)
+## vipor 0.4.5 2017-03-22 [1] CRAN (R 4.3.0)
+## viridis 0.6.4 2023-07-22 [1] CRAN (R 4.3.0)
+## viridisLite 0.4.2 2023-05-02 [1] CRAN (R 4.3.0)
+## withr 2.5.0 2022-03-03 [1] CRAN (R 4.3.0)
+## xfun 0.39 2023-04-20 [1] CRAN (R 4.3.0)
+## XML 3.99-0.14 2023-03-19 [1] CRAN (R 4.3.0)
+## xml2 1.3.5 2023-07-06 [1] CRAN (R 4.3.0)
+## xtable 1.8-4 2019-04-21 [1] CRAN (R 4.3.0)
+## XVector 0.40.0 2023-05-08 [1] Bioconductor
+## yaml 2.3.7 2023-01-23 [1] CRAN (R 4.3.0)
+## zlibbioc 1.46.0 2023-05-08 [1] Bioconductor
+##
+## [1] /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library
+##
+## ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
+
+Dra. Alejandra Medina Rivera
+09 de agosto de 2023
+ +¿Por qué usar Git en vez de solo renombrar los archivos?
+✅✅Por qué es mejor tener una filogenia del archivo.
+++Git es un sistema de control de versiones distribuido, gratuito y de código abierto, diseñado para manejar todo tipo de proyectos, desde los más pequeños hasta los más grandes, con rapidez y eficiencia.
+
++Git es fácil de aprender y ocupa poco espacio con un rendimiento rapidísimo. Supera a las herramientas SCM como Subversion, CVS, Perforce y ClearCase con características como la ramificación local barata, las cómodas áreas de preparación y los múltiples flujos de trabajo.
+
Con Git cada contribuidor tiene una copia del repositorio central, con todos los archivos y la historia de los cambios por los que han pasado.
+++Excuse me, do you have a moment to talk about version control?, Jennifer Bryan, 2017
+
⚠️ NO OLVIDES TENER INSTALADO Git, en caso de que aún no lo hayas instalado, lo puedes descargar en el siguiente enlace https://git-scm.com/downloads.
+Para conocer la localización y la versión de Git que tienes en tu computadora, corre el siguiente comando en la terminal: which git
y git --version
GitHub es una plataforma para guardar proyectos, hace uso de Git. Su principal utilidad es para generar código fuente de programas.
+⚠️ NO OLVIDES TENER UNA CUENTA EN GITHUB, en caso de que aún no lo hayas hecho, puedes ir la página de GitHub y seleccionar join
.
+Es indispensable tu usuario para los ejercicios que siguen.
También existen otras plataformas como Bitbucked y GitLab, las cuales funcionan de manera similar a GitHub.
+which git
(Mac y Linux)where git
(Windows)Git puede comunicarse con un servidor remoto usando uno de dos protocolos, HTTPS o SSH, y cada protocolo usa credenciales diferentes.
La recomendación actual de GitHub es usar HTTPS porque es la manera más fácil de configurar y tiene operabilidad en multiples redes y plataformas.
+credential.helper
almacene en caché su contraseña. (por tanto puedes configurar tu usuario y contraseña en tu equipo de uso)Usualmente cuando inicies un proyecto colaborativo con GitHub inicializa el ropositorio con un README. Copia el HTTPS URL para clonar el repositorio en la terminal git clone https://github.com/TU-USUARIO/TU-REPOSITORIO.git
.
Para usar HTTPS debes crear un token de acceso personal, PAT (PERSONAL ACCESS TOKEN), esa será tu credencial para HTTPS. Es una alternativa al uso de contraseñas para la autenticación en GitHub.
+Como precaución de seguridad, GitHub elimina automáticamente los tokens de acceso personales que no se han usado durante un año.
+¿Cómo crear un token?
+Ve a tu perfil de GitHub, dale click a la imagen de perfil (usualmente en la esquina superior derecha), y busca la opción de settings ó configuración según sea la configuración de idioma que tengas.
Da click a continuación en Developer settings ó Parámetros del desarrollador.
En la barra lateral izquierda da click en Tokens de acceso personal.
Haz click en Generar un nuevo token.
Asígna un nombre descriptivo a tu token.
Selecciona los alcances o permisos que deseas otorgarle a este token. Para usar tu token para acceder a repositorios desde la línea de comando, selecciona repo. (Recomendados: repo, user, workflow )
Finalmente haz click en generar token.
Listo, copia y pega tu token en el lugar dónde siempre lo puedas volver a copiar, ya que por razones de seguridad, una vez salgas de la página no podrás volver a ver el token.
Una vez que tengas un token, puedes ingresarlo en lugar de tu contraseña cuando realices operaciones de Git a través de HTTPS.
+El punto final es que una vez configurada una PAT, varios paquetes de R, incluidos usethis y gh, podrán trabajar con la API de GitHub en su nombre, de forma automática. Por lo tanto, una PAT configurada correctamente significa que todo esto funcionará a la perfección:
+- Operaciones HTTPS remotas a través de la línea de comando Git y, por lo tanto, a través de RStudio
+- Operaciones HTTPS remotas a través del paquete gert R y, por lo tanto, usethis
+- Operaciones de la API de GitHub a través del paquete gh R y, por lo tanto, usethis
+Probar el repositorio Clonado
+Después de hacer clone +Usa estos comandos para verificar tu repositorio y revisar desde dónde se está sincorinzando.
+ +Probemos haciendo un cambio en el README
+ +Qué pasó?
Ahora tenemos que decirle a git que queremos seguir los cambios de ese archivo
Vamos a commit los cambios y luego a subir (push) los mismos a GitHub
Recuerda tu TOKEN!!
+¿Cómo crear un token desde R?
+Puedes ir directamente a la página de GitHub a la parte para generar tu token de acceso personal mediante la siguiente función:
+ +Y con las opciones que se mencionaban anteriormente puedes configurar y crear tu PAT.
+Si lo que quieres es especificar tu PAT en RStudio, las siguientes funciones te serán útiles:
+ + +Para eliminar credenciales utiliza la función credentials::git_credential_forget()
Para lo que sigue a continuación, deberías tener esto:
+mi_repositorio
> Public > YES initialize this repository with a README > click
en el gran botón verde “Create repository”https://github.com/mi_usuario/mi_repositorio.git
. Da click en Create Project.Esto nos generará los siguientes elementos:
+Con este procedimiento ya no es necesario preocuparse por configurar controles remotos Git +y rastrear ramas en la línea de comandos.
+usethis::use_r("titulo_de_un_script")
y observa lo que sucede.PAUSA
+¿Cómo comento y doy push/pull desde RStudio?
+Con la flecha azul podemos hacer pull (RECUERDA HACERLO ANTES DE HACER UN PUSH), y con la flecha verde un push. Para poder comentar y hacer push debemos marcar con una flechita mediante un click
en las pequeñas cajas blancas de la columna Staged, damoc click en commit lo cual no abre la siguiente ventana.
Volvemos a dar click
en commit, y finalizamos con push (flecha verde).
dir.create("mis_imagenes")
en la consola de tu sesión de RStudio (la que está vinculada a tu repositorio). Ejecuta el siguiente código quitando los #
:install.packages("MASS")
+library (MASS)
+data(MASS::cats)
+# pdf("mis_imagenes/cats_plot.pdf")
+ggplot(cats, aes(x = Sex)) +
+ geom_bar(fill = "orange", color = "black") + theme_classic() +
+ xlab("Sexo") + ylab("Número de Gatos") + ggtitle("Gatos")
+# dev.off()
Usa uno de los proyectos que hayas generado en las sesiones anteriores, PERO, que no esté enlazado a GitHub. Ahora veremos como conectar un proyecto de R existente con GitHub. +Realiza los pasos que hicimos en GitHub primero, RStudio después pero asegurate de crear un repositorio con un nuevo nombre.
+Y LISTO!! usa un simple ctrl
+ c
, ó mv
ó click derecho
+ copiar
ó el método que prefieras para mover o copiar archivos. Copia los archivos de tu antigüo proyecto al proyecto nuevo.
+Solo haz commit y push y listo, lo que tenía en tu antigüo proyecto ya está enlazado a GitHub.
Supongamos que tenemos un proyecto de R existente en algún lugar de nuestra computadora.
+NOTA: +Para generar proyecto de RStudio desde la consola puedes utilizar el siguiente código:
+ +O en RStudio con File > New Project > Existing Directory
+Si su proyecto ya es un proyecto de RStudio, ejecútelo.
+¿Ya es un repositorio de Git? La presencia del panel de Git debería alertarlo. Si es así, ha terminado.
+Sino este es el primer camino a seguir:
+usethis::use_git
Si usaste RStudio o usethis, el proyecto debería reiniciarse en RStudio. Hazlo tu mismo si hizo git init. RStudio ahora debería tener un panel Git.
+Si usas el paquete usethis Y has configurado un token de acceso personal (PAT) de GitHub has esto en R:
+ +Esto creará un nuevo repositorio en GitHub, lo agregará como un control remoto, configurará una rama de seguimiento y lo abrirá en su navegador. Lea la ayuda de use_github()
para conocer sus argumentos y consejos sobre cómo configurar una PAT. Esto es extremadamente útil para una variedad de flujos de trabajo que llaman a la API de GitHub. Considere configurar esto si usa usethis, devtools o gh con regularidad.
Volviendo al tema de Proyecto existente, GitHub al final.
+Otra opción que se puede hacer para conectar un proyecto existen a GitHub es ir a hacer un repositorio a GitHub PERO ten en cuenta los siguientes cambios:
+Todo lo demás es igual a los pasos que hacíamos en GitHub primero, RStudio después…
+Ahora ve a tu proyecto de RStudio, has clic en los “dos cuadros de color púrpura y un cuadrado blanco” en el panel de Git. Has clic en “Agregar control remoto”. Pegue la URL aquí y elija un nombre remoto, casi con certeza el origin. Ahora “ADD”.
+ +Pasado esto deberiamos volver en el cuadro de diálogo “New Branch”. Ingresa “master” como el nombre de la rama y asegúrate de que la opción “Sync branch with remote” esté marcada. Haz clic en “Create”. En el siguiente cuadro de diálogo elije “overwrite”.
+ +Ahora solo haz commit/pull/push y cérciorate que FUNCIONE!!
+<<<<<<< HEAD:index.html
+<div id="footer">contact : email.support@github.com</div>
+=======
+<div id="footer">
+ please contact us at support@github.com
+</div>
+>>>>>>> issue-5:index.html
Andrés Arredondo Cruz
+09 de agosto de 2023
+ + +🙇️ Es la información complementaria que el desarrollador escribe sobre una
+función y que se accede con ?
seguido el nombre de una función actual de un
+paquete p.ej. ?unafuncion
.
📁 La documentación se almacena como un archivo .Rd (“R documentation) en la
+carpeta man/
.
🔎 La documentación usa una síntesis especial, que es distinta a la de r y +que esta ligeramente basada en LaTeX.
📄 Se puede renderizar como html, pdf o texto sin formato según se necesite.
En un paquete de r y en cualquier ecosistema de devtools no editamos un
+documento .Rd
manualmente. La documentación usa una síntesis parecida a LaTex
+que puede ser fácil de estropear. Por ventaja existen paquetes como roxigen2.
+Usar roxigen nos permite usar comentarios especiales sobre el inicio de la
+función, esto nos da un par de ventajas:
.Rd
Vamos a crear un función para nuestro paquete. Supongamos que para nuestro +paquete necesitamos una función que calcule la moda. Esta es una forma sencilla +de calcular la moda:
+getmoda <- function(serievector) {
+ uniqv <- unique(serievector)
+ uniqv[which.max(tabulate(match(serievector, uniqv)))]
+}
unique(serievector)
: Crea un vector que contiene únicamente los valores únicos
+de la serie de números serievector.match(serievector, uniqv)
: Encuentra la posición de cada valor de serievector
+en el vector único uniqv.tabulate(match(serievector, uniqv))
: Cuenta cuántas veces aparece cada valor en
+la serie serievector.which.max(tabulate(match(serievector, uniqv)))
: Encuentra el índice del valor
+máximo en el vector de frecuencias.uniqv[which.max(tabulate(match(serievector, uniqv)))]
: Devuelve el valor
+correspondiente al índice calculado, que es la moda.Creamos un ejemplo para ver que funcione:
+serie_numeros <- c(1, 2, 2, 2, 2, 3, 3, 4, 4, 4)
+resultado <- getmoda(serie_numeros)
+print(resultado)
## [1] 2
+Bien ahora si podemos podemos empezar a usar el paquete de roxygen para documentar nuestra función.. comencemos.
+Podemos insertar un esqueleto de comentarios de roxygen para ver su síntesis. Colocamos el cursor en algún lugar de la definición de nuestra función y buscamos en la pestaña Código > Insertar Roxygen Skeleton.
+#' Title
+#'
+#' @param serievector
+#'
+#' @return
+#' @export
+#'
+#' @examples
+
+getmoda <- function(serievector) {
+ uniqv <- unique(serievector)
+ uniqv[which.max(tabulate(match(serievector, uniqv)))]
+}
Ahora ya tenemos un esqueleto de la documentación que nos da una ventaja para su creación.
+Las líneas de comentarios de Roxygen siempre comienzan con #'
, el habitual para un comentario #
mas un '
Veamos los comentarios de uno por uno:
+Empezamos con el titulo. Se sugiere poner en el titulo las acciones principales que realiza la función en este caso por ejemplo podremos usar:
+#' @title Encontrar la Moda de una Serie de Números
+#'
+#' @param serievector
+#'
+#' @return
+#' @export
+#'
+#' @examples
+
+getmoda <- function(serievector) {
+ uniqv <- unique(serievector)
+ uniqv[which.max(tabulate(match(serievector, uniqv)))]
+}
Muy bien!.
+El siguiente comentario que podemos ver es @param
. Pero antes, vamos a añadir una pequeña descripción de la función y como usarla. Primero añadimos la pequeña descripción con @description
:
#' @title Encontrar la Moda de una Serie de Números
+#'
+#' @description Esta función lee una serie de números en forma de vector y
+#' encuentra el elemento que mas se repite, es decir la moda.
+#' @param serievector
+#'
+#' @return
+#' @export
+#'
+#' @examples
+
+getmoda <- function(serievector) {
+ uniqv <- unique(serievector)
+ uniqv[which.max(tabulate(match(serievector, uniqv)))]
+}
Ahora vamos a añadir el comentario @usage
que nos indica como puedes mandar a llamar la función.
#' @title Encontrar la Moda de una Serie de Números
+#'
+#' @description Esta función lee una serie de números en forma de vector y
+#' encuentra el elemento que mas se repite, es decir la moda.
+#' @usage getmoda(serievector)
+#' @param serievector
+#'
+#' @return
+#' @export
+#'
+#' @examples
+
+getmoda <- function(serievector) {
+ uniqv <- unique(serievector)
+ uniqv[which.max(tabulate(match(serievector, uniqv)))]
+}
Ahora si vamos a añadir una pequeña descripción de nuestros argumentos. Si tuviéramos mas de un parámetro en nuestra función podríamos llamar las veces que sea necesario el comentario de parámetro con @param
, veamoslo.
Ahora añadimos una pequeña descripción a nuestro único parámetro que es serievector:
+#' @title Encontrar la Moda de una Serie de Números
+#'
+#' @description Esta función lee una serie de números en forma de vector y
+#' encuentra el elemento que mas se repite, es decir la moda.
+#'
+#' @param serievector Es una serie de números en forma de un vector simple de r.
+#'
+#' @return
+#' @export
+#'
+#' @examples
+
+getmoda <- function(serievector) {
+ uniqv <- unique(serievector)
+ uniqv[which.max(tabulate(match(serievector, uniqv)))]
+}
Después, podemos añadir un comentario de detalles de la función con @details
.
+Por ejemplo, si en nuestro ejemplo tuviéramos ciertos valores no numéricos en nuestro vector de entrada, por ejemplo letras, ¿nuestra función podría leerlas?, o si le diéramos un vector sin caracteres ¿que pasaría?, veamos:
## [1] "d"
+
+## NULL
+Entonces, esto es un ejemplo de lo que podríamos poner en el comentario
+@details
. Hagamoslo describiendo esto. En details podemos agregar detalles
+un poco mas específicos que en la descripción de la función
#' @title Encontrar la Moda de una Serie de Números
+#'
+#' @description Esta función lee una serie de números en forma de vector y
+#' encuentra el elemento que mas se repite, es decir la moda.
+#'
+#' @param serievector Es una serie de números en forma de un vector simple de r.
+#'
+#' @details si tu vector de entrada puede ser interpretado alternando números y
+#' letras escritas entre comillas "". Si un vector esta vacío, dará como
+#' resultado un NULL.
+#' @return
+#' @export
+#'
+#' @examples
+
+getmoda <- function(serievector) {
+ uniqv <- unique(serievector)
+ uniqv[which.max(tabulate(match(serievector, uniqv)))]
+}
Ya casi terminamos de llenar nuestra documentación, pero antes vamos a ver algunos
+otros arrobas que pudieran ser importantes.
+El @import
e @importfrom
importan funciones de otros paquetes en caso
+de que las necesitemos, el primero importa todas las funciones del paquete que
+que solicites, y el segundo importa solo algunas funciones especificas.
+En nuestra función no necesitamos llamar funciones de otros paquetes puesto que
+todas las que usamos están en r base. Pero imaginemos que tu función, por
+ejemplo necesita leer un archivo .tsv
con la función read_tsv del paquete
+readr y después reconvertir la tabla resultante en un archivo con write.table
+pero solo necesitas esa función del paquete utils
, entonces haríamos:
#' @title Encontrar la Moda de una Serie de Números
+#'
+#' @description Esta función lee una serie de números en forma de vector y
+#' encuentra el elemento que mas se repite, es decir la moda.
+#'
+#' @param serievector Es una serie de números en forma de un vector simple de r.
+#'
+#' @details si tu vector de entrada puede ser interpretado alternando números y
+#' letras escritas entre comillas "". Si un vector esta vacío, dará como
+#' resultado un NULL.
+#' @import readr
+#' @importFrom utils write.table
+#' @return
+#' @export
+#'
+#' @examples
+
+getmoda <- function(serievector) {
+ uniqv <- unique(serievector)
+ uniqv[which.max(tabulate(match(serievector, uniqv)))]
+}
Así podemos importar las funciones que necesitemos de otros paquetes y se +incluirán en la documentación y se cargaran automáticamente al cargar tu paquete.
+:eyes::exclamation: Para un correcto funcionamiento de tu paquete y al estar los +paquetes necesarios incluidos en la documentación, no será necesario llamarlos +de la forma ``library(“apackage”)```.
+Entonces llegamos a la sección @return
. Esta descripción le servirá al
+usuario del paquete para conocer cual sera el resultado de la función, que puede
+ser un archivo, una tabla, un numero,etc. Entonces retomando la función que
+usamos al inicio, vamos a escribir una descripción corta del resultado de la
+función getmoda()
.
#' @title Encontrar la Moda de una Serie de Números
+#'
+#' @description Esta función lee una serie de números en forma de vector y
+#' encuentra el elemento que mas se repite, es decir la moda.
+#'
+#' @param serievector Es una serie de números en forma de un vector simple de r.
+#'
+#' @details si tu vector de entrada puede ser interpretado alternando números y
+#' letras escritas entre comillas "". Si un vector esta vacío, dará como
+#' resultado un NULL.
+#' @return El carácter con mas frecuencia de el vector de entrada.
+#' @export
+#'
+#' @examples
+
+getmoda <- function(serievector) {
+ uniqv <- unique(serievector)
+ uniqv[which.max(tabulate(match(serievector, uniqv)))]
+}
Por ultimo tenemos @export
que es el encargado de renderizar la documentación
+para que pueda aparecer en la ventana de Ayuda (abajo a la derecha). esta opción la dejamos para funciones principales que el usuario va a utilizar, aunque puede que existan alguna funciones internas que no queremos que el usuario vea. En ese caso vamos a usar
+@noRd
en lugar de este.
Antes de terminar podemos incluir ejemplos de como funciona nuestra función para +un mejor entendimiento, pongamos los que ya realizamos:
+#' @title Encontrar la Moda de una Serie de Números
+#'
+#' @description Esta función lee una serie de números en forma de vector y
+#' encuentra el elemento que mas se repite, es decir la moda.
+#'
+#' @param serievector Es una serie de números en forma de un vector simple de r.
+#'
+#' @details si tu vector de entrada puede ser interpretado alternando números y
+#' letras escritas entre comillas "". Si un vector esta vacío, dará como
+#' resultado un NULL.
+#' @return El carácter con mas frecuencia de el vector de entrada.
+#' @export
+#'
+#' @examples
+#' serie_números <- c(1, 2, 2, 2, 2, 3, 3, 4, 4, 4)
+#' resultado <- getmoda(serie_números)
+#' print(resultado)
+
+getmoda <- function(serievector) {
+ uniqv <- unique(serievector)
+ uniqv[which.max(tabulate(match(serievector, uniqv)))]
+}
Ahora si, una vez teniendo listo el bloque de comentarios para la documentación, vamos a ejecutar devtools::load_all()
para cargar nuestras funciones y hecho esto, ejecutamos
+devtools::document()
o presionamos Ctrl/Cmd + Shift + D para convertir los comentarios en archivo .Rd
y poder renderizarlo.
💯 Listo, tenemos nuestra documentación para una función. Así se verá cuando el paquete esté terminado.
+2023-07-29
+Les damos la bienvenida al Workshop Creando paquetes de R/Bioconductor para análisis transcriptómicos de célula única!
++
En los últimos años, la generación de datos transcriptómicos de célula única ha cobrado gran relevancia en la investigación biomédica; sin embargo, las herramientas necesarias para su análisis continúan en desarrollo.
+En este taller haremos un repaso a las herramientas estadísticas existentes para analizar datos transcriptómicos de célula única usando Bioconductor, cómo documentar tu análisis y revisaremos algunas herramientas para la interpretación de resultados.
+A la par, aprenderás cuáles son los pasos cruciales para desarrollar un paquete de R y algunas buenas prácticas para la generación de código. Con la integración de estas herramientas, tendrás la oportunidad de crear tu primer paquete y contribuir a la comunidad de desarrolladores.
+Consulta el calendario de este curso en: https://bit.ly/calendarcdsb2023
+Este material posee una licencia tipo Creative Commons Attribution-ShareAlike 4.0 International License.
+Para conocer más sobre esta licencia, visite http://creativecommons.org/licenses/by-sa/4.0/
+ +Instructora: Dra Laura Lucila Gómez Romero
+ +Al igual que con otras tecnologías, single-cell RNA-seq (scRNA-seq) puede presentar diferencias sistemáticas en la cobertura de las librerías de scRAN-seq.
+Típicamente se debe a diferencias técnicas en procesos como la captura de cDNA, la eficiencia de la amplificación por PCR.
+La normalización tiene como objetivo remover estas diferencias sistemáticas para que no interfieran cuando comparamos los perfiles de expresión entre células.
+Al normalizar los datos, las diferencias observadas entre poblaciones celulares o condiciones son debido a la biología y no por factores técnicos.
+La normalización es diferente a la corrección por lote . El proceso de normalización se lleva a cabo sin importar la estructura de los lotes y sólo considera sesgos técnicos, los cuales tienden a afectar los genes en una manera medianamente similar, por ejemplo: la eficiencia de la amplificación por PCR debido al contenido de GC o la longitud del gene.
+“Large single-cell RNA sequencing (scRNA-seq) projects usually need to generate data across multiple batches due to logistical constraints. However, the processing of different batches is often subject to uncontrollable differences, e.g., changes in operator, differences in reagent quality. This results in systematic differences in the observed expression in cells from different batches, which we refer to as “batch effects”. Batch effects are problematic as they can be major drivers of heterogeneity in the data, masking the relevant biological differences and complicating interpretation of the results” -OSCA
+“Normalization occurs regardless of the batch structure and only considers technical biases, while batch correction - as the name suggests - only occurs across batches and must consider both technical biases and biological differences. Technical biases tend to affect genes in a similar manner, or at least in a manner related to their biophysical properties (e.g., length, GC content), while biological differences between batches can be highly unpredictable” -OSCA
+Usaremos el dataset de Zeisel.
+Tipos celulares en cerebro de ratón (oligodendrocitos, microglias, neuronas, etc.)
Procesado con STRT-seq (similar a CEL-seq), un sistema de microfluido.
3005 células y 18441 genes.
Contiene UMIs.
## class: SingleCellExperiment
+## dim: 18441 3005
+## metadata(0):
+## assays(1): counts
+## rownames(18441): ENSMUSG00000029669 ENSMUSG00000046982 ... ENSMUSG00000064337 ENSMUSG00000065947
+## rowData names(2): featureType originalName
+## colnames(3005): 1772071015_C02 1772071017_G12 ... 1772066098_A12 1772058148_F03
+## colData names(10): tissue group # ... level1class level2class
+## reducedDimNames(0):
+## mainExpName: endogenous
+## altExpNames(2): ERCC repeat
+Eliminando células de baja calidad:
+Outliers bajos para log(cuenta total)
Outliers bajos para log(número de features detectados)
Outliers altos para el porcentaje de cuentas de conjuntos de genes especificados (mitocondriales, transcritos spike-in)
La normalización por escalamiento es la estrategia más simple y usada.
+Suposición: Cualquier sesgo específico en cada célula afecta a todos los genes de igual manera a través de escalar por el promedio esperado de cuentas para dicha célula.
+Se realiza de la siguiente manera:
+Estima un “size factor” para cada célula.
Este factor representa el sesgo relativo en esa célula.
Divide las cuentas de los genes de cada célula entre su correspondiente “size factor”.
\[ CuentasNormalizadas = Cuentas / Size factor\]
+Los valores de expresión normalizados pueden ser usados por análisis posteriores como clustering o reducción de dimensiones.
+Tamaño de biblioteca (Library Size): La suma total de las cuentas a tráves de todos los genes en una célula. Asumimos que este valor escala con cualquier sesgo específico en cada célula.
+\[Library Size_{cell} = \sum_{n=1}^{j} gene\] +Donde \(j\) es el número total de genes y \(gene\) es el número de cuentas por gen en la célula \(cell\).
+ +Para escalar los datos ocuparemos un factor de escalamiento llamado Library Size factor. Este valor se calcula a partir de escalar el tamaño de la biblioteca, tal que el promedio de los Library Size factor en todas las células sea igual a 1.
+\[ mean(Library Size factor) = 1 \]
+Lo que nos permite que los valores normalizados están en la misma escala y pueden ser útiles para la interpretación.
+ +# Estimar factores de normalización
+lib.sf.zeisel <- librarySizeFactors(sce.zeisel)
+# Examina la distribución de los tamaños de librerías
+# que acabamos de estimar
+summary(lib.sf.zeisel)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
+## 0.1754 0.5682 0.8669 1.0000 1.2758 4.0651
+
+
+¿A qué crees que se deba la variación en los factores de escalamiento?
+ + +La normalización por tamaño de biblioteca asume que no hay un sesgo de composición, es decir, si existe un grupo de genes cuya expresión aumenta entonces debe existir otro grupo de genes para los que la expresión disminuye en la misma magnitud, esto generalmente no es verdad para los experimentos de scRNA-seq.
Para algunos análisis exploratorios, la precisión de la normalización no es un punto mayor a considerar. La normalizacion por tamaño de biblioteca suele ser suficiente en algunas ocasiones donde se busca identificar clusters y los marcadores de los clusters, dado que el sesgo por composición normalmente no afecta la separación de los clusters, sólo la magnitud.
Sin embargo, esta falta de precisión podría afectar cuando se busca obtener estimaciones y estadísticas a nivel de genes individuales.
El problema causado por el sesgo de composición y cómo eliminarlo ha sido estudiado en bulk RNA-seq por métodos como DESeq2::estimateSizeFactorsFromMatrix()
y edgeR::calcNormFactors()
. En ambos casos, se asume que la mayoría de genes no estarán diferencialmente expresados entre las muestras (en nuestro caso células) y cualquier diferencia entre los genes no diferencialmente expresados representa un sesgo el cual se remueve utilzando el factor de normalización calculado.
Sin embargo, en single-cell RNA-seq se tienen muchas cuentas bajas y ceros debido a limitaciones en la tecnología, las cuales no necesariamente indican ausencia de expresión. Para este fenómento no están preparados los métodos utilizados en bulk RNA-seq.
+La normalización por deconvolución (o normalización scran):
+# Normalización por deconvolución (deconvolution)
+library("scran")
+# Pre-clustering (establece una semilla para obtener resultados reproducibles)
+set.seed(100)
+clust.zeisel <- quickCluster(sce.zeisel)
+
+# Calcula factores de tamaño para la deconvolución (deconvolution)
+deconv.sf.zeisel <-
+ calculateSumFactors(sce.zeisel, clusters = clust.zeisel, min.mean = 0.1)
El pre-clustering mejora la estimación de los size factors al normalizar células similares juntas.
+ +## Min. 1st Qu. Median Mean 3rd Qu. Max.
+## 0.1282 0.4859 0.8248 1.0000 1.3194 4.6521
+
+
+plot(lib.sf.zeisel,
+ deconv.sf.zeisel,
+ xlab = "Library size factor",
+ ylab = "Deconvolution size factor",
+ log = "xy",
+ pch = 16,
+ col=as.integer(factor(sce.zeisel$level1class))
+)
+abline(a = 0, b = 1, col = "red")
La normalización por deconvolución (deconvolution) mejora los resultados para análisis posteriores de una manera más precisa que los métodos para bulk RNA-seq.
+El método scran::calculateSumFactors
algunas veces calcula factores negativos o ceros lo cual altera la matriz de expresión normalizada. ¡Checa los factores que calculas! Si obtienes factores negativos intenta variar el número de clusters.
## Min. 1st Qu. Median Mean 3rd Qu. Max.
+## 0.1282 0.4859 0.8248 1.0000 1.3194 4.6521
+Una vez calculados los factores de normalización con computeSumFactors()
, podemos calcular las cuentas en escala logarítmica usando logNormCounts()
. ¿Qué es lo que hace esta función?
## class: SingleCellExperiment
+## dim: 18441 2815
+## metadata(0):
+## assays(1): counts
+## rownames(18441): ENSMUSG00000029669 ENSMUSG00000046982 ... ENSMUSG00000064337 ENSMUSG00000065947
+## rowData names(2): featureType originalName
+## colnames(2815): 1772071015_C02 1772071017_G12 ... 1772063068_D01 1772066098_A12
+## colData names(10): tissue group # ... level1class level2class
+## reducedDimNames(0):
+## mainExpName: endogenous
+## altExpNames(2): ERCC repeat
+# Normalization
+# Esto agrega el "sizeFactor" como parte del objeto
+# set.seed(100)
+# clust.zeisel <- quickCluster(sce.zeisel)
+sce.zeisel <- computeSumFactors(sce.zeisel, cluster=clust.zeisel, min.mean=0.1)
+
+sce.zeisel
## class: SingleCellExperiment
+## dim: 18441 2815
+## metadata(0):
+## assays(1): counts
+## rownames(18441): ENSMUSG00000029669 ENSMUSG00000046982 ... ENSMUSG00000064337 ENSMUSG00000065947
+## rowData names(2): featureType originalName
+## colnames(2815): 1772071015_C02 1772071017_G12 ... 1772063068_D01 1772066098_A12
+## colData names(11): tissue group # ... level2class sizeFactor
+## reducedDimNames(0):
+## mainExpName: endogenous
+## altExpNames(2): ERCC repeat
+
+## [1] "counts" "logcounts"
+¿Qué gen es más interesante?
+Gen X: el promedio de expresión en el tipo celular A: 50 y B: 10
Gen Y: el promedio de expresión en el tipo celular A: 1100 y B: 1000
## [1] 40
+
+## [1] 100
+
+## [1] 1.609438
+
+## [1] 0.09531018
+Exacto, Gen X
+Te invitamos a leer más sobre otras formas de normalizar, para empezar puedes consultar el siguiente curso del Sanger Institute.
+Si estas interesad@ en diferencias en el contenido total de RNA en cada célula, checa la normalización por spike-ins. La cual asume que los spike-ins fueron añadidos en un nivel constante en cada célula.
+Si tienes resultados donde el library size está asociado a tus datos a pesar de haber normalizado, checa la opción de downsample=TRUE
dentro de la función de logNormCounts()
.
Si te gustaría tener factores de normalización específicos para cada gen que tomen en cuenta la varianza, checa las funciones varianceStabilizingTransformation()
de DESeq2 o sctransform
de Seurat.
Este curso está basado en el libro Orchestrating Single Cell Analysis with Bioconductor de Aaron Lun, Robert Amezquita, Stephanie Hicks y Raphael Gottardo, además del curso de scRNA-seq para WEHI creado por Peter Hickey.
+Y en el material de la comunidadbioinfo/cdsb2020 con el permiso de Leonardo Collado-Torres.
+## [1] "2023-07-31 15:43:10 EDT"
+
+## user system elapsed
+## 376.387 15.493 789.255
+
+## ─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────
+## setting value
+## version R version 4.3.1 (2023-06-16)
+## os macOS Ventura 13.4.1
+## system aarch64, darwin20
+## ui RStudio
+## language (EN)
+## collate en_US.UTF-8
+## ctype en_US.UTF-8
+## tz America/New_York
+## date 2023-07-31
+## rstudio 2023.06.0+421 Mountain Hydrangea (desktop)
+## pandoc 3.1.1 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/ (via rmarkdown)
+##
+## ─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────
+## package * version date (UTC) lib source
+## AnnotationDbi * 1.62.2 2023-07-02 [1] Bioconductor
+## AnnotationFilter * 1.24.0 2023-05-08 [1] Bioconductor
+## AnnotationHub * 3.8.0 2023-05-08 [1] Bioconductor
+## beachmat 2.16.0 2023-05-08 [1] Bioconductor
+## beeswarm 0.4.0 2021-06-01 [1] CRAN (R 4.3.0)
+## Biobase * 2.60.0 2023-05-08 [1] Bioconductor
+## BiocFileCache * 2.8.0 2023-05-08 [1] Bioconductor
+## BiocGenerics * 0.46.0 2023-06-04 [1] Bioconductor
+## BiocIO 1.10.0 2023-05-08 [1] Bioconductor
+## BiocManager 1.30.21.1 2023-07-18 [1] CRAN (R 4.3.0)
+## BiocNeighbors 1.18.0 2023-05-08 [1] Bioconductor
+## BiocParallel 1.34.2 2023-05-28 [1] Bioconductor
+## BiocSingular 1.16.0 2023-05-08 [1] Bioconductor
+## BiocVersion 3.17.1 2022-12-20 [1] Bioconductor
+## biomaRt 2.56.1 2023-06-11 [1] Bioconductor
+## Biostrings 2.68.1 2023-05-21 [1] Bioconductor
+## bit 4.0.5 2022-11-15 [1] CRAN (R 4.3.0)
+## bit64 4.0.5 2020-08-30 [1] CRAN (R 4.3.0)
+## bitops 1.0-7 2021-04-24 [1] CRAN (R 4.3.0)
+## blob 1.2.4 2023-03-17 [1] CRAN (R 4.3.0)
+## bluster 1.10.0 2023-05-08 [1] Bioconductor
+## bookdown 0.34 2023-05-09 [1] CRAN (R 4.3.0)
+## bslib 0.5.0 2023-06-09 [1] CRAN (R 4.3.0)
+## cachem 1.0.8 2023-05-01 [1] CRAN (R 4.3.0)
+## cli 3.6.1 2023-03-23 [1] CRAN (R 4.3.0)
+## cluster 2.1.4 2022-08-22 [1] CRAN (R 4.3.1)
+## codetools 0.2-19 2023-02-01 [1] CRAN (R 4.3.1)
+## colorspace 2.1-0 2023-01-23 [1] CRAN (R 4.3.0)
+## cowplot 1.1.1 2020-12-30 [1] CRAN (R 4.3.0)
+## crayon 1.5.2 2022-09-29 [1] CRAN (R 4.3.0)
+## curl 5.0.1 2023-06-07 [1] CRAN (R 4.3.0)
+## DBI 1.1.3 2022-06-18 [1] CRAN (R 4.3.0)
+## dbplyr * 2.3.3 2023-07-07 [1] CRAN (R 4.3.0)
+## DelayedArray 0.26.6 2023-07-02 [1] Bioconductor
+## DelayedMatrixStats 1.22.0 2023-05-08 [1] Bioconductor
+## digest 0.6.33 2023-07-07 [1] CRAN (R 4.3.0)
+## dplyr * 1.1.2 2023-04-20 [1] CRAN (R 4.3.0)
+## dqrng 0.3.0 2021-05-01 [1] CRAN (R 4.3.0)
+## DropletUtils * 1.20.0 2023-05-08 [1] Bioconductor
+## edgeR 3.42.4 2023-06-04 [1] Bioconductor
+## ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.3.0)
+## EnsDb.Hsapiens.v86 * 2.99.0 2023-07-29 [1] Bioconductor
+## ensembldb * 2.24.0 2023-05-08 [1] Bioconductor
+## evaluate 0.21 2023-05-05 [1] CRAN (R 4.3.0)
+## ExperimentHub 2.8.1 2023-07-16 [1] Bioconductor
+## fansi 1.0.4 2023-01-22 [1] CRAN (R 4.3.0)
+## farver 2.1.1 2022-07-06 [1] CRAN (R 4.3.0)
+## fastmap 1.1.1 2023-02-24 [1] CRAN (R 4.3.0)
+## filelock 1.0.2 2018-10-05 [1] CRAN (R 4.3.0)
+## generics 0.1.3 2022-07-05 [1] CRAN (R 4.3.0)
+## GenomeInfoDb * 1.36.1 2023-07-02 [1] Bioconductor
+## GenomeInfoDbData 1.2.10 2023-06-08 [1] Bioconductor
+## GenomicAlignments 1.36.0 2023-05-08 [1] Bioconductor
+## GenomicFeatures * 1.52.1 2023-07-02 [1] Bioconductor
+## GenomicRanges * 1.52.0 2023-05-08 [1] Bioconductor
+## ggbeeswarm 0.7.2 2023-04-29 [1] CRAN (R 4.3.0)
+## ggplot2 * 3.4.2 2023-04-03 [1] CRAN (R 4.3.0)
+## ggrepel 0.9.3 2023-02-03 [1] CRAN (R 4.3.0)
+## glue 1.6.2 2022-02-24 [1] CRAN (R 4.3.0)
+## gridExtra 2.3 2017-09-09 [1] CRAN (R 4.3.0)
+## gtable 0.3.3 2023-03-21 [1] CRAN (R 4.3.0)
+## HDF5Array 1.28.1 2023-05-08 [1] Bioconductor
+## here 1.0.1 2020-12-13 [1] CRAN (R 4.3.0)
+## highr 0.10 2022-12-22 [1] CRAN (R 4.3.0)
+## hms 1.1.3 2023-03-21 [1] CRAN (R 4.3.0)
+## htmltools 0.5.5 2023-03-23 [1] CRAN (R 4.3.0)
+## httpuv 1.6.11 2023-05-11 [1] CRAN (R 4.3.0)
+## httr 1.4.6 2023-05-08 [1] CRAN (R 4.3.0)
+## igraph 1.5.0.1 2023-07-23 [1] CRAN (R 4.3.0)
+## interactiveDisplayBase 1.38.0 2023-05-08 [1] Bioconductor
+## IRanges * 2.34.1 2023-07-02 [1] Bioconductor
+## irlba 2.3.5.1 2022-10-03 [1] CRAN (R 4.3.0)
+## jquerylib 0.1.4 2021-04-26 [1] CRAN (R 4.3.0)
+## jsonlite 1.8.7 2023-06-29 [1] CRAN (R 4.3.0)
+## KEGGREST 1.40.0 2023-05-08 [1] Bioconductor
+## knitr 1.43 2023-05-25 [1] CRAN (R 4.3.0)
+## labeling 0.4.2 2020-10-20 [1] CRAN (R 4.3.0)
+## later 1.3.1 2023-05-02 [1] CRAN (R 4.3.0)
+## lattice 0.21-8 2023-04-05 [1] CRAN (R 4.3.1)
+## lazyeval 0.2.2 2019-03-15 [1] CRAN (R 4.3.0)
+## lifecycle 1.0.3 2022-10-07 [1] CRAN (R 4.3.0)
+## limma 3.56.2 2023-06-04 [1] Bioconductor
+## locfit 1.5-9.8 2023-06-11 [1] CRAN (R 4.3.0)
+## magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.3.0)
+## Matrix * 1.6-0 2023-07-08 [1] CRAN (R 4.3.0)
+## MatrixGenerics * 1.12.2 2023-07-02 [1] Bioconductor
+## matrixStats * 1.0.0 2023-06-02 [1] CRAN (R 4.3.0)
+## memoise 2.0.1 2021-11-26 [1] CRAN (R 4.3.0)
+## metapod 1.8.0 2023-04-25 [1] Bioconductor
+## mime 0.12 2021-09-28 [1] CRAN (R 4.3.0)
+## munsell 0.5.0 2018-06-12 [1] CRAN (R 4.3.0)
+## pillar 1.9.0 2023-03-22 [1] CRAN (R 4.3.0)
+## pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.3.0)
+## png 0.1-8 2022-11-29 [1] CRAN (R 4.3.0)
+## prettyunits 1.1.1 2020-01-24 [1] CRAN (R 4.3.0)
+## progress 1.2.2 2019-05-16 [1] CRAN (R 4.3.0)
+## promises 1.2.0.1 2021-02-11 [1] CRAN (R 4.3.0)
+## ProtGenerics 1.32.0 2023-05-08 [1] Bioconductor
+## purrr 1.0.1 2023-01-10 [1] CRAN (R 4.3.0)
+## R.methodsS3 1.8.2 2022-06-13 [1] CRAN (R 4.3.0)
+## R.oo 1.25.0 2022-06-12 [1] CRAN (R 4.3.0)
+## R.utils 2.12.2 2022-11-11 [1] CRAN (R 4.3.0)
+## R6 2.5.1 2021-08-19 [1] CRAN (R 4.3.0)
+## rappdirs 0.3.3 2021-01-31 [1] CRAN (R 4.3.0)
+## Rcpp 1.0.11 2023-07-06 [1] CRAN (R 4.3.0)
+## RCurl 1.98-1.12 2023-03-27 [1] CRAN (R 4.3.0)
+## restfulr 0.0.15 2022-06-16 [1] CRAN (R 4.3.0)
+## rhdf5 2.44.0 2023-05-08 [1] Bioconductor
+## rhdf5filters 1.12.1 2023-05-08 [1] Bioconductor
+## Rhdf5lib 1.22.0 2023-05-08 [1] Bioconductor
+## rjson 0.2.21 2022-01-09 [1] CRAN (R 4.3.0)
+## rlang 1.1.1 2023-04-28 [1] CRAN (R 4.3.0)
+## rmarkdown 2.23 2023-07-01 [1] CRAN (R 4.3.0)
+## rprojroot 2.0.3 2022-04-02 [1] CRAN (R 4.3.0)
+## Rsamtools 2.16.0 2023-06-04 [1] Bioconductor
+## RSQLite 2.3.1 2023-04-03 [1] CRAN (R 4.3.0)
+## rstudioapi 0.15.0 2023-07-07 [1] CRAN (R 4.3.0)
+## rsvd 1.0.5 2021-04-16 [1] CRAN (R 4.3.0)
+## rtracklayer 1.60.0 2023-05-08 [1] Bioconductor
+## S4Arrays 1.0.4 2023-05-14 [1] Bioconductor
+## S4Vectors * 0.38.1 2023-05-08 [1] Bioconductor
+## sass 0.4.7 2023-07-15 [1] CRAN (R 4.3.0)
+## ScaledMatrix 1.8.1 2023-05-08 [1] Bioconductor
+## scales 1.2.1 2022-08-20 [1] CRAN (R 4.3.0)
+## scater * 1.28.0 2023-04-25 [1] Bioconductor
+## scran * 1.28.2 2023-07-23 [1] Bioconductor
+## scRNAseq * 2.14.0 2023-04-27 [1] Bioconductor
+## scuttle * 1.9.4 2023-01-23 [1] Bioconductor
+## sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.3.0)
+## shiny 1.7.4.1 2023-07-06 [1] CRAN (R 4.3.0)
+## SingleCellExperiment * 1.22.0 2023-05-08 [1] Bioconductor
+## sparseMatrixStats 1.12.2 2023-07-02 [1] Bioconductor
+## statmod 1.5.0 2023-01-06 [1] CRAN (R 4.3.0)
+## stringi 1.7.12 2023-01-11 [1] CRAN (R 4.3.0)
+## stringr 1.5.0 2022-12-02 [1] CRAN (R 4.3.0)
+## SummarizedExperiment * 1.30.2 2023-06-06 [1] Bioconductor
+## tibble 3.2.1 2023-03-20 [1] CRAN (R 4.3.0)
+## tidyselect 1.2.0 2022-10-10 [1] CRAN (R 4.3.0)
+## utf8 1.2.3 2023-01-31 [1] CRAN (R 4.3.0)
+## vctrs 0.6.3 2023-06-14 [1] CRAN (R 4.3.0)
+## vipor 0.4.5 2017-03-22 [1] CRAN (R 4.3.0)
+## viridis 0.6.4 2023-07-22 [1] CRAN (R 4.3.0)
+## viridisLite 0.4.2 2023-05-02 [1] CRAN (R 4.3.0)
+## withr 2.5.0 2022-03-03 [1] CRAN (R 4.3.0)
+## xfun 0.39 2023-04-20 [1] CRAN (R 4.3.0)
+## XML 3.99-0.14 2023-03-19 [1] CRAN (R 4.3.0)
+## xml2 1.3.5 2023-07-06 [1] CRAN (R 4.3.0)
+## xtable 1.8-4 2019-04-21 [1] CRAN (R 4.3.0)
+## XVector 0.40.0 2023-05-08 [1] Bioconductor
+## yaml 2.3.7 2023-01-23 [1] CRAN (R 4.3.0)
+## zlibbioc 1.46.0 2023-05-08 [1] Bioconductor
+##
+## [1] /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library
+##
+## ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
+
+Instructora: Laura Gómez-Romero
+Contenido adaptado de: aquí
+El siguiente paso en el análisis de scRNA-seq usualmente consiste en identificar grupos de células “similares”.
+Por ejemplo: un análisis de clustering busca identificar células con un perfil transcriptómico similar al calcular distancias entre ellas.
+Si tuviéramos un dataset con dos genes podríamos hacer una gráfica de dos dimensiones para identificar clusters de células.
++ +
+Pero… tenemos decenas de miles de genes : Reducción de dimensionalidad
+Es posible porque la expresión de diferentes genes estaría correlacionada si estos genes son afectados por el mismo proceso biológico.
+Por lo tanto, no necesitamos almacenar información independiente para genes individuales. Podemos comprimir múltiples “features” (genes) en una única dimensión.
+Ventajas:
+library(scRNAseq)
+sce.zeisel <- ZeiselBrainData(ensembl = TRUE)
+
+# Estos datos contienen tipos celulares previamente anotados
+table(sce.zeisel$level1class)
##
+## astrocytes_ependymal endothelial-mural interneurons microglia oligodendrocytes
+## 224 235 290 98 820
+## pyramidal CA1 pyramidal SS
+## 939 399
+Estudio de tipos celulares en el cerebro de ratón (oligodendrocitos, microglia, neuronas, etc) procesados con el sistema STRT-seq (similar a CEL-Seq)
+Descripción aquí
+Zeisel, A. et al. Brain structure. Cell types in the mouse cortex and hippocampus revealed by single-cell RNA-seq. Science 347, 1138-1142 (2015)
+# Quality control
+# Descartar celulas con alto contenido mitocondrial o con alto porcentaje de spike-ins
+library(scater)
+is.mito <- which(rowData(sce.zeisel)$featureType == "mito")
+stats <- perCellQCMetrics(sce.zeisel,
+ subsets = list(Mt = is.mito)
+)
+qc <- quickPerCellQC(stats,
+ percent_subsets = c("altexps_ERCC_percent", "subsets_Mt_percent")
+)
+sce.zeisel <- sce.zeisel[, !qc$discard]
# normalization
+# encontramos unos clusters rápidos para las células y usamos esa información para calcular los factores de tamaño
+library(scran)
+set.seed(1000)
+clusters <- quickCluster(sce.zeisel)
+sce.zeisel <- computeSumFactors(sce.zeisel,
+ cluster = clusters
+)
+sce.zeisel <- logNormCounts(sce.zeisel)
# variance-modelling
+dec.zeisel <- modelGeneVarWithSpikes(
+ sce.zeisel,
+ "ERCC"
+)
+top.zeisel <- getTopHVGs(dec.zeisel, n = 2000)
library(BiocFileCache)
+bfc <- BiocFileCache()
+raw.path <- bfcrpath(bfc, file.path(
+ "http://cf.10xgenomics.com/samples",
+ "cell-exp/2.1.0/pbmc4k/pbmc4k_raw_gene_bc_matrices.tar.gz"
+))
+untar(raw.path, exdir = file.path(tempdir(), "pbmc4k"))
library(DropletUtils)
+library(Matrix)
+fname <- file.path(tempdir(), "pbmc4k/raw_gene_bc_matrices/GRCh38")
+sce.pbmc <- read10xCounts(fname, col.names = TRUE)
Dataset “Células mononucleares humanas de sangre periférica” de 10X Genomics
+Descripción aquí
+Zheng, G. X. Y. et al. Massively parallel digital transcriptional profiling of single cells. Nat. Commun. 8, 14049 (2017)
+# gene-annotation
+library(scater)
+rownames(sce.pbmc) <- uniquifyFeatureNames(
+ rowData(sce.pbmc)$ID, rowData(sce.pbmc)$Symbol
+)
+library(EnsDb.Hsapiens.v86)
+location <- mapIds(EnsDb.Hsapiens.v86,
+ keys = rowData(sce.pbmc)$ID,
+ column = "SEQNAME", keytype = "GENEID"
+)
+
+# cell-detection
+set.seed(100)
+e.out <- emptyDrops(counts(sce.pbmc))
+sce.pbmc <- sce.pbmc[, which(e.out$FDR <= 0.001)]
# quality-control
+stats <- perCellQCMetrics(sce.pbmc,
+ subsets = list(Mito = which(location == "MT"))
+)
+high.mito <- isOutlier(stats$subsets_Mito_percent,
+ type = "higher"
+)
+sce.pbmc <- sce.pbmc[, !high.mito]
+
+# normalization
+library(scran)
+set.seed(1000)
+clusters <- quickCluster(sce.pbmc)
+sce.pbmc <- computeSumFactors(sce.pbmc, cluster = clusters)
+sce.pbmc <- logNormCounts(sce.pbmc)
# variance modelling
+set.seed(1001)
+dec.pbmc <- modelGeneVarByPoisson(sce.pbmc)
+top.pbmc <- getTopHVGs(dec.pbmc, prop = 0.1)
PCA es el arma principal de la reducción de dimensionalidad.
+PCA descubre las combinaciones (lineales) de “features” que capturan la cantidad más grande de variación
+En un PCA, la primer combinación lineal (componente principal) se elige tal que permite capturar la mayor varianza a través de las células. El siguiente PC se elige tal que es “ortogonal” al primero y captura la cantidad más grande de la variación restante, y así sucesivamente…
+Podemos realizar reducción de dimensionalidad al aplicar PCA en la matriz de cuentas transformadas (log-counts matrix) y restringiendo los análisis posteriores a los primeros PCs (top PCs).
+Consideración: Los primeros PCs capturarón “batch effects” (efectos de lote) que afectan muchos genes en una manera coordinada
+## Aplicar PCA en la matriz de cuentas transformadas (log-counts matrix)
+library(scran)
+library(scater)
+set.seed(100)
+sce.zeisel <- runPCA(sce.zeisel,
+ subset_row = top.zeisel
+)
¿Estamos corriendo el análisis sobre todos los genes de nuestro dataset?
+Por default, runPCA() usa un método rápido aproximado que realiza simulaciones, por lo tanto, es necesario ‘configurar la semilla’ para obtener resultados reproducibles.
+Esta elección en análoga a la elección del número de highly variable genes (HGV). Elegir más PCs evitará descartar señal biológica a expensas de retener más ruido
+Ahora exploraremos algunas estrategias guiadas por los datos (data-driven) para hacer esta selección.
+library(PCAtools)
+percent.var <- attr(reducedDim(sce.zeisel), "percentVar")
+chosen.elbow <- PCAtools::findElbowPoint(percent.var)
+plot(percent.var, xlab = "PC", ylab = "Variance explained (%)")
+abline(v = chosen.elbow, col = "red")
Una heurística simple es elegir el número de PCs basado en el porcentaje de varianza explicado por PCs sucesivos.
+Esta es una aproximación heurística más sofisticada que usa el número de clusters como un proxy del número de subpoblaciones.
+choices <- getClusteredPCs(reducedDim(sce.zeisel))
+chosen.clusters <- metadata(choices)$chosen
+
+plot(choices$n.pcs, choices$n.clusters,
+ xlab = "Number of PCs", ylab = "Number of clusters"
+)
+abline(a = 1, b = 1, col = "red")
+abline(v = chosen.clusters, col = "grey80", lty = 2)
Supongamos que esperamos d subpoblaciones de células, en ese caso, necesitamos d-1 dimensiones para garantizar la separación de todas las subpoblaciones.
+Pero… en un escenario real realmente no sabes cuántas poblaciones hay…
+Ventaja: Es una solución pragmática que soluciona el equilibrio sesgo-varianza en los análisis posteriores (especialmente en el análisis de clustering).
+Desventaja: Hace suposiciones fuertes sobre la naturaleza de las diferencias biológicas entre los clusters, y de hecho supone la existencia de clusters, los cuales podrían no existir en algunos procesos biológicos como la diferenciación.
+set.seed(100)
+# Calcula y almacena todos los posibles PCs
+sce.zeisel <- runPCA(sce.zeisel, subset_row = top.zeisel)
+
+# Selecciona y almacena el conjunto reducido de PCs:
+# ... por el método del codo
+reducedDim(sce.zeisel, "PCA_elbow") <- reducedDim(
+ sce.zeisel, "PCA"
+)[, 1:chosen.elbow]
+# ... basado en la estructura poblacional
+reducedDim(sce.zeisel, "PCA_clusters") <- reducedDim(
+ sce.zeisel, "PCA"
+)[, 1:chosen.clusters]
Otra técnica de reducción de dimensiones consiste en conservar todos los PCs hasta que el % de variación explicado alcance algun límite (por ejemplo, basado en la estimación de la variación técnica).
+denoisePCA() automáticamente selecciona el número de PCs.
+Por default, denoisePCA() realiza algunas simulaciones, por lo tanto necesitamos ‘configurar la semilla’ para obtener resultados reproducibles.
+library(scran)
+set.seed(111001001)
+denoised.pbmc <- denoisePCA(sce.pbmc,
+ technical = dec.pbmc, subset.row = top.pbmc
+)
## [1] 3968 8
+La dimensionalidad del output es el límite inferior para el número de PCs requeridos para explicar toda la variación biológica. Lo que significa que cualquier número menor de PCs definitivamente descartará algún aspecto de la señal biológica.
+Esto no grantiza que los PCs retenidos capturen toda la señal biológica
+Esta técnica usualmente retiene más PCs que el método del punto del codo
+scran::denoisePCA() internamente limita el numero de PCs, por default 5-50, para evitar la selección de excesivamente pocos PCs (cuando el ruido técnico es alto relativo al ruido biológico) o excesivamente muchos PCs (cuando el ruido técnico es demasiado bajo).
+dec.pbmc2 <- modelGeneVar(sce.pbmc)
+denoised.pbmc2 <- denoisePCA(sce.pbmc,
+ technical = dec.pbmc2, subset.row = top.pbmc
+)
+dim(reducedDim(denoised.pbmc2))
## [1] 3968 5
+scran::denoisePCA() tiende a funcionar mejor cuando la relación media-varianza refleja el ruido técnico verdadero, i.e estimado por scran::modelGeneVarByPoisson() o scran::modelGeneVarWithSpikes() en vez de scran::modelGeneVar().
+El dataset PBMC está cerca de este límite inferior: el ruido técnico es alto relativo al ruido biológico.
+set.seed(001001001)
+denoised.zeisel <- denoisePCA(sce.zeisel,
+ technical = dec.zeisel, subset.row = top.zeisel
+)
+dim(reducedDim(denoised.zeisel))
## [1] 2815 50
+Los datos de cerebro de Zeisel están cerca de este límite superior: el ruido técnico es demasiado bajo
+ +## [1] 2815 50
+
+## [1] 2815 7
+
+## [1] 2815 15
+Los algoritmos de clustering, así como la mayoría de los algoritmos, operan fácilmente sobre 10-50 (a lo más) PCs, pero ese número es aún demasiado para la visualización.
+Por lo tanto, necesitamos estrategias adicionales para la reducción de dimensionalidad si queremos visualizar los datos.
+PCA es una técnica lineal, por lo tanto, no es eficiente para comprimir diferencias en más de 2 dimensiones en los primeros 2 PCs.
+Ventajas:
+Desventajas:
+set.seed(00101001101)
+sce.zeisel <- runTSNE(sce.zeisel, dimred = "PCA")
+plotReducedDim(sce.zeisel, dimred = "TSNE", colour_by = "level1class")
t-stochastic neighbour embedding (t-SNE) es la visualización por excelencia de datos de scRNA-seq. Intenta encontrar una representación (no-lineal) de los datos usando pocas dimensiones que preserve las distancias entre cada punto y sus vecinos en el espacio multi-dimensional
+ +¿Qué pasa si vuelves a correr runTSNE() sin especificar la semilla?
¿Qué pasa si especificas la semilla pero cambias el valor del parámetro perplexity?
set.seed(100)
+
+sce.zeisel <- runTSNE(sce.zeisel, dimred = "PCA", perplexity = 5)
+p1 <- plotReducedDim(sce.zeisel, dimred = "TSNE", colour_by = "level1class")
+
+sce.zeisel <- runTSNE(sce.zeisel, dimred = "PCA", perplexity = 20)
+p2 <- plotReducedDim(sce.zeisel, dimred = "TSNE", colour_by = "level1class")
+
+sce.zeisel <- runTSNE(sce.zeisel, dimred = "PCA", perplexity = 80)
+p3 <- plotReducedDim(sce.zeisel, dimred = "TSNE", colour_by = "level1class")
+
+library("patchwork")
+p1 + p2 + p3
El siguiente foro discute la selección de parámetros para t-SNE con cierta profundidad.
+Aún así: t-SNE es una herramienta probada para visualización general de datos de scRNA-seq y sigue siendo muy popular
+Uniform manifold approximation and project (UMAP) es una alternativa a t-SNE.
+Así como t-SNE, UMAP intenta encontrar una representación (no lineal) de pocas dimensiones de los datos que preserve las distancias entre cada punto y sus vecinos en el espacio multi-dimensional.
+t-SNE y UMAP están basados en diferentes teorías matemáticas.
+set.seed(1100101001)
+sce.zeisel <- runUMAP(sce.zeisel, dimred = "PCA")
+plotReducedDim(sce.zeisel,
+ dimred = "UMAP",
+ colour_by = "level1class"
+)
Comparado con t-SNE:
+¿Qué pasa si vuelves a correr runUMAP() sin especificar la semilla?
¿Qué pasa si especificas la semilla pero cambias el valor del parámetro n_neighbors?
+ +
+Igual que para t-SNE, es necesario configurar una semilla y diferentes valores para los parámetros cambiaron la visualización.
Si el valor para los parámetros n_neighbors o min_dist es demasiado bajo entonces el ruido aleatorio se interpretará como estructura de alta-resolución, si son demasiado altos entonces se perderá la estructura fina.
TIP: Trata un rango de valores para cada parámetro para asegurarte de que no comprometen ninguna de las conclusiones derivadas de la gráfica UMAP o t-SNE
+Recuerda:
+Reducción de dimensionalidad para la visualización de los datos necesariamente involucra descartar información y distorsionar las distancias entre las células.
No sobre interpretes las gráficas bonitas.
Las gráficas de t-SNE y UMAP son herramientas diagnóstico importantes, por ejemplo: para checar si dos clusters son realmente subclusters vecinos o si un cluster puede ser dividido en más de un cluster.
Es debatible cuál visualización, t-SNE o UMAP, es más útil o estéticamente agradable.
Está bien elegir aquella que funcione mejor para tu análisis (tomando en cuenta que tratarás la gráfica únicamente como una herramienta de visualización/diagnóstico y que no llegarás a ninguna conclusión fuerte basado únicamente en la gráfica).
## [1] "2023-07-31 15:46:07 EDT"
+
+## user system elapsed
+## 544.764 20.333 965.837
+
+## ─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────
+## setting value
+## version R version 4.3.1 (2023-06-16)
+## os macOS Ventura 13.4.1
+## system aarch64, darwin20
+## ui RStudio
+## language (EN)
+## collate en_US.UTF-8
+## ctype en_US.UTF-8
+## tz America/New_York
+## date 2023-07-31
+## rstudio 2023.06.0+421 Mountain Hydrangea (desktop)
+## pandoc 3.1.1 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/ (via rmarkdown)
+##
+## ─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────
+## package * version date (UTC) lib source
+## AnnotationDbi * 1.62.2 2023-07-02 [1] Bioconductor
+## AnnotationFilter * 1.24.0 2023-05-08 [1] Bioconductor
+## AnnotationHub * 3.8.0 2023-05-08 [1] Bioconductor
+## beachmat 2.16.0 2023-05-08 [1] Bioconductor
+## beeswarm 0.4.0 2021-06-01 [1] CRAN (R 4.3.0)
+## Biobase * 2.60.0 2023-05-08 [1] Bioconductor
+## BiocFileCache * 2.8.0 2023-05-08 [1] Bioconductor
+## BiocGenerics * 0.46.0 2023-06-04 [1] Bioconductor
+## BiocIO 1.10.0 2023-05-08 [1] Bioconductor
+## BiocManager 1.30.21.1 2023-07-18 [1] CRAN (R 4.3.0)
+## BiocNeighbors 1.18.0 2023-05-08 [1] Bioconductor
+## BiocParallel 1.34.2 2023-05-28 [1] Bioconductor
+## BiocSingular 1.16.0 2023-05-08 [1] Bioconductor
+## BiocVersion 3.17.1 2022-12-20 [1] Bioconductor
+## biomaRt 2.56.1 2023-06-11 [1] Bioconductor
+## Biostrings 2.68.1 2023-05-21 [1] Bioconductor
+## bit 4.0.5 2022-11-15 [1] CRAN (R 4.3.0)
+## bit64 4.0.5 2020-08-30 [1] CRAN (R 4.3.0)
+## bitops 1.0-7 2021-04-24 [1] CRAN (R 4.3.0)
+## blob 1.2.4 2023-03-17 [1] CRAN (R 4.3.0)
+## bluster 1.10.0 2023-05-08 [1] Bioconductor
+## bookdown 0.34 2023-05-09 [1] CRAN (R 4.3.0)
+## bslib 0.5.0 2023-06-09 [1] CRAN (R 4.3.0)
+## cachem 1.0.8 2023-05-01 [1] CRAN (R 4.3.0)
+## cli 3.6.1 2023-03-23 [1] CRAN (R 4.3.0)
+## cluster 2.1.4 2022-08-22 [1] CRAN (R 4.3.1)
+## codetools 0.2-19 2023-02-01 [1] CRAN (R 4.3.1)
+## colorspace 2.1-0 2023-01-23 [1] CRAN (R 4.3.0)
+## cowplot 1.1.1 2020-12-30 [1] CRAN (R 4.3.0)
+## crayon 1.5.2 2022-09-29 [1] CRAN (R 4.3.0)
+## curl 5.0.1 2023-06-07 [1] CRAN (R 4.3.0)
+## DBI 1.1.3 2022-06-18 [1] CRAN (R 4.3.0)
+## dbplyr * 2.3.3 2023-07-07 [1] CRAN (R 4.3.0)
+## DelayedArray 0.26.6 2023-07-02 [1] Bioconductor
+## DelayedMatrixStats 1.22.0 2023-05-08 [1] Bioconductor
+## digest 0.6.33 2023-07-07 [1] CRAN (R 4.3.0)
+## dplyr * 1.1.2 2023-04-20 [1] CRAN (R 4.3.0)
+## dqrng 0.3.0 2021-05-01 [1] CRAN (R 4.3.0)
+## DropletUtils * 1.20.0 2023-05-08 [1] Bioconductor
+## edgeR 3.42.4 2023-06-04 [1] Bioconductor
+## ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.3.0)
+## EnsDb.Hsapiens.v86 * 2.99.0 2023-07-29 [1] Bioconductor
+## ensembldb * 2.24.0 2023-05-08 [1] Bioconductor
+## evaluate 0.21 2023-05-05 [1] CRAN (R 4.3.0)
+## ExperimentHub 2.8.1 2023-07-16 [1] Bioconductor
+## fansi 1.0.4 2023-01-22 [1] CRAN (R 4.3.0)
+## farver 2.1.1 2022-07-06 [1] CRAN (R 4.3.0)
+## fastmap 1.1.1 2023-02-24 [1] CRAN (R 4.3.0)
+## filelock 1.0.2 2018-10-05 [1] CRAN (R 4.3.0)
+## FNN 1.1.3.2 2023-03-20 [1] CRAN (R 4.3.0)
+## generics 0.1.3 2022-07-05 [1] CRAN (R 4.3.0)
+## GenomeInfoDb * 1.36.1 2023-07-02 [1] Bioconductor
+## GenomeInfoDbData 1.2.10 2023-06-08 [1] Bioconductor
+## GenomicAlignments 1.36.0 2023-05-08 [1] Bioconductor
+## GenomicFeatures * 1.52.1 2023-07-02 [1] Bioconductor
+## GenomicRanges * 1.52.0 2023-05-08 [1] Bioconductor
+## ggbeeswarm 0.7.2 2023-04-29 [1] CRAN (R 4.3.0)
+## ggplot2 * 3.4.2 2023-04-03 [1] CRAN (R 4.3.0)
+## ggrepel * 0.9.3 2023-02-03 [1] CRAN (R 4.3.0)
+## glue 1.6.2 2022-02-24 [1] CRAN (R 4.3.0)
+## gridExtra 2.3 2017-09-09 [1] CRAN (R 4.3.0)
+## gtable 0.3.3 2023-03-21 [1] CRAN (R 4.3.0)
+## HDF5Array 1.28.1 2023-05-08 [1] Bioconductor
+## here 1.0.1 2020-12-13 [1] CRAN (R 4.3.0)
+## highr 0.10 2022-12-22 [1] CRAN (R 4.3.0)
+## hms 1.1.3 2023-03-21 [1] CRAN (R 4.3.0)
+## htmltools 0.5.5 2023-03-23 [1] CRAN (R 4.3.0)
+## httpuv 1.6.11 2023-05-11 [1] CRAN (R 4.3.0)
+## httr 1.4.6 2023-05-08 [1] CRAN (R 4.3.0)
+## igraph 1.5.0.1 2023-07-23 [1] CRAN (R 4.3.0)
+## interactiveDisplayBase 1.38.0 2023-05-08 [1] Bioconductor
+## IRanges * 2.34.1 2023-07-02 [1] Bioconductor
+## irlba 2.3.5.1 2022-10-03 [1] CRAN (R 4.3.0)
+## jquerylib 0.1.4 2021-04-26 [1] CRAN (R 4.3.0)
+## jsonlite 1.8.7 2023-06-29 [1] CRAN (R 4.3.0)
+## KEGGREST 1.40.0 2023-05-08 [1] Bioconductor
+## knitr 1.43 2023-05-25 [1] CRAN (R 4.3.0)
+## labeling 0.4.2 2020-10-20 [1] CRAN (R 4.3.0)
+## later 1.3.1 2023-05-02 [1] CRAN (R 4.3.0)
+## lattice 0.21-8 2023-04-05 [1] CRAN (R 4.3.1)
+## lazyeval 0.2.2 2019-03-15 [1] CRAN (R 4.3.0)
+## lifecycle 1.0.3 2022-10-07 [1] CRAN (R 4.3.0)
+## limma 3.56.2 2023-06-04 [1] Bioconductor
+## locfit 1.5-9.8 2023-06-11 [1] CRAN (R 4.3.0)
+## magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.3.0)
+## Matrix * 1.6-0 2023-07-08 [1] CRAN (R 4.3.0)
+## MatrixGenerics * 1.12.2 2023-07-02 [1] Bioconductor
+## matrixStats * 1.0.0 2023-06-02 [1] CRAN (R 4.3.0)
+## memoise 2.0.1 2021-11-26 [1] CRAN (R 4.3.0)
+## metapod 1.8.0 2023-04-25 [1] Bioconductor
+## mime 0.12 2021-09-28 [1] CRAN (R 4.3.0)
+## munsell 0.5.0 2018-06-12 [1] CRAN (R 4.3.0)
+## patchwork * 1.1.2 2022-08-19 [1] CRAN (R 4.3.0)
+## PCAtools * 2.12.0 2023-05-08 [1] Bioconductor
+## pillar 1.9.0 2023-03-22 [1] CRAN (R 4.3.0)
+## pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.3.0)
+## plyr 1.8.8 2022-11-11 [1] CRAN (R 4.3.0)
+## png 0.1-8 2022-11-29 [1] CRAN (R 4.3.0)
+## prettyunits 1.1.1 2020-01-24 [1] CRAN (R 4.3.0)
+## progress 1.2.2 2019-05-16 [1] CRAN (R 4.3.0)
+## promises 1.2.0.1 2021-02-11 [1] CRAN (R 4.3.0)
+## ProtGenerics 1.32.0 2023-05-08 [1] Bioconductor
+## purrr 1.0.1 2023-01-10 [1] CRAN (R 4.3.0)
+## R.methodsS3 1.8.2 2022-06-13 [1] CRAN (R 4.3.0)
+## R.oo 1.25.0 2022-06-12 [1] CRAN (R 4.3.0)
+## R.utils 2.12.2 2022-11-11 [1] CRAN (R 4.3.0)
+## R6 2.5.1 2021-08-19 [1] CRAN (R 4.3.0)
+## rappdirs 0.3.3 2021-01-31 [1] CRAN (R 4.3.0)
+## Rcpp 1.0.11 2023-07-06 [1] CRAN (R 4.3.0)
+## RCurl 1.98-1.12 2023-03-27 [1] CRAN (R 4.3.0)
+## reshape2 1.4.4 2020-04-09 [1] CRAN (R 4.3.0)
+## restfulr 0.0.15 2022-06-16 [1] CRAN (R 4.3.0)
+## rhdf5 2.44.0 2023-05-08 [1] Bioconductor
+## rhdf5filters 1.12.1 2023-05-08 [1] Bioconductor
+## Rhdf5lib 1.22.0 2023-05-08 [1] Bioconductor
+## rjson 0.2.21 2022-01-09 [1] CRAN (R 4.3.0)
+## rlang 1.1.1 2023-04-28 [1] CRAN (R 4.3.0)
+## rmarkdown 2.23 2023-07-01 [1] CRAN (R 4.3.0)
+## rprojroot 2.0.3 2022-04-02 [1] CRAN (R 4.3.0)
+## Rsamtools 2.16.0 2023-06-04 [1] Bioconductor
+## RSQLite 2.3.1 2023-04-03 [1] CRAN (R 4.3.0)
+## rstudioapi 0.15.0 2023-07-07 [1] CRAN (R 4.3.0)
+## rsvd 1.0.5 2021-04-16 [1] CRAN (R 4.3.0)
+## rtracklayer 1.60.0 2023-05-08 [1] Bioconductor
+## Rtsne 0.16 2022-04-17 [1] CRAN (R 4.3.0)
+## S4Arrays 1.0.4 2023-05-14 [1] Bioconductor
+## S4Vectors * 0.38.1 2023-05-08 [1] Bioconductor
+## sass 0.4.7 2023-07-15 [1] CRAN (R 4.3.0)
+## ScaledMatrix 1.8.1 2023-05-08 [1] Bioconductor
+## scales 1.2.1 2022-08-20 [1] CRAN (R 4.3.0)
+## scater * 1.28.0 2023-04-25 [1] Bioconductor
+## scran * 1.28.2 2023-07-23 [1] Bioconductor
+## scRNAseq * 2.14.0 2023-04-27 [1] Bioconductor
+## scuttle * 1.9.4 2023-01-23 [1] Bioconductor
+## sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.3.0)
+## shiny 1.7.4.1 2023-07-06 [1] CRAN (R 4.3.0)
+## SingleCellExperiment * 1.22.0 2023-05-08 [1] Bioconductor
+## sparseMatrixStats 1.12.2 2023-07-02 [1] Bioconductor
+## statmod 1.5.0 2023-01-06 [1] CRAN (R 4.3.0)
+## stringi 1.7.12 2023-01-11 [1] CRAN (R 4.3.0)
+## stringr 1.5.0 2022-12-02 [1] CRAN (R 4.3.0)
+## SummarizedExperiment * 1.30.2 2023-06-06 [1] Bioconductor
+## tibble 3.2.1 2023-03-20 [1] CRAN (R 4.3.0)
+## tidyselect 1.2.0 2022-10-10 [1] CRAN (R 4.3.0)
+## utf8 1.2.3 2023-01-31 [1] CRAN (R 4.3.0)
+## uwot 0.1.16 2023-06-29 [1] CRAN (R 4.3.0)
+## vctrs 0.6.3 2023-06-14 [1] CRAN (R 4.3.0)
+## vipor 0.4.5 2017-03-22 [1] CRAN (R 4.3.0)
+## viridis 0.6.4 2023-07-22 [1] CRAN (R 4.3.0)
+## viridisLite 0.4.2 2023-05-02 [1] CRAN (R 4.3.0)
+## withr 2.5.0 2022-03-03 [1] CRAN (R 4.3.0)
+## xfun 0.39 2023-04-20 [1] CRAN (R 4.3.0)
+## XML 3.99-0.14 2023-03-19 [1] CRAN (R 4.3.0)
+## xml2 1.3.5 2023-07-06 [1] CRAN (R 4.3.0)
+## xtable 1.8-4 2019-04-21 [1] CRAN (R 4.3.0)
+## XVector 0.40.0 2023-05-08 [1] Bioconductor
+## yaml 2.3.7 2023-01-23 [1] CRAN (R 4.3.0)
+## zlibbioc 1.46.0 2023-05-08 [1] Bioconductor
+##
+## [1] /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library
+##
+## ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
+
+José Antonio Ovando
+08 de agosto de 2023
+ +Ver las diapositivas originales aquí
+Usualmente usamos datos scRNA-seq para caracterizar la heterogeneidad entre células
Para hacer esto, usamos métodos como el clustering y la reducción de dimensionalidad
Esto involucra resumir las diferencias por gen en una sola medida de (dis)similitud entre un par de células
¿Cuáles genes deberíamos usar para calcular esta medida de (dis)similitud?
La elección de los features tiene un mayor impacto en qué tan similares decidimos que son las células
+Deseamos seleccionar los genes altamente variables (High Variable Genes HVGs). Genes con una variación incrementada en comparación con otros genes que están siendo afectados por ruido técnico u otra variación biológica que no es de nuestro interés.
+# Usemos datos de pbmc4k
+library(BiocFileCache)
+bfc <- BiocFileCache()
+raw.path <- bfcrpath(bfc, file.path(
+ "http://cf.10xgenomics.com/samples",
+ "cell-exp/2.1.0/pbmc4k/pbmc4k_raw_gene_bc_matrices.tar.gz"
+))
+untar(raw.path, exdir = file.path(tempdir(), "pbmc4k"))
+
+library(DropletUtils)
+library(Matrix)
+fname <- file.path(tempdir(), "pbmc4k/raw_gene_bc_matrices/GRCh38")
+sce.pbmc <- read10xCounts(fname, col.names = TRUE)
+sce.pbmc
## class: SingleCellExperiment
+## dim: 33694 737280
+## metadata(1): Samples
+## assays(1): counts
+## rownames(33694): ENSG00000243485 ENSG00000237613 ... ENSG00000277475 ENSG00000268674
+## rowData names(2): ID Symbol
+## colnames(737280): AAACCTGAGAAACCAT-1 AAACCTGAGAAACCGC-1 ... TTTGTCATCTTTAGTC-1 TTTGTCATCTTTCCTC-1
+## colData names(2): Sample Barcode
+## reducedDimNames(0):
+## mainExpName: NULL
+## altExpNames(0):
+Dataset “Células mononucleares humanas de sangre periférica” de 10X Genomics
+ +# Anotación de los genes
+library(scater)
+rownames(sce.pbmc) <- uniquifyFeatureNames(
+ rowData(sce.pbmc)$ID, rowData(sce.pbmc)$Symbol
+)
+library(EnsDb.Hsapiens.v86)
+location <- mapIds(EnsDb.Hsapiens.v86,
+ keys = rowData(sce.pbmc)$ID,
+ column = "SEQNAME", keytype = "GENEID"
+)
+
+# Detección de _droplets_ con células
+set.seed(100)
+e.out <- emptyDrops(counts(sce.pbmc))
+sce.pbmc <- sce.pbmc[, which(e.out$FDR <= 0.001)]
# Control de calidad
+stats <- perCellQCMetrics(sce.pbmc,
+ subsets = list(Mito = which(location == "MT"))
+)
+high.mito <- isOutlier(stats$subsets_Mito_percent,
+ type = "higher"
+)
+sce.pbmc <- sce.pbmc[, !high.mito]
+
+# Normalización de los datos
+library(scran)
+set.seed(1000)
+clusters <- quickCluster(sce.pbmc)
+sce.pbmc <- computeSumFactors(sce.pbmc, cluster = clusters)
+sce.pbmc <- logNormCounts(sce.pbmc)
library(scRNAseq)
+sce.416b <- LunSpikeInData(which = "416b")
+sce.416b$block <- factor(sce.416b$block)
Línea celular de células mieloides progenitoras inmortalizadas de ratón usando SmartSeq2 +https://osca.bioconductor.org/lun-416b-cell-line-smart-seq2.html 2
+# gene-annotation
+library(AnnotationHub)
+ens.mm.v97 <- AnnotationHub()[["AH73905"]]
+rowData(sce.416b)$ENSEMBL <- rownames(sce.416b)
+rowData(sce.416b)$SYMBOL <- mapIds(ens.mm.v97,
+ keys = rownames(sce.416b),
+ keytype = "GENEID", column = "SYMBOL"
+)
+rowData(sce.416b)$SEQNAME <- mapIds(ens.mm.v97,
+ keys = rownames(sce.416b),
+ keytype = "GENEID", column = "SEQNAME"
+)
+library(scater)
+rownames(sce.416b) <- uniquifyFeatureNames(
+ rowData(sce.416b)$ENSEMBL,
+ rowData(sce.416b)$SYMBOL
+)
# quality-control
+mito <- which(rowData(sce.416b)$SEQNAME == "MT")
+stats <- perCellQCMetrics(sce.416b, subsets = list(Mt = mito))
+qc <- quickPerCellQC(stats,
+ percent_subsets = c("subsets_Mt_percent", "altexps_ERCC_percent"),
+ batch = sce.416b$block
+)
+sce.416b <- sce.416b[, !qc$discard]
+
+# normalization
+library(scran)
+sce.416b <- computeSumFactors(sce.416b)
+sce.416b <- logNormCounts(sce.416b)
El enfoque más simple para cuantificar la variación per-feature es simplemente calcular la varianza de los log-counts
+# Ordenemos por los genes más interesantes para checar
+# los datos
+dec.pbmc[order(dec.pbmc$bio, decreasing = TRUE), ]
## DataFrame with 33694 rows and 6 columns
+## mean total tech bio p.value FDR
+## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
+## LYZ 1.94129 4.99949 0.824172 4.17532 1.44008e-250 2.83205e-246
+## S100A9 1.92243 4.48842 0.824008 3.66441 1.25315e-193 6.16112e-190
+## S100A8 1.68566 4.35303 0.814526 3.53850 6.89502e-185 2.71195e-181
+## HLA-DRA 2.09063 3.73576 0.819396 2.91636 7.27124e-125 1.19164e-121
+## CD74 2.90080 3.34898 0.786799 2.56218 6.72905e-105 6.96492e-102
+## ... ... ... ... ... ... ...
+## PTMA 3.83144 0.474774 0.733464 -0.258690 0.990672 1
+## HLA-B 4.49970 0.479899 0.748285 -0.268386 0.991626 1
+## EIF1 3.23434 0.471616 0.766366 -0.294750 0.994844 1
+## TMSB4X 6.07476 0.408677 0.719095 -0.310419 0.998006 1
+## B2M 5.94601 0.311508 0.690710 -0.379202 0.999875 1
+El coeficiente de variación de las cuentas al cuadrado (CV2) es una alternativa a la varianza de los log-counts
+👉 Se calcula usando las cuentas en lugar de los log-counts
🤓 CV es el ratio de la desviación estándar a la media y está muy relacionada con el parámetro de dispersión de la distribución binomial negativa usada en edgeR y DESeq2
# Ordenemos por los genes más interesantes para checar
+# los datos
+dec.cv2.pbmc[order(dec.cv2.pbmc$ratio,
+ decreasing = TRUE
+), ]
## DataFrame with 33694 rows and 6 columns
+## mean total trend ratio p.value FDR
+## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
+## GNG11 1.0300722 694.193 1.39084 499.117 0 0
+## PTK7 0.0710992 3664.385 17.32128 211.554 0 0
+## RNF208 0.0721694 3558.419 17.06755 208.490 0 0
+## TUBA8 0.0723199 3544.014 17.03246 208.074 0 0
+## ARHGAP6 0.0777368 3074.612 15.86022 193.857 0 0
+## ... ... ... ... ... ... ...
+## AC023491.2 0 NaN Inf NaN NaN NaN
+## AC233755.2 0 NaN Inf NaN NaN NaN
+## AC233755.1 0 NaN Inf NaN NaN NaN
+## AC213203.1 0 NaN Inf NaN NaN NaN
+## FAM231B 0 NaN Inf NaN NaN NaN
+Generalmente se usa la varianza de los log-counts
+Previamente, ajustamos una línea de tendencia a todos los genes endógenos y asumimos que la mayoría de los genes no están dominados por ruido técnico
En la práctica, todos los genes expresados tienen algún nivel de variabilidad biológica diferente de cero (e.g., transcriptional bursting)
Esto sugiere que nuestros estimados de los componentes técnicos estarán inflados probablemente
👉 Es mejor que pensemos estos estimados como una variación técnica más la variación biológica que no es interesante
🤔 Pero que tal si un grupo de genes a una abundancia particular es afectado por un proceso biológico? +E.g., fuerte sobre regulación de genes específicos de un tipo celular podrían conllevar a un enriquecimiento de HVGs en abundancias altas. Esto inflaría la tendencia y compromete la detección de los genes relevantes
¿Cómo podemos evitar este problema?
+Podemos revisar dos enfoques:
+dec.spike.416b <- modelGeneVarWithSpikes(
+ sce.416b,
+ "ERCC"
+)
+# Ordering by most interesting genes for
+# inspection.
+dec.spike.416b[order(dec.spike.416b$bio,
+ decreasing = TRUE
+), ]
## DataFrame with 46604 rows and 6 columns
+## mean total tech bio p.value FDR
+## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
+## Lyz2 6.61097 13.8497 1.57131 12.2784 1.48993e-186 1.54156e-183
+## Ccl9 6.67846 13.1869 1.50035 11.6866 2.21855e-185 2.19979e-182
+## Top2a 5.81024 14.1787 2.54776 11.6310 3.80015e-65 1.13040e-62
+## Cd200r3 4.83180 15.5613 4.22984 11.3314 9.46221e-24 6.08573e-22
+## Ccnb2 5.97776 13.1393 2.30177 10.8375 3.68706e-69 1.20193e-66
+## ... ... ... ... ... ... ...
+## Rpl5-ps2 3.60625 0.612623 6.32853 -5.71590 0.999616 0.999726
+## Gm11942 3.38768 0.798570 6.51473 -5.71616 0.999459 0.999726
+## Gm12816 2.91276 0.838670 6.57364 -5.73497 0.999422 0.999726
+## Gm13623 2.72844 0.708071 6.45448 -5.74641 0.999544 0.999726
+## Rps12l1 3.15420 0.746615 6.59332 -5.84670 0.999522 0.999726
+plot(dec.spike.416b$mean, dec.spike.416b$total,
+ xlab = "Mean of log-expression",
+ ylab = "Variance of log-expression"
+)
+fit.spike.416b <- metadata(dec.spike.416b)
+points(fit.spike.416b$mean, fit.spike.416b$var,
+ col = "red", pch = 16
+)
+curve(fit.spike.416b$trend(x),
+ col = "dodgerblue",
+ add = TRUE, lwd = 2
+)
set.seed(0010101)
+dec.pois.pbmc <- modelGeneVarByPoisson(sce.pbmc)
+# Ordering by most interesting genes for inspection.
+dec.pois.pbmc[order(dec.pois.pbmc$bio, decreasing = TRUE), ]
## DataFrame with 33694 rows and 6 columns
+## mean total tech bio p.value FDR
+## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
+## LYZ 1.94129 4.99949 0.625290 4.37420 0 0
+## S100A9 1.92243 4.48842 0.628977 3.85944 0 0
+## S100A8 1.68566 4.35303 0.669965 3.68307 0 0
+## HLA-DRA 2.09063 3.73576 0.594291 3.14146 0 0
+## CD74 2.90080 3.34898 0.418950 2.93003 0 0
+## ... ... ... ... ... ... ...
+## COX6A1 0.905824 0.579594 0.633121 -0.0535269 0.883561 0.968946
+## NEDD8 0.842764 0.558542 0.612707 -0.0541655 0.893881 0.978680
+## NDUFA1 0.863339 0.557910 0.619832 -0.0619222 0.920683 1.000000
+## SAP18 0.766806 0.515749 0.582464 -0.0667157 0.946979 1.000000
+## SUMO2 1.362264 0.612538 0.694805 -0.0822664 0.952613 1.000000
+plot(dec.pois.pbmc$mean, dec.pois.pbmc$total,
+ pch = 16, xlab = "Mean of log-expression",
+ ylab = "Variance of log-expression"
+)
+curve(metadata(dec.pois.pbmc)$trend(x),
+ col = "dodgerblue", add = TRUE
+)
Este dataset consiste de células de una línea celular de células inmortalizadas mieloides progenitoras de ratón utilizando SmartSeq2
+Una cantidad constante de spike-in ERCC RNA se agregó a cada lisado celular antes de la prepatación de la librería
+Descripción aquí
+Lun, A. T. L., Calero-Nieto, F. J., Haim-Vilmovsky, L., Göttgens, B. & Marioni, J. C. Assessing the reliability of spike-in normalization for analyses of single-cell RNA sequencing data. Genome Res. 27, 1795–1806 (2017)
+# calculando la variacion por bloque
+dec.block.416b <- modelGeneVarWithSpikes(sce.416b,
+ "ERCC",
+ block = sce.416b$block
+)
+dec.block.416b[order(
+ dec.block.416b$bio,
+ decreasing = TRUE
+), ]
## DataFrame with 46604 rows and 7 columns
+## mean total tech bio p.value FDR
+## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
+## Lyz2 6.61235 13.8619 1.58416 12.2777 0.00000e+00 0.00000e+00
+## Ccl9 6.67841 13.2599 1.44553 11.8143 0.00000e+00 0.00000e+00
+## Top2a 5.81275 14.0192 2.74571 11.2734 3.89855e-137 8.43398e-135
+## Cd200r3 4.83305 15.5909 4.31892 11.2719 1.17783e-54 7.00721e-53
+## Ccnb2 5.97999 13.0256 2.46647 10.5591 1.20380e-151 2.98405e-149
+## ... ... ... ... ... ... ...
+## Gm12816 2.91299 0.842574 6.67730 -5.83472 0.999989 0.999999
+## Gm5786 2.90717 0.879485 6.71686 -5.83738 0.999994 0.999999
+## Rpl9-ps4 3.26421 0.807057 6.64932 -5.84226 0.999988 0.999999
+## Gm13623 2.72788 0.700296 6.63875 -5.93845 0.999998 0.999999
+## Rps12l1 3.15425 0.750775 6.70033 -5.94955 0.999995 0.999999
+## per.block
+## <DataFrame>
+## Lyz2 6.35652:13.3748:2.08227:...:6.86819:14.3490:1.08605:...
+## Ccl9 6.68726:13.0778:1.65923:...:6.66956:13.4420:1.23184:...
+## Top2a 5.34891:17.5972:3.91642:...:6.27659:10.4411:1.57501:...
+## Cd200r3 4.60115:15.7870:5.55587:...:5.06496:15.3948:3.08197:...
+## Ccnb2 5.56701:15.4150:3.46931:...:6.39298:10.6362:1.46362:...
+## ... ...
+## Gm12816 2.86995:0.624143:7.43036:...:2.95604:1.061004:5.92424:...
+## Gm5786 2.96006:0.902872:7.49911:...:2.85427:0.856098:5.93462:...
+## Rpl9-ps4 3.60690:0.543276:7.36805:...:2.92151:1.070839:5.93058:...
+## Gm13623 2.83129:0.852901:7.39442:...:2.62447:0.547692:5.88308:...
+## Rps12l1 3.14399:0.716670:7.57246:...:3.16452:0.784881:5.82819:...
+Al calcular tendencias específicas por batch se tomarán en cuenta las diferencias en la relación media-varianza entre batches
+Se deben obtener estimados de los componentes biológico y técnico para cada gene específicos de cada batch, los cuales se promedian entre los batches para crear una única lista de HVGs
+Hasta ahora hemos ordenado los genes del más al menos interesantemente variable
+¿Qué tanto debemos de bajar en la lista para seleccionar nuestros HVGs?
+Para responder esta pregunta debemos tomar en cuenta lo siguiente: elegir un subset más grande:
+Es difícil determinar el balance óptimo porque el rudio en un contexto podría ser una señal útil en otro contexto
+Discutiremos algunas estrategias para seleccionar HVGs
+La estrategia más simple es seleccionar los top-X genes con los valores más grandes para la métrica relevante de varianza, por ejemplo, la varianza biológica más grande calculada con scran::modelGeneVar()
Pro: El usuario puede controlar directamente el número de HVGs
+Contra: ¿Qué valor de X se debe usar?
+# Works with modelGeneVar() output
+hvg.pbmc.var <- getTopHVGs(dec.pbmc, n = 1000)
+str(hvg.pbmc.var)
## chr [1:1000] "LYZ" "S100A9" "S100A8" "HLA-DRA" "CD74" "CST3" "TYROBP" "IGKC" "CCL5" "S100A4" "HLA-DRB1" "NKG7" ...
+# Works with modelGeneVarWithSpikes() output
+hvg.416b.var <- getTopHVGs(dec.spike.416b, n = 1000)
+str(hvg.416b.var)
## chr [1:1000] "Lyz2" "Ccl9" "Top2a" "Cd200r3" "Ccnb2" "Gm10736" "Hbb-bt" "Mcm5" "Cenpa" "Birc5" "Cd200r4" "Mcm6" ...
+# Also works with modelGeneCV2() but note `var.field`
+hvg.pbmc.cv2 <- getTopHVGs(dec.cv2.pbmc,
+ var.field = "ratio", n = 1000
+)
+str(hvg.pbmc.cv2)
## chr [1:1000] "GNG11" "PTK7" "RNF208" "TUBA8" "ARHGAP6" "PDZK1IP1" "CTB-55O6.8" "NDST2" "GUCY1B3" "MAP1LC3A" "HYI" ...
+Asume que, por ejemplo, no más del 5% de los genes están diferencialmente expresados entre las células de nuestra población:
+Normalmente no conocemos el número de genes diferencialmente expresados desde antes, por lo tanto, solo hemos cambiado un número arbitrario por otro número arbitrario
+RECOMENDACIÓN: Si decides utilizar los top-X HGVs, elige un valor de X y procede con el resto del análisis con la intención de regresar más adelante y probar otros valores, en vez de dedicarle mucho esfuerzo a encontrar el valor óptimo
+Establece un límite fijo en alguna métrica de significancia estadística. Por ejemplo: algunos de los métodos reportan un p-valor para cada gene, entonces selecciona todos los genes con un p-valor ajustado menor que 0.05
+Recuerda que las pruebas estadísticas siempre dependen del tamaño de la muestra
+Ventajas: +* Fácil de implementar +* Menos predecible que la estrategia de los top-X
+Desventajas: +* Podría priorizar genes con significancia estadística fuerte en vez de significancia biológica fuerte
+# Works with modelGeneVar() output
+hvg.pbmc.var.2 <- getTopHVGs(dec.pbmc, fdr.threshold = 0.05)
+str(hvg.pbmc.var.2)
## chr [1:629] "LYZ" "S100A9" "S100A8" "HLA-DRA" "CD74" "CST3" "TYROBP" "IGKC" "CCL5" "S100A4" "HLA-DRB1" "NKG7" ...
+# Works with modelGeneVarWithSpikes() output
+hvg.416b.var.2 <- getTopHVGs(dec.spike.416b,
+ fdr.threshold = 0.05
+)
+str(hvg.416b.var.2)
## chr [1:2568] "Lyz2" "Ccl9" "Top2a" "Cd200r3" "Ccnb2" "Gm10736" "Hbb-bt" "Mcm5" "Cenpa" "Birc5" "Cd200r4" "Mcm6" ...
+# Also works with modelGeneCV2() but note `var.field`
+hvg.pbmc.cv2.2 <- getTopHVGs(dec.cv2.pbmc,
+ var.field = "ratio", fdr.threshold = 0.05
+)
+str(hvg.pbmc.cv2.2)
## chr [1:1706] "GNG11" "PTK7" "RNF208" "TUBA8" "ARHGAP6" "PDZK1IP1" "CTB-55O6.8" "NDST2" "GUCY1B3" "MAP1LC3A" "HYI" ...
+Selecciona todos los genes con una varianza biológica positiva
+Este es un extremo del equilibrio sesgo-varianza que minimiza el sesgo con el costo de maximizar el ruido
+Si seguimos esta aproximación, estamos:
+Funciona mejor si tenemos datasets altamente heterogeneos que contienen muchos tipos celulares diferentes
+# Works with modelGeneVar() output
+hvg.pbmc.var.3 <- getTopHVGs(dec.pbmc, var.threshold = 0)
+str(hvg.pbmc.var.3)
## chr [1:12852] "LYZ" "S100A9" "S100A8" "HLA-DRA" "CD74" "CST3" "TYROBP" "IGKC" "CCL5" "S100A4" "HLA-DRB1" "NKG7" ...
+# Works with modelGeneVarWithSpikes() output
+hvg.416b.var.3 <- getTopHVGs(dec.spike.416b,
+ var.threshold = 0
+)
+str(hvg.416b.var.3)
## chr [1:11087] "Lyz2" "Ccl9" "Top2a" "Cd200r3" "Ccnb2" "Gm10736" "Hbb-bt" "Mcm5" "Cenpa" "Birc5" "Cd200r4" "Mcm6" ...
+# Also works with modelGeneCV2() but note `var.field` and
+# value of `var.threshold`
+hvg.pbmc.cv2.3 <- getTopHVGs(dec.cv2.pbmc,
+ var.field = "ratio", var.threshold = 1
+)
+str(hvg.pbmc.cv2.2)
## chr [1:1706] "GNG11" "PTK7" "RNF208" "TUBA8" "ARHGAP6" "PDZK1IP1" "CTB-55O6.8" "NDST2" "GUCY1B3" "MAP1LC3A" "HYI" ...
+Para este ejercicio tendrás que repetir la gráfica que muestra la tendencia de la relación media-varianza (ejeX: media de la expresión, ejeY: varianza de la expresión) incluyendo la línea de tendencia obtenida con alguna de las funciones vistas en la primer parte de la clase (modelGeneVar, modelGeneVarWithSpikes, modelGeneCV2). En esta gráfica, deberás colorear los puntos que corresponden a los HGVs obtenidos con algunos de los enfoques revisados
+RESPUESTA
+ + +Una estrategia contundente es usar sets predefinidos de genes de interés. No hay vergüenza en aprovechar el conocimiento biológivo previo
+Sin embargo, limita nuestra capacidad de descubrir aspectos nuevos o inesperados de la variación. Por lo tanto, considera esta como una estrategia complementaria a otros tipo de estrategias de selección de HGVs
+También podrías eliminar listas pre-definidas de genes:
+Sin embargo, tampoco hay que pecar de prevacido, espera a que estos genes demuestren ser problemáticos para removerlos
+# Elegimos el 10% de los genes con con componente biologico de variacion mayor
+dec.pbmc <- modelGeneVar(sce.pbmc)
+chosen <- getTopHVGs(dec.pbmc, prop = 0.1)
+str(chosen)
## chr [1:1285] "LYZ" "S100A9" "S100A8" "HLA-DRA" "CD74" "CST3" "TYROBP" "IGKC" "CCL5" "S100A4" "HLA-DRB1" "NKG7" ...
+Después de esto tenemos varias opciones para imponer nuestra selección de HGVs durante el resto del análisis:
+## class: SingleCellExperiment
+## dim: 1285 3968
+## metadata(1): Samples
+## assays(2): counts logcounts
+## rownames(1285): LYZ S100A9 ... ZNF770 TMEM106C
+## rowData names(2): ID Symbol
+## colnames(3968): AAACCTGAGACAGACC-1 AAACCTGAGGCATGGT-1 ... TTTGTCAGTTAAGACA-1 TTTGTCATCCCAAGAT-1
+## colData names(3): Sample Barcode sizeFactor
+## reducedDimNames(0):
+## mainExpName: NULL
+## altExpNames(0):
+PRO: Te aseguras de que los métodos posteriores sólo usen estos genes para sus cálculos
+CONTRA: Los genes no-HGVs son eliminados del nuevo objeto SingleCellExperiment, lo cual hace menos conveniente la interrogación del dataset completo sobre genes que no son HGVs
+# Example of specifying HVGs in a downstream function
+# Performing PCA only on the chosen HVGs.
+library(scater)
+sce.pbmc <- runPCA(sce.pbmc, subset_row = chosen)
+sce.pbmc
## class: SingleCellExperiment
+## dim: 33694 3968
+## metadata(1): Samples
+## assays(2): counts logcounts
+## rownames(33694): RP11-34P13.3 FAM138A ... AC213203.1 FAM231B
+## rowData names(2): ID Symbol
+## colnames(3968): AAACCTGAGACAGACC-1 AAACCTGAGGCATGGT-1 ... TTTGTCAGTTAAGACA-1 TTTGTCATCCCAAGAT-1
+## colData names(3): Sample Barcode sizeFactor
+## reducedDimNames(1): PCA
+## mainExpName: NULL
+## altExpNames(0):
+Mantiene el objeto SingleCellExperiment original y especifica los genes para usar en funciones posteriores mediante un argumento adicional como subset_row
+PRO: Es útil si el análisis usa varios conjuntos de HGVs en diferentes pasos
+CONTRA: Podría ser inconveniente especificar repetidamente el mismo conjunto de HGVs en diferentes pasos
+# Add the full SCE to the subsetted data SCE
+altExp(sce.pbmc.hvg, "original") <- sce.pbmc
+sce.pbmc.hvg
## class: SingleCellExperiment
+## dim: 1285 3968
+## metadata(1): Samples
+## assays(2): counts logcounts
+## rownames(1285): LYZ S100A9 ... ZNF770 TMEM106C
+## rowData names(2): ID Symbol
+## colnames(3968): AAACCTGAGACAGACC-1 AAACCTGAGGCATGGT-1 ... TTTGTCAGTTAAGACA-1 TTTGTCATCCCAAGAT-1
+## colData names(3): Sample Barcode sizeFactor
+## reducedDimNames(0):
+## mainExpName: NULL
+## altExpNames(1): original
+
+## class: SingleCellExperiment
+## dim: 33694 3968
+## metadata(1): Samples
+## assays(2): counts logcounts
+## rownames(33694): RP11-34P13.3 FAM138A ... AC213203.1 FAM231B
+## rowData names(2): ID Symbol
+## colnames(3968): AAACCTGAGACAGACC-1 AAACCTGAGGCATGGT-1 ... TTTGTCAGTTAAGACA-1 TTTGTCATCCCAAGAT-1
+## colData names(3): Sample Barcode sizeFactor
+## reducedDimNames(1): PCA
+## mainExpName: NULL
+## altExpNames(0):
+Utilizando el sistema de “experimento alternartivo” en la clase SingleCellExperiment
+PRO: Evita algunos problemas a largo plazo cuando el dataset original no está sincronizado con el conjunto filtrado por HVGs
+CONTRA: Ralentiza todos los análisis subsecuentes
+Para CEL-Seq2:
+scran::modelGeneVarWithSpikes()
Para 10X:
+scran::modelGeneVarByPoisson()
Si quieres irte por el lado de conservar demasiados genes:
+scran::getTopHVGs(dec, var.threshold=0)
Y realiza una comparación rápida con, por lo menos, el top-1000 HVGs
+Regresa al paso de selección de HVG para eliminar genes problemáticos tantas veces como sea necesario
+## [1] "2023-07-31 15:44:00 EDT"
+
+## user system elapsed
+## 422.764 17.018 839.118
+
+## ─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────
+## setting value
+## version R version 4.3.1 (2023-06-16)
+## os macOS Ventura 13.4.1
+## system aarch64, darwin20
+## ui RStudio
+## language (EN)
+## collate en_US.UTF-8
+## ctype en_US.UTF-8
+## tz America/New_York
+## date 2023-07-31
+## rstudio 2023.06.0+421 Mountain Hydrangea (desktop)
+## pandoc 3.1.1 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/ (via rmarkdown)
+##
+## ─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────
+## package * version date (UTC) lib source
+## AnnotationDbi * 1.62.2 2023-07-02 [1] Bioconductor
+## AnnotationFilter * 1.24.0 2023-05-08 [1] Bioconductor
+## AnnotationHub * 3.8.0 2023-05-08 [1] Bioconductor
+## beachmat 2.16.0 2023-05-08 [1] Bioconductor
+## beeswarm 0.4.0 2021-06-01 [1] CRAN (R 4.3.0)
+## Biobase * 2.60.0 2023-05-08 [1] Bioconductor
+## BiocFileCache * 2.8.0 2023-05-08 [1] Bioconductor
+## BiocGenerics * 0.46.0 2023-06-04 [1] Bioconductor
+## BiocIO 1.10.0 2023-05-08 [1] Bioconductor
+## BiocManager 1.30.21.1 2023-07-18 [1] CRAN (R 4.3.0)
+## BiocNeighbors 1.18.0 2023-05-08 [1] Bioconductor
+## BiocParallel 1.34.2 2023-05-28 [1] Bioconductor
+## BiocSingular 1.16.0 2023-05-08 [1] Bioconductor
+## BiocVersion 3.17.1 2022-12-20 [1] Bioconductor
+## biomaRt 2.56.1 2023-06-11 [1] Bioconductor
+## Biostrings 2.68.1 2023-05-21 [1] Bioconductor
+## bit 4.0.5 2022-11-15 [1] CRAN (R 4.3.0)
+## bit64 4.0.5 2020-08-30 [1] CRAN (R 4.3.0)
+## bitops 1.0-7 2021-04-24 [1] CRAN (R 4.3.0)
+## blob 1.2.4 2023-03-17 [1] CRAN (R 4.3.0)
+## bluster 1.10.0 2023-05-08 [1] Bioconductor
+## bookdown 0.34 2023-05-09 [1] CRAN (R 4.3.0)
+## bslib 0.5.0 2023-06-09 [1] CRAN (R 4.3.0)
+## cachem 1.0.8 2023-05-01 [1] CRAN (R 4.3.0)
+## cli 3.6.1 2023-03-23 [1] CRAN (R 4.3.0)
+## cluster 2.1.4 2022-08-22 [1] CRAN (R 4.3.1)
+## codetools 0.2-19 2023-02-01 [1] CRAN (R 4.3.1)
+## colorspace 2.1-0 2023-01-23 [1] CRAN (R 4.3.0)
+## cowplot 1.1.1 2020-12-30 [1] CRAN (R 4.3.0)
+## crayon 1.5.2 2022-09-29 [1] CRAN (R 4.3.0)
+## curl 5.0.1 2023-06-07 [1] CRAN (R 4.3.0)
+## DBI 1.1.3 2022-06-18 [1] CRAN (R 4.3.0)
+## dbplyr * 2.3.3 2023-07-07 [1] CRAN (R 4.3.0)
+## DelayedArray 0.26.6 2023-07-02 [1] Bioconductor
+## DelayedMatrixStats 1.22.0 2023-05-08 [1] Bioconductor
+## digest 0.6.33 2023-07-07 [1] CRAN (R 4.3.0)
+## dplyr * 1.1.2 2023-04-20 [1] CRAN (R 4.3.0)
+## dqrng 0.3.0 2021-05-01 [1] CRAN (R 4.3.0)
+## DropletUtils * 1.20.0 2023-05-08 [1] Bioconductor
+## edgeR 3.42.4 2023-06-04 [1] Bioconductor
+## ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.3.0)
+## EnsDb.Hsapiens.v86 * 2.99.0 2023-07-29 [1] Bioconductor
+## ensembldb * 2.24.0 2023-05-08 [1] Bioconductor
+## evaluate 0.21 2023-05-05 [1] CRAN (R 4.3.0)
+## ExperimentHub 2.8.1 2023-07-16 [1] Bioconductor
+## fansi 1.0.4 2023-01-22 [1] CRAN (R 4.3.0)
+## farver 2.1.1 2022-07-06 [1] CRAN (R 4.3.0)
+## fastmap 1.1.1 2023-02-24 [1] CRAN (R 4.3.0)
+## filelock 1.0.2 2018-10-05 [1] CRAN (R 4.3.0)
+## generics 0.1.3 2022-07-05 [1] CRAN (R 4.3.0)
+## GenomeInfoDb * 1.36.1 2023-07-02 [1] Bioconductor
+## GenomeInfoDbData 1.2.10 2023-06-08 [1] Bioconductor
+## GenomicAlignments 1.36.0 2023-05-08 [1] Bioconductor
+## GenomicFeatures * 1.52.1 2023-07-02 [1] Bioconductor
+## GenomicRanges * 1.52.0 2023-05-08 [1] Bioconductor
+## ggbeeswarm 0.7.2 2023-04-29 [1] CRAN (R 4.3.0)
+## ggplot2 * 3.4.2 2023-04-03 [1] CRAN (R 4.3.0)
+## ggrepel 0.9.3 2023-02-03 [1] CRAN (R 4.3.0)
+## glue 1.6.2 2022-02-24 [1] CRAN (R 4.3.0)
+## gridExtra 2.3 2017-09-09 [1] CRAN (R 4.3.0)
+## gtable 0.3.3 2023-03-21 [1] CRAN (R 4.3.0)
+## HDF5Array 1.28.1 2023-05-08 [1] Bioconductor
+## here 1.0.1 2020-12-13 [1] CRAN (R 4.3.0)
+## highr 0.10 2022-12-22 [1] CRAN (R 4.3.0)
+## hms 1.1.3 2023-03-21 [1] CRAN (R 4.3.0)
+## htmltools 0.5.5 2023-03-23 [1] CRAN (R 4.3.0)
+## httpuv 1.6.11 2023-05-11 [1] CRAN (R 4.3.0)
+## httr 1.4.6 2023-05-08 [1] CRAN (R 4.3.0)
+## igraph 1.5.0.1 2023-07-23 [1] CRAN (R 4.3.0)
+## interactiveDisplayBase 1.38.0 2023-05-08 [1] Bioconductor
+## IRanges * 2.34.1 2023-07-02 [1] Bioconductor
+## irlba 2.3.5.1 2022-10-03 [1] CRAN (R 4.3.0)
+## jquerylib 0.1.4 2021-04-26 [1] CRAN (R 4.3.0)
+## jsonlite 1.8.7 2023-06-29 [1] CRAN (R 4.3.0)
+## KEGGREST 1.40.0 2023-05-08 [1] Bioconductor
+## knitr 1.43 2023-05-25 [1] CRAN (R 4.3.0)
+## labeling 0.4.2 2020-10-20 [1] CRAN (R 4.3.0)
+## later 1.3.1 2023-05-02 [1] CRAN (R 4.3.0)
+## lattice 0.21-8 2023-04-05 [1] CRAN (R 4.3.1)
+## lazyeval 0.2.2 2019-03-15 [1] CRAN (R 4.3.0)
+## lifecycle 1.0.3 2022-10-07 [1] CRAN (R 4.3.0)
+## limma 3.56.2 2023-06-04 [1] Bioconductor
+## locfit 1.5-9.8 2023-06-11 [1] CRAN (R 4.3.0)
+## magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.3.0)
+## Matrix * 1.6-0 2023-07-08 [1] CRAN (R 4.3.0)
+## MatrixGenerics * 1.12.2 2023-07-02 [1] Bioconductor
+## matrixStats * 1.0.0 2023-06-02 [1] CRAN (R 4.3.0)
+## memoise 2.0.1 2021-11-26 [1] CRAN (R 4.3.0)
+## metapod 1.8.0 2023-04-25 [1] Bioconductor
+## mime 0.12 2021-09-28 [1] CRAN (R 4.3.0)
+## munsell 0.5.0 2018-06-12 [1] CRAN (R 4.3.0)
+## pillar 1.9.0 2023-03-22 [1] CRAN (R 4.3.0)
+## pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.3.0)
+## png 0.1-8 2022-11-29 [1] CRAN (R 4.3.0)
+## prettyunits 1.1.1 2020-01-24 [1] CRAN (R 4.3.0)
+## progress 1.2.2 2019-05-16 [1] CRAN (R 4.3.0)
+## promises 1.2.0.1 2021-02-11 [1] CRAN (R 4.3.0)
+## ProtGenerics 1.32.0 2023-05-08 [1] Bioconductor
+## purrr 1.0.1 2023-01-10 [1] CRAN (R 4.3.0)
+## R.methodsS3 1.8.2 2022-06-13 [1] CRAN (R 4.3.0)
+## R.oo 1.25.0 2022-06-12 [1] CRAN (R 4.3.0)
+## R.utils 2.12.2 2022-11-11 [1] CRAN (R 4.3.0)
+## R6 2.5.1 2021-08-19 [1] CRAN (R 4.3.0)
+## rappdirs 0.3.3 2021-01-31 [1] CRAN (R 4.3.0)
+## Rcpp 1.0.11 2023-07-06 [1] CRAN (R 4.3.0)
+## RCurl 1.98-1.12 2023-03-27 [1] CRAN (R 4.3.0)
+## restfulr 0.0.15 2022-06-16 [1] CRAN (R 4.3.0)
+## rhdf5 2.44.0 2023-05-08 [1] Bioconductor
+## rhdf5filters 1.12.1 2023-05-08 [1] Bioconductor
+## Rhdf5lib 1.22.0 2023-05-08 [1] Bioconductor
+## rjson 0.2.21 2022-01-09 [1] CRAN (R 4.3.0)
+## rlang 1.1.1 2023-04-28 [1] CRAN (R 4.3.0)
+## rmarkdown 2.23 2023-07-01 [1] CRAN (R 4.3.0)
+## rprojroot 2.0.3 2022-04-02 [1] CRAN (R 4.3.0)
+## Rsamtools 2.16.0 2023-06-04 [1] Bioconductor
+## RSQLite 2.3.1 2023-04-03 [1] CRAN (R 4.3.0)
+## rstudioapi 0.15.0 2023-07-07 [1] CRAN (R 4.3.0)
+## rsvd 1.0.5 2021-04-16 [1] CRAN (R 4.3.0)
+## rtracklayer 1.60.0 2023-05-08 [1] Bioconductor
+## S4Arrays 1.0.4 2023-05-14 [1] Bioconductor
+## S4Vectors * 0.38.1 2023-05-08 [1] Bioconductor
+## sass 0.4.7 2023-07-15 [1] CRAN (R 4.3.0)
+## ScaledMatrix 1.8.1 2023-05-08 [1] Bioconductor
+## scales 1.2.1 2022-08-20 [1] CRAN (R 4.3.0)
+## scater * 1.28.0 2023-04-25 [1] Bioconductor
+## scran * 1.28.2 2023-07-23 [1] Bioconductor
+## scRNAseq * 2.14.0 2023-04-27 [1] Bioconductor
+## scuttle * 1.9.4 2023-01-23 [1] Bioconductor
+## sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.3.0)
+## shiny 1.7.4.1 2023-07-06 [1] CRAN (R 4.3.0)
+## SingleCellExperiment * 1.22.0 2023-05-08 [1] Bioconductor
+## sparseMatrixStats 1.12.2 2023-07-02 [1] Bioconductor
+## statmod 1.5.0 2023-01-06 [1] CRAN (R 4.3.0)
+## stringi 1.7.12 2023-01-11 [1] CRAN (R 4.3.0)
+## stringr 1.5.0 2022-12-02 [1] CRAN (R 4.3.0)
+## SummarizedExperiment * 1.30.2 2023-06-06 [1] Bioconductor
+## tibble 3.2.1 2023-03-20 [1] CRAN (R 4.3.0)
+## tidyselect 1.2.0 2022-10-10 [1] CRAN (R 4.3.0)
+## utf8 1.2.3 2023-01-31 [1] CRAN (R 4.3.0)
+## vctrs 0.6.3 2023-06-14 [1] CRAN (R 4.3.0)
+## vipor 0.4.5 2017-03-22 [1] CRAN (R 4.3.0)
+## viridis 0.6.4 2023-07-22 [1] CRAN (R 4.3.0)
+## viridisLite 0.4.2 2023-05-02 [1] CRAN (R 4.3.0)
+## withr 2.5.0 2022-03-03 [1] CRAN (R 4.3.0)
+## xfun 0.39 2023-04-20 [1] CRAN (R 4.3.0)
+## XML 3.99-0.14 2023-03-19 [1] CRAN (R 4.3.0)
+## xml2 1.3.5 2023-07-06 [1] CRAN (R 4.3.0)
+## xtable 1.8-4 2019-04-21 [1] CRAN (R 4.3.0)
+## XVector 0.40.0 2023-05-08 [1] Bioconductor
+## yaml 2.3.7 2023-01-23 [1] CRAN (R 4.3.0)
+## zlibbioc 1.46.0 2023-05-08 [1] Bioconductor
+##
+## [1] /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library
+##
+## ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
+
+Zheng, G. X. Y. et al. Massively parallel digital transcriptional profiling of single cells. Nat. Commun. 8, 14049 (2017).↩︎
Lun, A. T. L., Calero-Nieto, F. J., Haim-Vilmovsky, L., Göttgens, B. & Marioni, J. C. Assessing the reliability of spike-in normalization for analyses of single-cell RNA sequencing data. Genome Res. 27, 1795–1806 (2017)).↩︎