En este tutorial, construimos un proceso completo para el análisis de secuenciación de ARN unicelular utilizando Scanpy. Comenzamos instalando las bibliotecas necesarias y cargando el conjunto de datos PBMC 3k, luego realizamos control de calidad, filtrado y normalización para preparar los datos para el análisis posterior. Luego identificamos genes altamente variables, realizamos PCA para reducir la dimensionalidad y construimos un gráfico de vecindad para generar incrustaciones UMAP y grupos de Leiden. A través del descubrimiento y visualización de genes marcadores, exploramos cómo los grupos se corresponden con poblaciones de células biológicas e implementamos una estrategia de anotación simple basada en reglas para inferir tipos de células.
importar sys importar subproceso importar importlib def pip_install(*paquetes): subproceso.check_call([sys.executable, “-m”, “pip”, “install”, “-q”, *packages]) requerido = [
“scanpy”,
“anndata”,
“leidenalg”,
“igraph”,
“harmonypy”,
“seaborn”
]
pip_install(*required) importar os importar advertencias advertencias.filterwarnings(“ignorar”) importar numpy como np importar pandas como pd importar matplotlib.pyplot como plt importar scanpy como sc importar anndata como ad sc.settings.verbosity = 2 sc.settings.set_figure_params(dpi=110, facecolor=”white”, frameon=False) np.random.seed(42) print(“Versión Scanpy:”, sc.__version__) adata = sc.datasets.pbmc3k() adata.var_names_make_unique() print(“\nInitial AnnData:”) print(adata) adata.layers[“counts”] = adata.X.copia() adata.var[“mt”] = adata.var_names.str.upper().str.startswith(“MT-“) sc.pp.calculate_qc_metrics(adata, qc_vars=[“mt”]percent_top=Ninguno, log1p=False, inplace=True) print(“\nResumen de control de calidad:”) display( adata.obs[[“n_genes_by_counts”, “total_counts”, “pct_counts_mt”]].describe().T )
Instalamos todas las dependencias necesarias e importamos las bibliotecas informáticas científicas básicas necesarias para el análisis. Configuramos los ajustes de Scanpy, inicializamos el entorno y cargamos el conjunto de datos RNA-seq unicelular PBMC 3k. Luego calculamos métricas de control de calidad, incluido el porcentaje de genes mitocondriales, los recuentos totales y la cantidad de genes detectados, para cada célula.
higo, ejes = plt.subplots(1, 3, tamaño de higo=(15, 4)) sc.pl.violin(adata, [“n_genes_by_counts”]fluctuación = 0,4, ax = ejes[0]mostrar=Falso) sc.pl.violin(adata, [“total_counts”]fluctuación = 0,4, ax = ejes[1]mostrar=Falso) sc.pl.violin(adata, [“pct_counts_mt”]fluctuación = 0,4, ax = ejes[2]show=False) plt.tight_layout() plt.show() sc.pl.scatter(adata, x=”total_counts”, y=”n_genes_by_counts”, color=”pct_counts_mt”) adata = adata[adata.obs[“n_genes_by_counts”] >= 200].copiar() adata = adata[adata.obs[“n_genes_by_counts”] <= 5000].copiar() adata = adata[adata.obs["pct_counts_mt"] < 10].copy() sc.pp.filter_genes(adata, min_cells=3) print("\nDespués del filtrado:") print(adata) sc.pp.normalize_total(adata, target_sum=1e4) sc.pp.log1p(adata) adata.raw = adata.copy() sc.pp.highly_variable_genes( adata, sabor="seurat", min_mean=0.0125, max_mean=3, min_disp=0.5 ) print("\nGenes altamente variables seleccionados:", int(adata.var["highly_variable"].sum())) sc.pl.highly_variable_genes(adata) adata = adata[:, adata.var["highly_variable"]].Copiar()
Visualizamos métricas de control de calidad utilizando gráficos para verificar la distribución de los recuentos de genes y el contenido mitocondrial. Aplicamos pasos de filtrado para eliminar células y genes de baja calidad que no cumplen con los umbrales de expresión básicos. Luego normalizamos los datos, aplicamos una transformación logarítmica e identificamos genes altamente variables para el análisis posterior.
sc.pp.regress_out(adatos, [“total_counts”, “pct_counts_mt”]) sc.pp.scale(adata, max_value=10) sc.tl.pca(adata, svd_solver=”arpack”) sc.pl.pca_variance_ratio(adata, log=True) sc.pl.pca(adata, color=Ninguno) sc.pp.neighbors(adata, n_neighbors=12, n_pcs=30, metric=”euclidean”) sc.tl.umap(adata, min_dist=0.35, spread=1.0) sc.tl.leiden(adata, resolución=0.6, key_added=”leiden”) print(“\nRecuentos de clústeres:”) display(adata.obs[“leiden”].value_counts().sort_index().rename(“cells_per_cluster”).to_frame()) sc.pl.umap(adata, color=[“leiden”]legend_loc=”on data”, title=”PBMC 3k – Clústeres de Leiden”) sc.tl.rank_genes_groups(adata, groupby=”leiden”, método=”wilcoxon”) sc.pl.rank_genes_groups(adata, n_genes=20, sharey=False) marcador_table = sc.get.rank_genes_groups_df(adata, group=None) print(“\nTop filas de marcadores:”) display(marker_table.head(20))
Eliminamos los factores de confusión técnicos y escalamos el conjunto de datos para prepararlo para la reducción de dimensionalidad. Realizamos análisis de componentes principales para capturar la estructura de varianza principal del conjunto de datos. Luego construimos el gráfico de vecindad, calculamos incrustaciones de UMAP, realizamos agrupaciones de Leiden e identificamos genes marcadores para cada grupo.
top_markers_per_cluster = ( marcador_table.groupby(“grupo”) .head(10) .loc[:, [“group”, “names”, “logfoldchanges”, “pvals_adj”]].reset_index(drop=True) ) print(“\nTop 10 marcadores por grupo:”) display(top_markers_per_cluster) candidato_marcadores = [
“IL7R”, “LTB”, “MALAT1”, “CCR7”,
“NKG7”, “GNLY”, “PRF1”,
“MS4A1”, “CD79A”, “CD79B”,
“LYZ”, “S100A8”, “FCER1A”, “CST3”,
“PPBP”, “FCGR3A”, “LGALS3”, “CTSS”,
“CD3D”, “TRBC1”, “TRAC”
]
marcadores_candidatos = [g for g in candidate_markers if g in adata.var_names]
if marcadores_candidatos: sc.pl.dotplot( adata, var_names=marcadores_candidatos, groupby=”leiden”, standard_scale=”var”, dendrogram=True ) sc.pl.matrixplot( adata, var_names=marcadores_candidatos, groupby=”leiden”, standard_scale=”var”, dendrogram=True ) cluster_marker_reference = { “T_cells”: [“IL7R”, “LTB”, “CCR7”, “CD3D”, “TRBC1”, “TRAC”]”NK_celdas”: [“NKG7”, “GNLY”, “PRF1″]”células_B”: [“MS4A1”, “CD79A”, “CD79B”]”Monocitos”: [“LYZ”, “FCGR3A”, “LGALS3”, “CTSS”, “S100A8”, “CST3″]”células_dendriticas”: [“FCER1A”, “CST3″]”Plaquetas”: [“PPBP”]
}
Examinamos los genes marcadores más importantes detectados para cada grupo y resumimos los marcadores principales. Visualizamos patrones de expresión genética en grupos utilizando diagramas de puntos y diagramas de matriz para marcadores de células inmunes conocidos. También definimos un mapeo de referencia de genes marcadores asociados con los principales tipos de células inmunes para su posterior anotación.
referencia_disponible = {tipo de celda: [g for g in genes if g in adata.var_names]
para tipo de célula, genes en cluster_marker_reference.items() } disponible_referencia = {k: v para k, v en disponible_reference.items() si len(v) > 0} para tipo de célula, genes en disponible_reference.items(): sc.tl.score_genes(adata, gene_list=genes, score_name=f”{celltype}_score”, use_raw=False) score_cols = [f”{ct}_score” for ct in available_reference.keys()]
cluster_scores = adata.obs.groupby(“leiden”)[score_cols].mean() display(cluster_scores) cluster_to_celltype = {} para el clúster en cluster_scores.index: mejor = cluster_scores.loc[cluster].idxmax().replace(“_score”, “”) cluster_to_celltype[cluster] = mejor adata.obs[“cell_type”] =adata.obs[“leiden”].map(cluster_to_celltype).astype(“categoría”) print(“\nAsignación de tipo de celda a grupo:”) display(pd.DataFrame.from_dict(cluster_to_celltype, orient=”index”, columns=[“assigned_cell_type”])) sc.pl.umap(adata,color=[“leiden”, “cell_type”]legend_loc=”on data”, wspace=0.45 ) sc.tl.rank_genes_groups(adata, groupby=”cell_type”, método=”wilcoxon”) sc.pl.rank_genes_groups(adata, n_genes=15, sharey=False) celltype_markers = sc.get.rank_genes_groups_df(adata, group=None) print(“\nMarcadores principales por tipo de celda anotado:”) display( celltype_markers.groupby(“group”).head(8)[[“group”, “names”, “logfoldchanges”, “pvals_adj”]]) cluster_prop = (adata.obs[“cell_type”]
.value_counts(normalize=True) .mul(100) .round(2) .rename(“porcentaje”) .to_frame() ) print(“\nProporciones de tipo de celda (%):”) display(cluster_prop) plt.figure(figsize=(7, 4)) cluster_prop[“percent”].sort_values().plot(kind=”barh”) plt.xlabel(“Porcentaje de celdas”) plt.ylabel(“Tipo de celda”) plt.title(“Composición estimada del tipo de celda”) plt.tight_layout() plt.show() output_dir = “scanpy_pbmc3k_outputs” os.makedirs(output_dir, exist_ok=True) adata.write(os.path.join(output_dir, “pbmc3k_scanpy_advanced.h5ad”)) marcador_table.to_csv(os.path.join(output_dir, “cluster_markers.csv”), index=False) celltype_markers.to_csv(os.path.join(output_dir, “celltype_markers.csv”), index=False) cluster_scores.to_csv(os.path.join(output_dir, “cluster_score_matrix.csv”)) print(f”\nSalidas guardadas en: {output_dir}”) print(“Archivos:”) for f in sorted(os.listdir(output_dir)): print(” -“, f) resumen = { “n_cells_final”: int(adata.n_obs), “n_genes_final”: int(adata.n_vars), “n_clusters”: int(adata.obs[“leiden”].nunique()), “clústeres”: ordenados(adata.obs[“leiden”].unique().tolist()), “tipos_celdas”: ordenados(adata.obs[“cell_type”].unique().tolist()), } print(“\nResumen del análisis:”) para k, v en resumen.items(): print(f”{k}: {v}”)
Calificamos cada célula utilizando conjuntos de genes marcadores conocidos y asignamos tipos de células probables a grupos según patrones de expresión. Visualizamos los tipos de células anotados en la incrustación de UMAP y realizamos análisis diferenciales de expresión génica en las poblaciones de células previstas. Además, calculamos proporciones de tipo de celda, generamos visualizaciones resumidas y guardamos el conjunto de datos procesados y los resultados del análisis para futuras investigaciones.
En conclusión, desarrollamos un flujo de trabajo completo de un extremo a otro para analizar datos transcriptómicos unicelulares utilizando Scanpy. Realizamos preprocesamiento, agrupamiento, análisis de genes marcadores y anotaciones de tipo celular, y visualizamos la estructura de datos utilizando UMAP y gráficos de expresión génica. Al guardar el objeto AnnData procesado y los resultados del análisis, creamos un conjunto de datos reutilizable para una mayor interpretación biológica y modelado avanzado. Este flujo de trabajo demuestra cómo Scanpy permite un análisis unicelular escalable y reproducible a través de un canal de Python estructurado y modular.
Consulte los códigos completos aquí. Además, no dude en seguirnos en Twitter y no olvide unirse a nuestro SubReddit de más de 120.000 ML y suscribirse a nuestro boletín. ¡Esperar! estas en telegrama? Ahora también puedes unirte a nosotros en Telegram.
La publicación Una guía de codificación para construir un canal completo de análisis de secuenciación de ARN unicelular utilizando Scanpy para visualización de agrupaciones y anotación de tipos de células apareció por primera vez en MarkTechPost.