单细胞数量太多可以抽样也可以
我分享过 对单细胞表达矩阵做gsea分析的代码，是不同单细胞亚群两两之间差异分析后，对基因进行排序，非常正常的gsea分析。以及 单细胞转录组数据的批量GSVA代码大放送，是根据单细胞亚群分组后使用AverageExpression得到一个简单的表达量矩阵后进行gsva分析，把2万多个基因的表达量矩阵转换为几十或者上百个 通路的基因集打分矩阵。
但是有同学提问，它的单细胞表达量矩阵是五万到十万个细胞，并不想预先拆分成为单细胞亚群分组，所以没办法使用AverageExpression得到一个简单的表达量矩阵，想直接对全部的单细胞矩阵进行gsva，但是矩阵每次都会内存溢出，大家也可以尝试下面的代码：
# write.csv(t(as.matrix(sce.all@assays$RNA@counts)), file = "sce.all.csv")
一般来说都会报错，除非你的机器配置惊人：
Error in asMethod(object) :
Cholmod error 'problem too large' at file ../Core/cholmod_dense.c, line 102
其实也有一个简单的解决方案，就是直接抽取十分之一细胞即可，代码如下所示：
sce.tmp = sce.all[, sample(1:ncol(sce.all),round(ncol(sce.all)/10)) ]
tmp=as.matrix(sce.tmp@assays$RNA@counts)
write.csv(t(tmp), file = "sce.all.csv")
其实， 如果有单细胞亚群注释信息，我们仍然是推荐拆分成为不同单细胞亚群，如下所示：
load(file = 'phe-by-markers.Rdata')
table(phe$celltype)
sce.all$celltype = phe$celltype
#拆分为 多个 seurat子对象
sce.all.list <- SplitObject(sce.all , split.by = "celltype")
sce.all.list
names(sce.all.list)
for (i in names(sce.all.list)) {
mat = sce.all.list[[i]]@assays$RNA@counts
write.csv(t(as.matrix( mat )), file = paste0(i,".sce.all.csv"))
}
一般来说，哪怕是五万到十万个细胞的矩阵，拆分成为不同单细胞亚群，每个亚群都是少于一万个细胞的，就可以很容易转变为真正的矩阵存储在R里面啦。
但是，如果是下面的这样的情况就麻烦大了 ，总共是三百多万个单细胞，每个亚群仍然是有几十万个单细胞，比普通人一个单细胞转录组项目的细胞数量还有多一个数量级。
为什么要输出csv文件呢
其实是因为有一些后续分析步骤在其它编程语言里面完成，比如在Python里面的转录因子分析。大家可以再次复习一下前面的笔记：pyscenic的转录因子分析结果展示之5种可视化 ，回顾了一下 单细胞转录因子分析之SCENIC流程 ，需要重新认识了 使用pyscenic做转录因子分析 后的结果。
如果是多个单细胞亚群各自的csv文件，就需要写一个脚本接受输入输出文件了，在Linux环境里面写一个 Python脚本 （ csv2loom.py ）把 csv格式的表达量矩阵 转为 .loom 文件 ，代码如下所示：
import os, sys
os.getcwd()
os.listdir(os.getcwd())
input_csv = sys.argv[1]
output_loom = sys.argv[2]
import loompy as lp;
import numpy as np;
import scanpy as sc;
x=sc.read_csv(input_csv);
row_attrs = {"Gene": np.array(x.var_names),};
col_attrs = {"CellID": np.array(x.obs_names)};
lp.create(output_loom,x.X.transpose(),row_attrs,col_attrs);
有两个这个Python脚本，就需要在Linux环境下面运行：
conda activate pyscenic
# 如果是单个文件：
python csv2loom.py tmp.csv.gz tmp.loom
# 如果是多个文件，可以走下面的代码：
ls *csv.gz|while read id;do( python csv2loom.py $id ${id%%.*}.loom );done
把每个单细胞亚群的csv格式的表达量矩阵批量转变为loom格式后走 使用pyscenic做转录因子分析 的流程。
学徒作业
对pbmc3k这个经典的单细胞表达量矩阵，根据单细胞亚群注释信息，拆分成为不同的csv格式的表达量矩阵后，独立走 使用pyscenic做转录因子分析 流程，然后跟整个矩阵的 使用pyscenic做转录因子分析 流程的结果对比！
大家先安装这个数据集对应的包，并且对它进行降维聚类分群，参考前面的例子：人人都能学会的单细胞聚类分群注释 ，而且每个亚群找高表达量基因，都存储为Rdata文件。
# 0.安装R包 ----
# devtools::install_github('caleblareau/BuenColors')
# utils::install.packages(pkgs = "ggstatsplot")
# InstallData("pbmc3k")
library(SeuratData) #加载seurat数据集
getOption('timeout')
options(timeout=10000)
#InstallData("pbmc3k")
data("pbmc3k")
sce <- pbmc3k.final
library(Seurat)
table(Idents(sce))
p1=DimPlot(sce,label = T)
这个数据集已经进行了不同单细胞亚群的标记，接下来就靠大家啦！
期待大家完成作业哦！