基于单细胞测序的蛋白活性推断（PISCES）分析流程笔记，凎！

作者: 喂鸡组的阿翔 | 来源:发表于2022-05-05 14:49 被阅读0次

写在前面

感谢califano-lab团队将代码无私的分享出来=。=！
本流程需要一点点R语言与linux基础，我都是在生信技能树学的=。=！
流程中做的是自己的项目，所以没有把结果图放上来；
有一说一，我想设置付费阅读来着，但是我简书钻不够不让弄=。=！；码字不易，如果对你有帮助的话打个赏吧

What is PISCES?

单细胞蛋白活性推断（The pipeline for Protein Activity Inference in Single Cells，PISCES）是基于网络调控分析单细胞基因表达谱的方法。PISCES将高可变性的、表达谱数据伴有技术噪音的单细胞测序数据转换为稳定且可重复性良好的蛋白活性表达谱。PISCES主要基于两个算法：A.精确的细胞网络重构算法ARACNe；B.基于富集调控分析对蛋白质活性虚拟推断算法metaVIPER。
ARACNe算法是应用最广泛的一种基于基因表达数据推断转录交互作用的算法。而metaVIPER算法通过ARACNe推断的一个既定蛋白的靶向调控基因（如其转录因子等）的表达数据作为此蛋白活性高低的依据。一般来说，PISCES可以从单细胞表达数据精确的评估多至6000个以上的蛋白活性，显著提高了针对单个细胞测序中部分不可检测到的mRNA的生物作用及其相关基因产物的分析能力。

PISCES workflow

1.加载PISCES 需要用到的功能包:

rm(list = ls())
options(warn = -1)
options(stringsAsFactors = F)

library(viper)
library(Seurat)
library(ggplot2)
library(ggpubr)
library(pheatmap)
library(RColorBrewer)
library(MUDAN)
library(umap)
library(igraph)
library(biomaRt)
library(uwot)
library(Matrix)
library(PISCES)

2.加载单细胞数据。

此教程中默认使用基因的symbol进行分析，也可以使用ENSMBL ID；
这里加载的单细胞数据是使用标准Seurat process进行质控后，并进一步使用SCTransform方法进行整合后的Seurat Object，详见Seurat官网整合数据教程；
PISCES老一版的官方教程中使用的是CPM转换后的标准化数据，新一版教程则使用的是SCT转换后的标准化数据；

#加载质控并SCTransform整合后的Seurat Object
load(file = "sce_af_integrate.Rdata")

3.对细胞进行分群

由于后续使用ARACNe的运行需要消耗大量的计算资源与时间，PISCES作者推荐基于细胞分群参数的设定来进行ARACNe计算；
为了生成准确的、鲁棒性好的ARACNe network，ARACNe需要输入表达矩阵中细胞的大部分转录结构相同的数据。对于单细胞转录组数据而言，这需要在生成ARACNe network之前将数据中的细胞进行clustering。这个cluster可以通过多种方式获取：任何一种流行的用于单细胞聚类分群的方法都可以，也可以是简单的通过前三十个主成分进行的简单聚类分群。PISCES包中的Clustering方法有：Partition Around Medioids (PAM), Multi-Way K-Means, and Louvain with Resolution Optimization。PISCES作者使用的是基于Seurat与PISCES R Package 对数据进行的two step optimize resolution cluster：所有的clustering step均分两步完成。Seurat中FindNeighbors与FindClusters函数使用的是Louvain算法，这种算法的缺陷是会导致过度分群。因此，在0.01 ~ 1.0的分辨率（resolution）范围内进行聚类以0.01为间隔，并在每个分辨率值上评估聚类质量，以在此范围内选择最佳的聚类方式。对于每个分辨率值（resolution），将cluster中的细胞数再次取样至1000，并计算这1000个细胞及其cluster标签的silhouette score。对于基因表达数据，Pearson correlation被用于细胞距离矩阵（就是用CorDist函数算出来的dist.mat）被用于计算silhouette score。对于VIPER推断出来的蛋白活性数据，VIPER包中的ViperSimilarity函数计算出的distance metric被用于计算silhouette score。这个过程会针对这1000个细胞随机进行100次，然后得出一个针对一个resolution的mean and standard deviation of average silhouette score。选择使平均silhouette score最大化的最高resolution值作为对数据进行聚类而不过度聚类的最佳resolution。
Clustering完成后就可以产生meta-cells用于输入ARACNe：将cluster中距离最近的10个细胞的reads相加后进一步re-normalizing,生成一个具有250个sample的矩阵用于后续ARACNe。
PISCES 默认数据集中没有已知的细胞类型。如果在数据集中的不同细胞类型已经进行了定义和注释，那么cell-type specific networks可以基于细胞注释得出。然而，由于无监督下的（无细胞定义与注释）PISCES 计算可以进一步确认实验的设计是否有问题并可能进一步得出新的生物学发现，因此PISCES 作者推荐无监督下的PISCES 计算。

# Seurat clustering
sce.combined.sct <- RunPCA(sce.combined.sct, verbose = FALSE)
sce.combined.sct <- RunUMAP(sce.combined.sct, dims = 1:30, verbose = FALSE)
sce.combined.sct <- FindNeighbors(sce.combined.sct, dims = 1:30, verbose = FALSE)
sce.combined.sct <- FindClusters(sce.combined.sct, method = "igraph" ,verbose = FALSE)
# PISCES clustering
sce.combined.sct <- CorDist(sce.combined.sct, pca.feats = 10)#运行时间长
sce.combined.sct <- LouvainKRange(sce.combined.sct, kmax = 100)#运行时间长

*我自己使用的单细胞数据集有60000+个细胞，因此上述代码中LouvainKRange()函数运行过程中会报错，简单检索后猜测是计算内容的矢量超出了R语言的上限（2^31-1）,尝试减少细胞数量可以顺利运行，但咱不能无缘无故把质控后保留的细胞给删了；由于报错的内容和igraph这个包相关，进一步检索后发现有可能是LouvainKRange()这个函数在运行过程中把Seurat Object中的稀疏矩阵转换为了正常矩阵交给了igraph包中的graph_from_adjacency_matrix()函数。尝试过后发现将交给graph_from_adjacency_matrix()函数的对象转换为稀疏矩阵即可解决，代码如下很简单：

trace('LouvainClust', edit = T, where = asNamespace("PISCES"))
#在弹出的对话框中将这一段代码：
“graph.obj <- igraph::graph_from_adjacency_matrix(adj.mat)”
#改为
“graph.obj <- igraph::graph_from_adjacency_matrix(as(adj.mat, "sparseMatrix"))”
#保存后再次运行即可

这里我们就得到了基于PISCES的cell clusters,下一步使用MetaCells()函数为每一个cluster生成一个具有250个cells的表达矩阵用于后续ARACNe-network的构建；

4.构建ARACNe-network需要使用的cluster表达矩阵,下面是官方给的代码：

gexp.dist <- sce.combined.sct@assays$SCT@misc$dist.mat
pisces.metacell.mats <- MetaCells(as.matrix(sce.combined.sct@assays$RNA@counts),
                                  dist.mat = gexp.dist,
                                  clust.vect = sce.combined.sct@assays$SCT@misc$pisces.cluster)

如果你的数据是通过SCT转换后并使用Seurat包中的IntegrateData()函数整合的，那么上述代码就需要修改，因为你的pisces.cluster数据在Seurat Object中的存储对象不在SCT这个位置，integrated这个位置（具体需要你对Seurat Object这个对象的存储方式进行学习，网上资料很多，例如“生信技能树=。=！”），代码如下：

gexp.dist <- sce.combined.sct@assays$SCT@misc$dist.mat
pisces.metacell.mats <- MetaCells(as.matrix(sce.combined.sct@assays$RNA@counts),
                                  dist.mat = gexp.dist,
                                  clust.vect = sce.combined.sct@assays$integrated@misc$pisces.cluster)

5.保存Pisces-cluster结果并写出用于ARACNe-network构建

dir.create('pisces-clust_aracne-mats/')#当前工作路径下创建pisces-clust_aracne-mats文件夹
#将cluster表达矩阵以RDS的格式存储到pisces-clust_aracne-mats文件夹中
for (mcn in names(pisces.metacell.mats)) {
  saveRDS(pisces.metacell.mats[[mcn]], 
          file = paste('pisces-clust_aracne-mats/pisces-', 
                       mcn, '-arac.rds', sep = ''))
}
#你的PISCES clustering得到了几群细胞，这里你就会得到几个表达矩阵，每个表达矩阵中含有250个细胞样本

#再次读取表达矩阵，这里演示假如我得到了3群细胞的表达矩阵
pisces1 <- readRDS('pisces-clust_aracne-mats/pisces-1-arac.rds')
pisces2 <- readRDS('pisces-clust_aracne-mats/pisces-2-arac.rds')
pisces3 <- readRDS('pisces-clust_aracne-mats/pisces-3-arac.rds')
#存储为CSV格式
write.csv(pisces1,file = "pisces-1-arac.csv")
write.csv(pisces2,file = "pisces-2-arac.csv")
write.csv(pisces3,file = "pisces-3-arac.csv")

*将保存好的csv格式的文件修改为txt格式的文件用于后续ARACNe-network构建

6.ARACNe-network generation

从这里开始followARACNe-AP的tutorial，这部分内容可以基于windows平台的Terminal或者linux平台完成，我这里使用的是linux平台（租用的服务器），方法大同小异。以下使用的linux命令及其基础有不会的检索后都能查到。

先下载基于ARACNe-AP作者放在Github上的ARACNe-AP压缩包，使用fileZilla软件将解压后的文件上传至linux服务器并放在后续ARACNe工作目录下并解压；
分别下载并安装JAVA和ANT，并设置相应环境变量
下载JAVA-JDK：https://www.oracle.com/java/technologies/downloads/#java8
下载ANTs:https://ant.apache.org/bindownload.cgi
Xshell或MAC-Terminal登录linux服务器

mkdir JAVA_JDK#创建JAVA安装位置
mkdir ANTs    #创建ANTs安装位置

使用FileZilla将下载好的JDK与ANTs安装文件分别传送至服务器的“/JAVA_JDK”与“/ANTs”文件夹中
安装软件:分别进入到/JAVA_JDK”与“/ANTs”文件夹中，解压文件

tar -zxvf /your/pathway/to/apache-ant-1.9.14-bin.tar.gz
tar -zxvf /your/pathway/to/jdk-8u131-linux-x64.tar.gz

配置环境文件，进入到自己账户的目录下（我是租用共享服务器）

ls #查看是否有“.profile”文件
vim .profile #使用文本编辑器打开并编辑“.profile”文件

按“i”进入插入模式分别为JAVA与ANT增加环境变量,编辑一下内容加入到“.profile”文件

export JAVA_HOME=/YOUR/PATHWAY/TO/jdk1.8.0_211 
export CLASSPATH=$:CLASSPATH:$JAVA_HOME/lib/ 
export PATH=$PATH:$JAVA_HOME/bin

export ANT_HOME=/YOUR/PATHWAY/TO/apache-ant-1.9.15/
export PATH=$ANT_HOME/bin:$PATH

进入ARACNe工作目录，检查JAVA与ANT是否安装成功, 如果安装成功输入以下命令会弹出详细版本号

source .profile
java -version
ant -version

构建JAVA需要的aracne.jar文件

ant main

构建ARACNe-network需要两个文件：表达矩阵（Matrix）与基因列表(list)，注意表达矩阵中使用ensembl ID，那么基因列表也需要是ensembl ID, 如果矩阵中使用gene symble,那么基因列表也需要是gene symble。这两个文件都放在ARACNe-AP-master.zip解压后的test文件夹中，运行参数中的文件名一定要与你文件夹中的文件名一致，否则java在调用这个文件的时候会报错
表达矩阵：我们上一步已经得到：“pisces-1-arac.txt、pisces-2-arac.txt、pisces-3-arac.txt”；
基因列表：我们使用PISCES作者给的："tfs-hugo.txt", COTFs-"cotfs-hugo.txt",Signaling Proteins-"sig-hugol.txt"，因为是三个基因集（文件可以在PISCES作者的github上下载：https://github.com/califano-lab/PISCES/tree/master/data，他这里除了有我们上述提到的3个基因集，还有一个surface基因集，本教程演示过程没有使用），因此我们针对每一个cluster的表达矩阵都需要跑三遍ARACNe流程，再将最后的三个文件合并成一个，即得到这个cluster的ARACNe-network 。

#Step1 Calculate threshold with a fixed seed
java -Xmx5G -jar dist/aracne.jar -e test/pisces-1-arac.txt  -o outputFolder --tfs test/tfs.txt --pvalue 1E-8 --seed 1 \--calculateThreshold
#Step2: Run ARACNe on a single bootstrap
java -Xmx5G -jar dist/aracne.jar -e test/pisces-1-arac.txt  -o outputFolder --tfs test/tfs.txt --pvalue 1E-8 --seed 1
#Step3: Run 100 reproducible bootstraps
for i in {1..100}
do
java -Xmx5G -jar dist/aracne.jar -e test/pisces-1-arac.txt  -o outputFolder --tfs test/tfs.txt --pvalue 1E-8 --seed $i
done
#Step4: Consolidate bootstraps in the output folder
java -Xmx5G -jar dist/aracne.jar -o outputFolder --consolidate
#Step5: Run a single ARACNE with no bootstrap and no DPI
java -Xmx5G -jar dist/aracne.jar -e test/pisces-1-arac.txt  -o outputFolder --tfs test/tfs.txt --pvalue 1E-8 --seed 1 \ --nobootstrap --nodpi
#Step6: Consolidate bootstraps without Bonferroni correction
java -Xmx5G -jar dist/aracne.jar -o outputFolder --consolidate --nobonferroni

最后在当前工作目录下进入“outputFolder", 里面的”network.txt“就是我们要的结果，我们将一个cluster针对tfs-hugo.txt、cotfs-hugo.txt、sig-hugol.txt的三次运行结果合并在一起，并使用文本打开network.txt文件，将其表头“Regulator Target MI pvalue”删除，保存后修改文件后缀为tsv，就可以将其输入RegProcess()函数了。最终得到的文件就是这个cluster的ARACNe-network，几个cluster就有几个ARACNe-network。
如果你的cluster比较多，上述构建network的流程可以写一个含有嵌套循环的bash脚本在linux中执行，这样不管你有多少个cluster，针对一个基因集只运行一次上述代码即可，如下：
先在Linux中使用vim 命令编辑一个文本文件并命名为："aracen3.sh",文本内容如下：

#!/bin/bash
for j in {1..14}
do
java -Xmx5G -jar dist/aracne.jar -e test/pisces-$j-arac.txt  -o outputFolder_$j --tfs test/sig-hugo.txt --pvalue 1E-8 --seed 1 \--calculateThreshold

java -Xmx5G -jar dist/aracne.jar -e test/pisces-$j-arac.txt  -o outputFolder_$j --tfs test/sig-hugo.txt --pvalue 1E-8 --seed 1

for i in {1..100} 
do 
java -Xmx5G -jar dist/aracne.jar -e test/pisces-$j-arac.txt  -o outputFolder_$j --tfs test/sig-hugo.txt --pvalue 1E-8 --seed $i 
done

java -Xmx5G -jar dist/aracne.jar -o outputFolder_$j --consolidate

java -Xmx5G -jar dist/aracne.jar -e test/pisces-$j-arac.txt  -o outputFolder_$j --tfs test/sig-hugo.txt --pvalue 1E-8 --seed 1 \ --nobootstrap --nodpi

java -Xmx5G -jar dist/aracne.jar -o outputFolder_$j --consolidate --nobonferroni
done

chmod +x aracen3.sh #给脚本一个执行权限
nohup ./aracen3.sh & #挂后台运行
jobs -l #查看脚本运行情况及进程号，关闭你的Xshell或MAC-Terminal后再次登录服务器该命令就看不到此进程了
pgrep java #查看脚本是否还在运行，关闭你的Xshell或MAC-Terminal后再次登录服务器该命令依旧可以查看进程
cat nohups.out #查看aracen3.sh运行过程中返回的信息

7.Protein activity based clustering analysis

回到R语言，将你构建的ARACNe-network文件放在R语言当前工作路径下，通过RegProcess()函数将每一个cluster的ARACNe-network文件（pisces-1-net-final.tsv）与其相应的表达矩阵基因表达矩阵（MetaCells函数生成的表达矩阵）进行综合并生成一个调节子对象文件（pisces-1-net-pruned.rds）。

#这里演示还是假如我们有三个cluster
RegProcess('pisces-1-net-final.tsv', mat1, out.dir = 'tutorial/', out.name = 'pisces-1-net-')
RegProcess('pisces-2-net-final.tsv', mat2, out.dir = 'tutorial/', out.name = 'pisces-2-net-')
RegProcess('pisces-3-net-final.tsv', mat3, out.dir = 'tutorial/', out.name = 'pisces-3-net-')

#RegProcess()函数生成的文件保存为rds格式，先使用readRDS()函数读进来。
c1.net <- readRDS('pisces-1-net-pruned.rds')
c2.net <- readRDS('pisces-2-net-pruned.rds')
c3.net <- readRDS('pisces-3-net-pruned.rds')
#使用list()函数构建net.list文件
net.list<-list(c1.net,c2.net,c3.net)

## viper and protein-activity based clustering
## net.list here would be a list of networks generated from ARACNe
sce.combined.sct <- AddPISCESAssay(sce.combined.sct)#将Seurat对象里的count矩阵拿出来命名为PISCES，然后放回到Seurat对象里去，并设置为当前默认的active.assay
sce.combined.sct <- CPMTransform(sce.combined.sct)#针对PISCESassay进行CPMTransform()、GESTransform()标准化。
sce.combined.sct <- GESTransform(sce.combined.sct)
sce.combined.sct <- PISCESViper(sce.combined.sct, net.list)#蛋白活性推断

推断出最终的VIPER矩阵（蛋白活性矩阵），就可以对细胞重新进行clustering。VIPER结果通常允许分辨出更小的、转录上更不同的cluster。然后这些cluster可以被用于鉴定master regulator protein。这可以通过几种方式实现，其中Bootstrapped Mann Whitney-U test的鲁棒性最好，Cluster-specific Stouffer integration、a data-wide ANOVA、Kruskal-Wallis test这几种方法同样也可再PISCES包中找到并使用。

sce.combined.sct <- CorDist(sce.combined.sct, pca.feats = 10)#根据蛋白活性矩阵计算细胞间距离
# 降维
sce.combined.sct <- MakeUMAP(sce.combined.sct)
sce.combined.sct <- MakeMDS(sce.combined.sct)

再次clustering and MR analysis:重新分群，对于VIPER推断出来的蛋白活性数据，VIPER包中的ViperSimilarity函数计算出的distance metric被用于计算silhouette score。这个过程会针对这1000个细胞随机进行100次，然后得出一个针对一个resolution的mean and standard deviation of average silhouette score。选择使平均silhouette score最大化的最高resolution值作为对数据进行聚类而不过度聚类的最佳resolution。

sce.combined.sct <- LouvainKRange(sce.combined.sct, kmin = 5, kmax = 100)
sce.combined.sct <- MWUMrs(sce.combined.sct)#找每个分群的master regulator

#作图
plot.df <- data.frame('MDS1' = sce.combined.sct@assays$PISCES@misc$mds[,1],
                      'MDS2' = sce.combined.sct@assays$PISCES@misc$mds[,2],
                      'UMAP1' = sce.combined.sct@assays$PISCES@misc$umap[,1],
                      'UMAP2' = sce.combined.sct@assays$PISCES@misc$umap[,2],
                      'Cluster' = as.factor(sce.combined.sct@assays$PISCES@misc$pisces.cluster))

umap.plot <- ggplot(plot.df, aes(UMAP1, UMAP2)) + geom_point(aes(color = Cluster)) +
  ggtitle(paste(tissue.title, ' (UMAP)', sep = '')) + plot.theme
print(umap.plot)

mds.plot <- ggplot(plot.df, aes(MDS1, MDS2)) + geom_point(aes(color = Cluster)) +
  ggtitle(paste(tissue.title, ' (MDS)', sep = '')) + plot.theme
print(mds.plot)

参考资料：

[1].https://github.com/califano-lab/PISCES
[2].Ding H, Douglass EF Jr, Sonabend AM, Mela A, Bose S, Gonzalez C, Canoll PD, Sims PA, Alvarez MJ, Califano A. Quantitative assessment of protein activity in orphan tissues and single cells using the metaVIPER algorithm. Nat Commun. 2018 Apr 16;9(1):1471. doi: 10.1038/s41467-018-03843-3. PMID: 29662057; PMCID: PMC5902599.
[3].PISCES: A pipeline for the Systematic, Protein Activity-based Analysis of Single Cell RNA Sequencing Data | bioRxiv
[4].Obradovic A, Chowdhury N, Haake SM, Ager C, Wang V, Vlahos L, Guo XV, Aggen DH, Rathmell WK, Jonasch E, Johnson JE, Roth M, Beckermann KE, Rini BI, McKiernan J, Califano A, Drake CG. Single-cell protein activity analysis identifies recurrence-associated renal tumor macrophages. Cell. 2021 May 27;184(11):2988-3005.e16. doi: 10.1016/j.cell.2021.04.038. Epub 2021 May 20. PMID: 34019793; PMCID: PMC8479759.
[5].https://github.com/califano-lab/ARACNe-AP

基于单细胞测序的蛋白活性推断（PISCES）分析流程笔记，凎！

写在前面

What is PISCES?

PISCES workflow

1.加载PISCES 需要用到的功能包:

2.加载单细胞数据。

3.对细胞进行分群

4.构建ARACNe-network需要使用的cluster表达矩阵,下面是官方给的代码：

5.保存Pisces-cluster结果并写出用于ARACNe-network构建

6.ARACNe-network generation

7.Protein activity based clustering analysis

参考资料：

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

单细胞挖掘思路

单细胞流程