El perfilado del receptor de células T (TCR) y la Indexación Celular de Transcriptomas y Epítopos mediante Secuenciación (CITE-Seq) son técnicas fundamentales en la investigación de célula única, que ofrecen perspectivas sin precedentes sobre el sistema inmunitario adaptativo y la heterogeneidad celular. El perfilado de TCR permite profundizar en el repertorio y la diversidad de las poblaciones de células T, destacando la especificidad y singularidad de sus respuestas. Por otro lado, CITE-Seq facilita la evaluación simultánea de datos transcriptómicos y de expresión proteica en células individuales, generando una representación integral de los estados celulares. En este módulo, los participantes explorarán las profundas implicaciones del perfilado de TCR en la comprensión de las respuestas inmunitarias y las sinergias que puede lograr cuando se combina con CITE-Seq. Comenzaremos con conceptos y teorías fundamentales, para luego avanzar rápidamente hacia aplicaciones prácticas utilizando herramientas computacionales avanzadas. A través de este enfoque práctico, los asistentes dominarán las particularidades del perfilado de TCR y de CITE-Seq, adquiriendo herramientas valiosas para sus investigaciones inmunológicas y de célula única.
## Función para ejecutar comandos de terminal en Google Colab con el kernel de R
shell_call <- function(command, ...) {
result <- system(command, intern = TRUE, ...) # Ejecuta el comando de terminal y guarda la salida
cat(paste0(result, collapse = "\n")) # Imprime la salida en la consola
}
## Función para cargar múltiples paquetes de R
loadPackages = function(pkgs){
myrequire = function(...){
suppressWarnings(suppressMessages(suppressPackageStartupMessages(require(...))))
}
ok = sapply(pkgs, require, character.only=TRUE, quietly=TRUE) # Verifica si cada paquete está instalado
if (!all(ok)){
message("Faltan los siguientes paquetes: ", paste(pkgs[!ok], collapse=", ")) # Muestra los paquetes faltantes
}
}
## Descargar el código para agregar R2U (sistema rápido de instalación de paquetes)
download.file("https://github.com/eddelbuettel/r2u/raw/master/inst/scripts/add_cranapt_jammy.sh",
"add_cranapt_jammy.sh")
Sys.chmod("add_cranapt_jammy.sh", "0755") # Cambiar permisos para permitir ejecución
## Ejecutar el código para configurar R2U
shell_call("./add_cranapt_jammy.sh")
## Habilitar el gestor de paquetes BSPM para instalar paquetes a nivel del sistema
bspm::enable()
options(bspm.version.check=FALSE)
## Eliminar el código de configuración para mantener limpio el entorno de trabajo
shell_call("rm add_cranapt_jammy.sh")
## Definir la lista de paquetes a instalar
cranPkgs2Install = c("dplyr", "ggpubr", "Seurat", "cowplot",
"Rtsne", "hdf5r", "patchwork")
## Instalar todos los paquetes sin pedir confirmación al usuario
install.packages(cranPkgs2Install, ask=FALSE, update=TRUE, quietly=TRUE)
## Para simplificar la carga de paquetes, usamos la función loadPackages().
## Pero si no tienes la función, puedes usar 'library(nombre_del_paquete)'
pkgs = c("Seurat", "dplyr", "patchwork") # Lista de paquetes a cargar
loadPackages(pkgs) # Carga los paquetes usando la función definida previamente
## Definir el directorio donde se guardarán y accederán los datos
## IMPORTANTE: El usuario debe cambiar "scw01" por el nombre real de su directorio de trabajo
mydir <- "/content"
En este cuaderno (notebook), profundizaremos en el análisis de datos unicelulares multimodales. Nos centraremos en una muestra de cáncer de pulmón de células no pequeñas (CPNM) procesada mediante la tecnología de perfil inmunológico 10X 5'. Esta tecnología avanzada captura las secuencias de ARN y de receptores de células T de cada célula, lo que proporciona una visión completa del panorama celular..
Puede descargar el conjunto de datos desde el sitio web de 10x Genomics aquí.
El enfoque de secuenciación de TCR en célula única permite un análisis de alta resolución de la diversidad de los receptores de células T, identificando clonotipos específicos y sus funciones en el sistema inmunológico. A diferencia de la secuenciación masiva, que combina información de múltiples células, el método de célula única preserva la identidad individual de cada célula T, permitiendo la correlación entre el perfil transcriptómico y el repertorio de TCR.
En este cuaderno, repasaremos brevemente los pasos clave del análisis de célula única. Dado que este módulo se centra en el análisis de la secuencia de TCR, proporcionaremos una visión general de los principales procesos. Para explicaciones más detalladas de cada paso, consulte los otros módulos disponibles.
# Descargar una matriz de genes y códigos de barras filtrada de 10X Genomics
# Este comando usa curl para descargar un archivo desde internet. (-O) Guarda el archivo con el mismo nombre original.
shell_call("curl -O https://cf.10xgenomics.com/samples/cell-vdj/5.0.0/vdj_v1_hs_nsclc_multi_5gex_t_b/vdj_v1_hs_nsclc_multi_5gex_t_b_count_filtered_feature_bc_matrix.tar.gz")
## Extraer el archivo descargado (descomprimir el conjunto de datos)
# tar -xf le indica a tar que extraiga (-x) el archivo especificado (-f) y lo descomprima.
shell_call("tar -xf /content/vdj_v1_hs_nsclc_multi_5gex_t_b_count_filtered_feature_bc_matrix.tar.gz")
# Leer la matriz filtrada de características y códigos de barras en un objeto tipo matriz
# Este comando carga los datos de un archivo 10X Genomics
# ubicado en el directorio especificado y los asigna a la variable counts.
counts <- Read10X(paste0(mydir, "/filtered_feature_bc_matrix/"))
A continuación, utilizaremos la matriz de cuentas y los metadatos para crear un objeto Seurat. Este objeto actúa como un contenedor integral que alberga no solo los datos de cuentas sin procesar y los metadatos, sino también los resultados de análisis posteriores, como el PCA (Análisis de Componentes Principales) y los resultados de agrupamiento (clustering). Este contenedor centralizado permite una manipulación y un análisis de datos optimizados, proporcionando una estructura organizada para nuestro conjunto de datos scRNA-Seq.
Al crear un objeto Seurat, podemos gestionar y analizar nuestros datos de forma eficiente, lo que facilita la realización de operaciones complejas y la visualización de los resultados.
# Crear un objeto Seurat usando la matriz de cuentas crudas
seurat.raw <- CreateSeuratObject(counts = counts)
# Mostrar el contenido del objeto Seurat
seurat.raw
El primer paso de nuestro análisis es familiarizarnos con el conjunto de datos almacenado en el objeto Seurat. Seurat proporciona potentes herramientas que nos permiten explorar y visualizar nuestros datos de cuentas junto con los metadatos asociados. Esto nos permite comprender a fondo nuestros datos de scRNA-Seq.
Al aprovechar las capacidades de Seurat, podemos realizar diversos análisis exploratorios, como:
Al realizar estas exploraciones iniciales, sentamos las bases para análisis más avanzados, como la reducción de dimensionalidad, el agrupamiento (clustering) y el análisis de expresión diferencial. Este paso fundamental es crucial para comprender las particularidades de nuestro conjunto de datos y tomar decisiones informadas durante todo el proceso de análisis.
# ¿Cuántas células y genes tenemos actualmente?
print(paste0("El número de genes es ", dim(seurat.raw)[1], " y el número de células es ", dim(seurat.raw)[2]))
# La función print() muestra el resultado en la consola.
# La función dim() devuelve las dimensiones de un objeto, generalmente una matriz o un data frame.
# En los objetos Seurat, la matriz de cuentas tiene genes como filas y células como columnas.
# paste0() concatena (une) cadenas de texto sin espacios.
# Las frases dentro de los paréntesis (") también se imprimirán.
# Ver una porción (slice) de la matriz de cuentas
# Recuerda: filas = genes, columnas = células/códigos de barras
GetAssayData(seurat.raw, slot = "counts")[8:10,13:14]
NOTA: Los puntos ('.') indican un valor cero. La tabla de conteos se almacena en formato de matriz dispersa, que almacena explícitamente sólo valores distintos de cero para ahorrar espacio.
## Mostrar las columnas de metadatos disponibles en el objeto Seurat
# ¿Qué columnas de metadatos hay disponibles?
print(colnames(seurat.raw@meta.data))
# (@) Accede al slot meta.data del objeto seurat.raw.
# El slot meta.data usualmente contiene metadatos de las células como:
# el tipo celular, ID de muestra u otras anotaciones.
# colnames() recupera los nombres de columna del data frame almacenado en meta.data
# Crear un gráfico de violín que muestre la distribución del número de UMIs por célula
options(repr.plot.width=7, repr.plot.height=7) # Define el tamaño del gráfico
VlnPlot(seurat.raw, features = c("nCount_RNA"),y.max=2e4)
# Crea un gráfico de violín para el número de moléculas de ARN (nCount_RNA) por célula
# y.max=2e4: establece el valor máximo del eje y a 20,000
Ya hemos observado que Seurat calculó automáticamente dos de las principales métricas de control de calidad (CC): el número de UMI (Identificadores Moleculares Únicos) y el número de genes detectados por célula. La siguiente métrica importante de CC a considerar es el porcentaje de genes mitocondriales. Los genes mitocondriales pueden indicar estrés celular o apoptosis, por lo que es crucial monitorizar sus niveles de expresión.
Para calcular el porcentaje de genes mitocondriales, utilizaremos el método PercentageFeatureSet de Seurat. Este método calcula el porcentaje de UMI que se originan a partir de genes que coinciden con un patrón específico. Para los genes mitocondriales, normalmente buscamos genes que empiezan por "MT-" (prefijo común para genes mitocondriales).
# Los genes mitocondriales humanos comienzan con "MT-",
# así que calculamos el porcentaje que coincide con el patrón "^MT-"
# Esta función calcula el porcentaje de cuentas para los genes que coinciden con ese patrón.
# También crea una nueva columna de metadatos llamada percent.mt con esos valores.
seurat.raw[["percent.mt"]] <- PercentageFeatureSet(seurat.raw, pattern = "^MT-")
# Ahora podemos ver que el % de expresión mitocondrial ha sido calculado para cada célula
head(seurat.raw$percent.mt) # Muestra las primeras filas del porcentaje mitocondrial
# Gráfico de violín para tres métricas de calidad: UMIs, número de genes, % mitocondrial
options(repr.plot.width=12, repr.plot.height=6)
VlnPlot(seurat.raw, features = c("nCount_RNA", "nFeature_RNA", "percent.mt"))
# Podemos visualizar las tres métricas de calidad celular juntas utilizando el método VlnPlot de Seurat
Suele ser útil visualizar estas métricas de control de calidad en conjunto, ya que las células con valores atípicos (outliers) en múltiples dimensiones tienen mayor probabilidad de ser de baja calidad. FeatureScatter de Seurat crea un diagrama de dispersión de dos columnas dadas a partir de nuestros metadatos.
# Gráfico de dispersión: UMIs totales vs porcentaje de genes mitocondriales
options(repr.plot.width=6, repr.plot.height=6)
FeatureScatter(seurat.raw, feature1 = "nCount_RNA", feature2 = "percent.mt")
# Visualiza UMIs vs % de genes mitocondriales
# Una vez visualizadas las métricas, podemos definir los umbrales para filtrar
# Usamos la función subset de R
seurat.raw <- subset(
seurat.raw,
subset =
nFeature_RNA > 200 & # Elimina células con menos de 200 genes
nCount_RNA > 400 & # Elimina células con menos de 400 UMIs
nFeature_RNA < 6000 & # Elimina células con más de 6000 genes (posibles dobles)
percent.mt < 40) # Elimina células con más del 40% de ARN mitocondrial
Tras eliminar las células no deseadas del conjunto de datos, el siguiente paso es normalizar los datos. La normalización es esencial para ajustar las diferencias en la profundidad de secuenciación entre las células, garantizando así comparaciones significativas.
1.Método de Normalización: LogNormalize
2.Normalización por Expresión Total:
3.Multiplicación por un Factor de Escala:
4.Transformación Logarítmica:
5.Almacenamiento de datos:
# Normalizar los datos usando LogNormalization
seurat.raw <- NormalizeData(seurat.raw,
normalization.method = "LogNormalize", scale.factor = 10000)
# Identificar los genes más variables
seurat.raw <- FindVariableFeatures(seurat.raw, selection.method = "vst", nfeatures = 2000)
# Escalar los datos
seurat.raw <- ScaleData(seurat.raw, features = VariableFeatures(seurat.raw), do.scale = T, do.center = T)
# Ejecutar PCA (Análisis de Componentes Principales)
seurat.raw <- RunPCA(seurat.raw, features = VariableFeatures(seurat.raw))
# Encontrar vecinos celulares
seurat.raw <- FindNeighbors(seurat.raw, dims = 1:20, k.param = 20)
# Ejecutar UMAP con reducción PCA
seurat.raw <- RunUMAP(seurat.raw, dims = 1:20, reduction = "pca", seed.use = 1)
# Encontrar clústeres con baja resolución (0.1 es un buen punto de partida para identificar todas las células T)
seurat.raw <- FindClusters(seurat.raw, resolution = 0.1)
## Visualizar UMAP con etiquetas de la agrupación (clúster)
DimPlot(seurat.raw, reduction = "umap", label = T,group.by = "seurat_clusters")
Observe el gráfico de características de los marcadores de tipo celular común para determinar qué grupo corresponde a las células T.
## Mapa de expresión de genes marcadores específicos para identificar tipos celulares
FeaturePlot(seurat.raw, features = c("CD3D", "CD8A", "GNLY",
"CD79A", "FCGR3A", "KRT7"), min.cutoff = "q1")
En esta sección leeremos archivos generados por el software 10X CellRanger que describen las secuencias TCR encontradas en células individuales.
## Descargar archivos con curl
shell_call("curl -O https://cf.10xgenomics.com/samples/cell-vdj/5.0.0/vdj_v1_hs_nsclc_multi_5gex_t_b/vdj_v1_hs_nsclc_multi_5gex_t_b_vdj_t_all_contig_annotations.csv")
shell_call("curl -O https://cf.10xgenomics.com/samples/cell-vdj/5.0.0/vdj_v1_hs_nsclc_multi_5gex_t_b/vdj_v1_hs_nsclc_multi_5gex_t_b_vdj_t_clonotypes.csv")
## Leer archivo con la información de TCR por célula
tcr <- read.csv(paste0(mydir, "/vdj_v1_hs_nsclc_multi_5gex_t_b_vdj_t_all_contig_annotations.csv"))
## Leer archivo con la información de clonotipos (varias células pueden compartir el mismo clonotipo o secuencia TCR)
clono <- read.csv(paste0(mydir, "/vdj_v1_hs_nsclc_multi_5gex_t_b_vdj_t_clonotypes.csv"))
## Eliminar el sufijo "-1" de cada código de barra (barcode)
tcr$barcode <- gsub("-1", "", tcr$barcode)
## Filtrar para mantener solo la primera línea de cada barcode, así analizamos solo un clonotipo por célula
tcr <- tcr[!duplicated(tcr$barcode), ]
## También eliminamos el sufijo "-1" en los nombres de las células del objeto Seurat
seurat.raw <- RenameCells(seurat.raw, new.names = gsub("-1", "", colnames(seurat.raw)))
## Mantener solo las columnas de código de barras (barcode) y clonotipo. Renombrar "raw_clonotype_id" a "clonotype_id"
tcr <- tcr[, c("barcode", "raw_clonotype_id")]
names(tcr)[names(tcr) == "raw_clonotype_id"] <- "clonotype_id"
## Unir las tablas de TCR y clonotipos para obtener la secuencia de aminoácidos del TCR por célula
tcr <- merge(tcr, clono[, c("clonotype_id", "cdr3s_aa")])
## Reordenar columnas para tener código de barras (barcode) primero, establecerlos como nombre de filas (rownames) y eliminar columna extra de códigos de barras (barcodes)
tcr <- tcr[, c(2, 1, 3)]
rownames(tcr) <- tcr[, 1]
tcr[, 1] <- NULL
## Agregar la información de TCR como metadatos al objeto Seurat
seurat.raw <- AddMetaData(object = seurat.raw, metadata = tcr)
## Verificar que la información de TCR se haya agregado correctamente
head(seurat.raw@meta.data)
## Crear un histograma de frecuencias de clones a partir de la tabla clono
## ¿Cuál sería un buen umbral para distinguir clones expandidos de los que no lo están?
barplot(table(clono$frequency), xlab = "Número de células por clon")
# La función barplot() crea un gráfico de barras.
# La función table() crea una tabla de contingencia de los recuentos en cada valor único de clono$frequency.
# clono$frequency accede a la columna de frecuencia en el marco de datos clono.
## Identificar clones expandidos, marcarlos en los metadatos y etiquetarlos en el UMAP
## ¿Cuántos clones expandidos hay? Prueba cambiar el umbral y comparar los resultados
Nexpand = 1 # Este es el umbral para considerar un clon como expandido
length(which(clono$frequency > 1))
# which(clono$frequency > 1): Devuelve los índices de las filas donde la frecuencia es mayor que 1.
# length(which(clono$frequency > 1)): Cuenta el número de elementos (clones) que tienen una frecuencia mayor que 1. Esto proporciona el número de clones expandidos.
expanded_clones = clono$clonotype_id[clono$frequency>1]
# clono$clonotype_id[clono$frequency > 1]: Selecciona los valores de clonotype_id cuya frecuencia correspondiente es mayor que 1.
# expanded_clones: Almacena los ID de clonotipo de los clones expandidos en una nueva variable.
# Agregue una nueva columna de metadatos que indique qué células son parte de un clon expandido
seurat.raw = AddMetaData(seurat.raw,metadata = rep("no",ncol(seurat.raw)),col.name="TCR_expanded")
# Este comando añade una nueva columna de metadatos al objeto Seurat seurat.raw
# en estos nuevos metadatos, cada célula tiene un valor "no" y luego se cambia el valor de las células expandidas a "sí".
seurat.raw@meta.data$TCR_expanded[seurat.raw@meta.data$clonotype_id %in% expanded_clones] = "yes"
#This command updates the TCR_expanded metadata column.
# seurat.raw@meta.data$TCR_expanded: Accede a la columna TCR_expanded en los metadatos.
# seurat.raw@meta.data$clonotype_id %in% expanded_clones: Comprueba qué celdas tienen valores de clonotype_id en la lista expanded_clones.
# = "yes": Establece el valor "yes" para las células que pertenecen a clones expandidos.
# Mire el UMAP para ver si los clones expandidos tienen una expresión genética similar
DimPlot(seurat.raw, reduction = "umap", group.by = "TCR_expanded", label = T)
# ¿Las células T expandidas expresan genes diferentes a los de las células T no expandidas?
# Para responder a esta pregunta, primero separamos las células T del resto del objeto.
# Revisa tu agrupamiento para ver cuál parece contener células T.
Idents(seurat.raw) = "seurat_clusters" # Esto significa que las operaciones posteriores utilizarán las identidades de agrupamiento (clúster)
# asignadas a cada célula para seurat_clusters
seurat.t = subset(seurat.raw, idents = "1") #es el nuevo objeto Seurat que contiene solo las células del grupo 1.
seurat.t # Verifique el número de muestras (células) para ver con cuántas células T tenemos que trabajar
Idents(seurat.t) = "TCR_expanded" # Este comando establece la clase de identidad activa en el objeto seurat.t en TCR_expanded
deg_expanded = FindMarkers(seurat.t,ident.1="yes",ident.2="no",logfc.threshold = 0.25,min.pct = 0.1)
# Este comando identifica genes con expresión diferencial entre dos grupos de células.
# Los marcados como "sí" en los metadatos TCR_expanded y los marcados como "no".
# logfc.threshold: Umbral mínimo de cambio log2 para identificar genes con expresión diferencial.
# min.pct: Porcentaje mínimo de células en las que se detecta el gen.
# Visualice la expresión de los principales genes DE para cada subconjunto de células anotadas.
top30genes <- deg_expanded %>% filter(avg_log2FC > 0) %>% top_n(30, avg_log2FC)
# deg_expanded %>% filter(avg_log2FC > 0): Filtra los genes con expresión diferencial (deg_expanded)
# para incluir solo aquellos con un cambio de log2 promedio positivo (avg_log2FC > 0).
# top_n(30, avg_log2FC): De los genes filtrados, selecciona los 30 genes principales con el mayor cambio de log2 promedio.
genes <- rownames(top30genes)
seurat.t <- ScaleData(seurat.t, features = genes, do.center = T, do.scale = T)
# Escala los datos de expresión de los genes seleccionados (features = genes) en el objeto seurat.t.
# do.center = T: Centra los datos restando la expresión media de cada gen.
# do.scale = T: Escala los datos dividiendo por la desviación estándar de cada gen.
DoHeatmap(seurat.t, features = genes)
# Crea un mapa de calor de los datos de expresión