提高单细胞分析准确度的工具之一：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
如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程
生信爆款入门-第8期（线上直播4周，马拉松式陪伴，带你入门）你的生物信息入门课
数据挖掘学习班第6期（线上直播3周，马拉松式陪伴，带你入门） 医学生/医生首选技能提高课
生信技能树的2019年终总结 你的生物信息成长宝藏
看完记得顺手点个“在看”哦！
长按扫码可关注