参考：https://mp.weixin.qq.com/s/x4HOHNwikvuTCRK4kQRjRA?token=971255988&lang=zh_CN
一般而言，R分析单细胞使用Seurat，python分析单细胞使用Scanpy，都是很好得工作。可是有些时候，我们希望两者之间进行转化，或者更多的情况是可以自由切换进行数据分析。关于seurat转h5ad我们之前已经讲过了不再赘述，这里我们介绍4种途径（4种包）将h5ad转化为seurat。这样做的目的是让大家的选择性更高，一种不行立马换，不用因为解决报错耽误时间！因为我们没有scanpy构建的数据，以及考虑到一些包的更新，方便转化过程中一些error的解决，所以我们按照官网流程走了一遍scanpy分析，流程在文后，开头先上重点内容吧！
接下来我们看看具体的方法：
ps：推荐指数仅仅关乎使用简易程度和出bug频率，无他！
一、SeuratDisk: 推荐指数⭐⭐⭐⭐
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 <- meta
rownames(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)
二、anndataR: 推荐指数⭐⭐⭐
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_umap
colnames(a) <- c("umap_1","umap_2")
rownames(a) <- adata$obs_names
A <- 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)
三、schard: 推荐指数⭐⭐⭐⭐⭐
#https://github.com/cellgeni/schard
devtools::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)
四、sceasy: 推荐指数⭐⭐
#如果可以就用，不行就放弃使用别的方法
# 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 pd
import scanpy as sc
sc.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 natrix
adata = 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 QC
sc.pp.filter_cells(adata, min_genes=200)#each cell expres 200 genes at least
sc.pp.filter_genes(adata, min_cells=3) #each gene must expres in 3cells at least
# mt gene percentage
adata.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 = adata
adata = 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')
