玩转单细胞(16)：Scanpy单细胞h5ad数据转化为Seurat对象

getwd()setwd("/home/ks_ts/data_analysis/Scanpy分析/")

library(SeuratDisk)Convert('./pbmc3k.h5ad', dest = "h5seurat", overwrite = TRUE)seuratObj <- LoadH5Seurat("./pbmc3k.h5seurat", meta.data = FALSE, misc = FALSE)
meta <- read.csv('./meta.csv',header = T)seuratObj@meta.data <- metarownames(seuratObj@meta.data) <- meta$X
Idents(seuratObj) <- 'leiden'DimPlot(seuratObj)
df_markers <- FindAllMarkers(seuratObj, logfc.threshold = 0.2, min.pct = 0.2, only.pos = T)

devtools::install_github("scverse/anndataR")library(anndataR)# adata.raw.to_adata().write("./pbmc3k_withoutX.h5ad")adata <- anndataR::read_h5ad("./pbmc3k_withoutX.h5ad", to = "HDF5AnnData")seuratObj2 <- adata$to_Seurat()

a = adata$obsm$X_umapcolnames(a) <- c("umap_1","umap_2")rownames(a) <- adata$obs_namesA <- Seurat::CreateDimReducObject(as.matrix(a), key="umap",assay = 'RNA')
seuratObj2@reductions[['umap']] <- A

Idents(seuratObj2) <- 'leiden'DimPlot(seuratObj2)
# df_markers <- FindAllMarkers(seuratObj2, logfc.threshold = 0.2, min.pct = 0.2, only.pos = T)

#https://github.com/cellgeni/scharddevtools::install_github("cellgeni/schard")
seuratObj3 = schard::h5ad2seurat('pbmc3k.h5ad',use.raw=T)
Idents(seuratObj3) <- 'leiden'DimPlot(seuratObj3)df_markers <- FindAllMarkers(seuratObj3, logfc.threshold = 0.2, min.pct = 0.2, only.pos = T)

#如果可以就用，不行就放弃使用别的方法# devtools::install_github("cellgeni/sceasy")# library(sceasy)# library(reticulate)# use_condaenv('sceasy')# loompy <- reticulate::import('loompy')# # ad_path <- "./pbmc3k.h5ad"# sceasy::convertFormat(ad_path, from="anndata", to="seurat", outFile="file.rds")

pip install scanpy -i https://pypi.tuna.tsinghua.edu.cn/simple

import pandas as pdimport scanpy as scsc.settings.verbosity = 3  # verbosity: errors (0), warnings (1), info (2), hints (3)sc.logging.print_header()sc.settings.set_figure_params(dpi=80, facecolor="white")#scanpy==1.9.8 anndata==0.9.2 umap==0.5.7 numpy==1.24.4 scipy==1.10.1 pandas==2.0.3 scikit-learn==1.3.2 statsmodels==0.14.1 pynndescent==0.5.13

#read natrixadata = sc.read_10x_mtx(    "./filtered_gene_bc_matrices/hg19/",  # the directory with the `.mtx` file    var_names="gene_symbols",  # use gene symbols for the variable names (variables-axis index)    cache=True,  # write a cache file for faster subsequent reading)

#basic QCsc.pp.filter_cells(adata, min_genes=200)#each cell expres 200 genes at leastsc.pp.filter_genes(adata, min_cells=3) #each gene must expres in 3cells at least

# mt gene percentageadata.var["mt"] = adata.var_names.str.startswith("MT-")sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], percent_top=None, log1p=False, inplace=True)adata

adata = adata[adata.obs.n_genes_by_counts < 2500, :]adata = adata[adata.obs.pct_counts_mt < 5, :].copy()

sc.pp.normalize_total(adata, target_sum=1e4)

sc.pp.log1p(adata)sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)

adata.raw = adataadata = adata[:, adata.var.highly_variable]

#Regress out effects of total counts per cell and the percentage of mitochondrial genes expressed. sc.pp.regress_out(adata, ["total_counts", "pct_counts_mt"])

sc.pp.scale(adata, max_value=10)sc.tl.pca(adata, svd_solver="arpack")

sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)

sc.tl.umap(adata)

sc.pl.umap(adata, color='leiden', legend_loc='on data', title='', frameon=False)

new_cluster_names = [    'CD4 T', 'CD14 Monocytes',    'B', 'CD8 T',     'NK', 'FCGR3A Monocytes',    'Dendritic', 'Megakaryocytes']adata.rename_categories('leiden', new_cluster_names)

marker_genes = ['IL7R', 'CD79A', 'MS4A1','CD8A', 'CD8B', 'LYZ', 'CD14',                'LGALS3', 'S100A8', 'GNLY', 'NKG7', 'KLRB1',                  'FCGR3A', 'MS4A7', 'FCER1A', 'CST3', 'PPBP']sc.pl.dotplot(adata, marker_genes, groupby='leiden')

adata.write('./pbmc3k.h5ad')adata.obs.to_csv('./meta.csv')