Seurat-Scanpy互通 | R-Python互通
Tips:
Atlas级别的数据，需要做加速处理，比如常见的marker鉴定，不要用seurat默认的函数，非常慢。
presto包的算法复杂度极低，再大的数据基本一分钟内就能出结果。
markers_rna <- presto:::wilcoxauc.Seurat(X = all.colon, group_by = 'Major', assay = 'data', seurat_assay = 'integrated') # markers_motifs <- presto:::wilcoxauc.Seurat(X = seuset.flt, group_by = 'celltype', assay = 'data', seurat_assay = 'chromvar')
单细胞数据分析发展到现在，已经分化为两大阵营，Seurat和Scanpy，一个代表了R，有丰富的多组学package支持，一个代表了Python，为大数据和AI而生，都有其各自独特的优势，不可被取代。
作为资深用户就很烦了，懂一个是不行的，必须两个都用，所以R和Python语言你必须都精通，否则非常影响项目的数据分析进展。
我也是最先接触Python的，但后面就专心与R，他们语法存在一定的混淆和冲突，所以我只能放弃一个。经过整个博士的训练，我自认为已经精通R了，没有任何问题能把我难倒，Debug小王子。
但现在我不得不捡起Python，因为Scanpy是在是不得不用，以下就记录一些Python中对应的R的核心功能，以便实时查阅。
Seurat对象转Scanpy对象h5ad
【只能导出一个layer，RNA或者ATAC，好处就是metadata都在】默认是导出scale.data的slot，需要指定一下。
library(sceasy)
library(reticulate)
sceasy::convertFormat(seuset.flt.allen.D21, from="seurat", to="anndata", assay = "RNA",
outFile='keyRdata/seuset.flt.allen.D21.RNA.h5ad')
sceasy::convertFormat(escc.seuset, from="seurat", to="anndata", assay = "RNA", main_layer = "count",
outFile='escc.seuset_RNA.h5ad')
Scanpy读取h5ad
adata_ref = sc.read_h5ad(os.path.join(work_dir, '/home/zz950/projects/ApcKO_multiomics/keyRdata/seuset.flt.allen.D21.RNA.h5ad'))
Seurat读取h5ad
【不要用SeuratDisk的LoadH5Seurat，就用最原始的anndata，只需要count matrix和metadata即可】
library("Seurat")
library("anndata")
data <- read_h5ad("Kang_2022_Stomach_Reference.h5ad")
stomach.ref <- CreateSeuratObject(counts = t(data$X), meta.data = data$obs)
print(str(stomach.ref))
Scanpy的数据结构
adata.var adata.obs adata.X adata.raw
查看原始数据
pd.DataFrame.sparse.from_spmatrix(gc.X[0:10,0:10])
Scanpy基本操作
UMAP改颜色
import matplotlib as mpl sc.pl.umap(adata, color=['Sox9', 'Ly6a', 'Tacstd2'], palette="Set2", color_map=mpl.cm.Reds)
import seaborn as sns
blue = sns.color_palette('light:blue', as_cmap = True)
Continuous color schemes are given with the color_map argument, categorical schemes are given with palette. All the scatter plots (scatter, pca, tsne, umap, etc...) share these arguments.
Scanpy obs metadata操作
- 创建一个新的df
- 变换cellName
- 行或列取子集
- 按行或列合并两个df
- 某一列的字符串split
- 取两个vector交集
- 正则匹配
创建一个新的df
import pandas as pd
ref_df = pd.DataFrame(index=a)
ref_df["rawName"] = b
# or
d = {'col1': [1, 2], 'col2': [3, 4]}
df = pd.DataFrame(data=d, index=a)
修改列名
chromsizes.columns = ["Chromosome", "Start", "End"]
按行取子集
adata = adata[:, var_names]
按列取子集
类似seurat和R里的%in%函数，非常常用
adata_subset = adata[adata.obs['celltype'].isin(['a', 'b'])]
按行合并两个df
类似R里的cbind，按索引合并
merged_df = adata.obs.merge(ref_anno, left_index=True, right_index=True)
按列拼装两个df
类似R里的rbind
df1.append(df2)
某一列的字符串split
类似R里的strsplit函数，提取自己需要的信息 【这个python比R的语法简单】
[i.split("-")[0] for i in adata_ref.obs.index]
取两个vector交集
var_names = adata_ref.var_names.intersection(adata.var_names)
元素计数，类似R的table函数
pd.value_counts(gene_sets.term)
human和mouse基因大小转换
fetal_genes[0].str.upper().values
R float的matrix转化为integer
df_float_to_int <- function(counts){
counts<-apply(counts,2,function(x) {storage.mode(x) <- 'integer'; x})
return(counts)
}
参考：
scanpy pipeline： http://localhost:17435/notebooks/projects/MRK60/single-cell-org/Python-visualization.ipynb
