软件介绍

STREAM (Single-cell Trajectories Reconstruction, Exploration And Mapping) 软件是一个伪时间分析软件，通过一系列的特征选择，降纬（MLLE），轨迹推断（EPG）等步骤，逐步构建细胞发育轨迹，从而帮助我们更好的理解发育进化过程。

下载地址

github 地址：https://github.com/pinellolab/STREAM

安装本地版有两种方法，docker或者bioconda，我因为要安装在集群上，不具有管理员权限，所以，建议使用bioconda

conda 安装的方法就不说了，可以参考前文小白也看得懂的 ETE-toolkit 构建进化树流程

在终端配置 conda 的安装通道：

  1. conda config --add channels defaults
  2. conda config --add channels bioconda
  3. conda config --add channels conda-forge

创建一个stream的虚拟空间

  1. conda create -n stream python=3.6 stream jupyter
  2. conda activate stream

使用命令行进行一键化分析

输入文件的格式

该流程需要输入三种文件（都是采用 tab 分割的 tsv 文件），分别是细胞表达矩阵，细胞注释名称以及细胞的颜色定义，当然，后面两个非必须，但是，还是提供比较好。

细胞表达矩阵示例：

Name AAACCTGGTATCAGTC-1_1 AACCGCGCATTCCTGC-1_1 AACGTTGAGCCCAATT-1_1
GNLY 0.142776811271601 0.181584363054937 0.152404882566299
PAEP 0.28157594020077 0.26124097727991 0.33343093647266
MT1G -0.0174682667681339 1.30917646946733 0.193603016988391
HBB …… …… ……
CD74 …… …… ……

表达矩阵可以是原始的 counts 值，也可以是经过标准化后的数据

细胞列表示例：
MDSC
MDSC
Dendritic.cells.activated
MDSC
MDSC

这里不需要输入标题栏，细胞名称和细胞 barcode 一一对应

颜色示例：
DC2 #66c2a5
MDSC #fc8d62
Dendritic.cells.activated #8da0cb

你还可以选择一些有意识的基因出来做后续分析，这里不再展示。

可能很多人要问，如何从 Seurat 对象中获取该信息呢？翻一下官方文档可以知道：

  1. library(Seurat)
  2. Seurat.obj<-readRDS("seurat.rds")
  3. #获取表达矩阵前，应该将默认矩阵转换为整合多样本的矩阵
  4. DefaultAssay(Seurat.obj) <- "integrated"
  5. expMat<-GetAssayData(Seurat.obj,slot="data")#data就是标准化后的数据，你也可以用counts
  6. expMat<-as.data.frame(expMat)
  7. expMat$name<-rownames(expMat)
  8. expMat<-expMat[,c(ncols(expMat),1:(ncols(expMat)-1))]
  9. write.table(expMat,"data_Nestorowa.tsv",sep="\t",quote=FALSE,row.names=FALSE)
  10. #获取细胞信息就简单了
  11. cell<-as.data.frame(as.character(Idents(Seurat.obj)))
  12. write.table(cell,"cell_label.tsv",sep="\t",quote=FALSE,row.names=FALSE,col.names=FALSE)
  13. #映射颜色
  14. cell_c<-data.frame(cell=c("DC2","MDSC","Dendritic.cells.activated"),col=c("#66c2a5","#fc8d62","#8da0cb"))
  15. write.table(cell_c,"cell_label_color.tsv",sep="\t",quote=FALSE,col.names=FALSE)

命令行参数

  1. -m, --matrix
  2. input file name. Matrix is in .tsv or tsv.gz format in which each row represents a unique gene and each column is one cell. (default: None)
  3. -l, --cell_labels
  4. file name of cell labels (default: None)
  5. -c, --cell_labels_colors
  6. file name of cell label colors (default: None)
  7. -s, --select_features
  8. LOESS, PCA, all: Select variable genes using LOESS or top principal components using PCA or keep all the gene (default: LOESS)
  9. --TG
  10. detect transition genes automatically
  11. --DE
  12. detect DE genes automatically
  13. --LG
  14. etect leaf genes automatically
  15. -g, --gene_list
  16. genes to visualize, it can either be filename which contains all the genes in one column or a set of gene names separated by comma (default: None)
  17. -p, --use_precomputed
  18. use precomputed data files without re-computing structure learning part
  19. --log2
  20. perform log2 transformation
  21. --norm
  22. normalize data based on library size
  23. --atac
  24. indicate scATAC-seq data
  25. --n_processes
  26. Specify the number of processes to use. (default, all the available cores).
  27. --loess_frac
  28. The fraction of the data used in LOESS regression (default: 0.1)
  29. --pca_first_PC
  30. keep first PC
  31. --pca_n_PC
  32. The number of selected PCs (default: 15)
  33. --n_processes
  34. Specify the number of processes to use. The default uses all the cores available
  35. --lle_neighbours
  36. LLE neighbour percent (default: 0.1)
  37. --lle_components
  38. Number of components for LLE space (default: 3)
  39. --clustering
  40. Clustering method used for seeding the intial structure, choose from 'ap','kmeans','sc'.
  41. --damping
  42. Affinity Propagation: damping factor (default: 0.75)
  43. --n_clusters
  44. Number of clusters for spectral clustering or kmeans
  45. --EPG_n_nodes
  46. Number of nodes for elastic principal graph (default: 50)
  47. --EPG_lambda
  48. lambda parameter used to compute the elastic energy (default: 0.02)
  49. --EPG_mu
  50. mu parameter used to compute the elastic energy (default: 0.1)
  51. --EPG_trimmingradius
  52. maximal distance of point from a node to affect its embedment (default: Inf)
  53. --EPG_alpha
  54. positive numeric, alpha parameter of the penalized elastic energy (default: 0.02)
  55. --disable_EPG_optimize
  56. disable optimizing branching
  57. --EPG_collapse
  58. Collapsing small branches
  59. --EPG_collapse_mode
  60. the mode used to collapse branches. It can be 'PointNumber','PointNumber_Extrema', 'PointNumber_Leaves','EdgesNumber' or 'EdgesLength' (default:'PointNumber')
  61. --EPG_collapse_par
  62. the control parameter used for collapsing small branches
  63. --EPG_shift
  64. shift branching point
  65. --EPG_shift_mode
  66. the mode to use to shift the branching points 'NodePoints' or 'NodeDensity' (default: NodeDensity)
  67. --EPG_shift_DR
  68. positive numeric, the radius used when computing point density if EPG_shift_mode is 'NodeDensity' (default:0.05)
  69. --EPG_shift_maxshift
  70. positive integer, the maximum distance (number of edges) to consider when exploring the branching point neighborhood (default:5)
  71. --disable_EPG_ext
  72. disable extending leaves with additional nodes
  73. --EPG_ext_mode
  74. the mode used to extend the graph. It can be 'QuantDists', 'QuantCentroid' or 'WeigthedCentroid'. (default: QuantDists)
  75. --EPG_ext_par
  76. the control parameter used for contribution of the different data points when extending leaves with nodes (default: 0.5)
  77. --DE_zscore_cutoff
  78. Differentially Expressed Genes z-score cutoff (default: 2)
  79. --DE_logfc_cutoff
  80. Differentially Expressed Genes log fold change cutoff (default: 0.25)
  81. --TG_spearman_cutoff
  82. Transition Genes Spearman correlation cutoff (default: 0.4)
  83. --TG_logfc_cutoff
  84. Transition Genes log fold change cutoff (default: 0.25)
  85. --LG_zscore_cutoff
  86. Leaf Genes z-score cutoff (default: 1.5)
  87. --LG_pvalue_cutoff
  88. Leaf Genes p value cutoff (default: 1e-2)
  89. --umap
  90. Whether to use UMAP for visualization (default: No)
  91. -r
  92. root node for subwaymap_plot and stream_plot (default:None)
  93. --stream_log_view
  94. use log2 scale for y axis of stream_plot
  95. --for_web
  96. Output files for website
  97. -o, --output_folder
  98. Output folder (default: None)
  99. --new
  100. file name of data to be mapped (default: None)
  101. --new_l
  102. filename of new cell labels (default: None)
  103. --new_c
  104. filename of new cell label colors (default: None)

使用默认参数运行：

  1. stream -m data_Nestorowa.tsv -l cell_label.tsv -c cell_label_color.tsv

建议添加的参数：

  1. stream -m data_Nestorowa.tsv -l cell_label.tsv -c cell_label_color.tsv  --n_processes 4 --umap

这里因为使用的是标准化后的数据，所以跳过了标准化过程：否则添加--norm

进一步可视化基因：

  1. stream -m data_Nestorowa.tsv -l cell_label.tsv -c cell_label_color.tsv -g gene_list.tsv -p

p参数可以避免重复计算，节省时间，默认会把所有节点都作为起始节点（根节点）计算一遍，你可以自己定义根节点-r

探索 Marker 基因：

  1. stream -m data_Nestorowa.tsv -l cell_label.tsv -c cell_label_color.tsv --TG --DG --LG -p

结果

个人对于构建轨迹也不是太理解，就目前的理解来看，比较有用的结果是

  • flat_tree.pdf
  • subway_map.pdf
  • std_vs_means.pdf

使用STREAM对单细胞数据进行伪时间分析 - 图1

这个图是按照默认参数降纬跑的，这里需要注意的是，蓝线必须和点拟合的比较好才行，否则需要适当的调整参数

使用STREAM对单细胞数据进行伪时间分析 - 图2

这张图应该是说，该细胞群由 S0 作为时间起点，向 S1，S2，S3 和 S4 分别进化出四支，就看我们选择哪一个作为根节点了。

使用STREAM对单细胞数据进行伪时间分析 - 图3

基于先验知识，激活的基因座细胞有可能作为发育的起点，所以我们选择 S3 作为根节点，构建出如上的时间序列。

以上，是采用一键化的方式计算，我们还可以使用交互的方式进行，但是因为我的集群没有可视化界面，操作起来就不方便，就没法演示了。具体方法可以参考
STREAM_scRNA-seq.ipynb

可视化网页版链接：http://stream.pinellolab.org