提高单细胞分析准确度的工具之一：Self-assembling manifolds

目前单细胞分析中使用Seurat的居多，该流程中可以影响下游分群的因素有很多，比如前期的质控、Feature select、PCA数量以及resolution的值等等，每一个环节研究者都可以通过多次的调试找到合适的参数值。

Feature select这一步对下游的影响很重要，seurat默认挑选前2000个高变异的Features作为下游PCA计算的输入数据，根据不同的data，我们可以尝试不同的Feature select数量，也可以将scanpy和seurat的var feature整合后进入下游分析，甚至有软件比如scry和M2Drop等专门用于Feature select这一步。但这里小编要说的是另一个软件Self-assembling manifolds(SAM)，其可以独自完成一个完整的单细胞分析的流程。小编这里只挑出SAM的Feature select这一步的结果，整合到Seurat流程中做分析。

软件的开发者在文章《Self-assembling manifolds in single-cell RNA sequencing data》中说（下图），SAM在多项评估中得分最高，其次为SAM的top rank feature作为输入的Seurat的得分仅次于SAM，这两者和Seurat自身的默认参数（Seurat）和优化参数的Seurat分析(Seurat*)结果有明显的差距。

首先是SAM的安装，有两种方式，一是SAM被封装进scanpy中，可通过scanpy调用，二是独立使用安装使用SAM，这里写的是第二种，参考github地址：https://github.com/atarashansky/self-assembling-manifold

conda create -n SAM python=3.6
conda activate SAM
#pip安装
pip install sam-algorithm
#或者github上安装最新版
git clone https://github.com/atarashansky/self-assembling-manifold.git
cd self-assembling-manifold
python setup.py install
#其他安装
pip install matplotlib
conda install -c conda-forge -c plotly jupyter ipywidgets plotly=4.0.0 colorlover ipyevents

安装好之后，可以打开jupyter notebook进行交互式的分析，此部分读者可以查看github的教程。小编喜欢直接写成python脚本分析并导出feature select信息，最后整合到seurat流程中，脚本如下

import sys
from samalg import SAM
a=sys.argv[1]
sam=SAM()
sam.load_data(a)
sam.preprocess_data(sum_norm='cell_median',thresh=0.01)
sam.run(projection='tsne')
a= sam.adata.var
for i in  range(len(a)) :
        print(a[i-1:i])

执行脚本，输入为csv格式的表达矩阵，可以用seurat质控好后导出的count矩阵：

python sam.py   count.from.seurat.csv  |  grep  -E "False|True"   >>  top.var.feature.txt

将结果导入到R中，整合进seurat的流程中：

tmp.sam.genes<- read.table("top.var.feature.txt",header = F)
tmp.sam.genes<-tmp.sam.genes[order(tmp.sam.genes[,4],decreasing = T),]
tmp.sam.top500<- tmp.sam.genes[1:500,1] 
#在seurat的RunPCA输入SAM选择的top 500 features
seu.obj <- RunPCA(seu.obj, features =tmp.sam.top500, npcs = 50, verbose = FALSE)

这里选取的top Features的数量也要根据具体的数据进行调试。

如果是喜欢用scanpy的同学，可以直接在scanpy里调用SAM：

import scanpy as sc
import scanpy.external as sce
sam = sce.tl.sam(adata, inplace=True)  #adata is your AnnData object

有兴趣的同学可以用自己的数据试一下，最后附上这篇文章的地址：https://elifesciences.org/articles/48994