单细胞处理单细胞流程

单细胞数据分析之蛋白活性推断篇

2023-03-14  本文已影响0人  单细胞空间交响乐

作者，Evil Genius

世事无常

315打假日

2月26日获知家里电信诈骗，到今日过去了17天，从一开始的震惊，到冷静也仅用了3天， 我特别感谢那些帮助了我的人，很多人无偿捐助了我。

很多人都还是学生，将来都会走向社会，进入岗位，其中有一些人也会遇到很大的挫折，我希望大家遇到挫折的时候可以想起我，我这么倒霉的情况下，依然要相信，生活还是很美好的，大多数人的挫折比起我来，也就不再是挫折了。

今天我们来分享一个关于蛋白活性推断的内容，最近一段时间因为一篇文章的发表，运用基因表达来推断蛋白活性，文章在Single-cell protein activity analysis identifies recurrence-associated renal tumor macrophages，杂志 Cell，顶刊，其中就用到了单细胞转录组数据来推断蛋白活性，其中用到的软件是viper,2021年5月的一个软件，值得关注。

推断原理

VIPER(Virtual Inference of Protein-activity by Enriched Regulon analysis)算法允许在单个样本的基础上，从基因表达谱数据计算蛋白质活性推断。它利用最直接受特定蛋白质调控的基因表达，如转录因子(TF)的靶标，作为其活性的准确推断手段。

viper实现了一种专门用于估计调控活性的算法，该算法考虑了调节子的作用模式， regulator-target gene相互作用的可信度和每个靶基因调控的多效性。

VIPER在这个包中提供了两种推断方法:多样本版本(msVIPER)设计用于基于多个样本或表达谱的基因表达特征，以及单样本版本(VIPER)，它在逐个样本的基础上估计相对蛋白质活性，从而允许将典型的基因表达矩阵(即多个样本中的多个mRNA)转换为蛋白质活性矩阵，表示每个样本中每个蛋白质的相对活性。

看一下实例代码

安装，其中bcellViper提供了示例数据和需要的调控网络作为参考
if (!requireNamespace("BiocManager", quietly=TRUE))
+ install.packages("BiocManager")
BiocManager::install("mixtools")
BiocManager::install("bcellViper")
BiocManager::install("viper")
Getting started
library(viper)
Generating the regulon object

需要输入两个文件

调控文件通常是由 ARACNe的输出文件产生的，我们来看一下分析过程：
# 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)

这样的话，抽取每个cluster的矩阵进行下游推断即可，如果cluster的细胞数过多，则需要下采用或者合并。

ARACNe-network generation

从这里开始ARACNe-AP的教程，这部分内容可以基于windows平台的Terminal或者linux平台完成，这里使用的是linux平台，方法大同小异。

mkdir JAVA_JDK#创建JAVA安装位置
mkdir ANTs    #创建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

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


source .profile
java -version
ant -version

ant main

#Step1 Calculate threshold with a fixed seed
java -Xmx5G -jar dist/aracne.jar -e 单细胞cluster矩阵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 单细胞cluster矩阵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 单细胞cluster矩阵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 单细胞cluster矩阵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

Protein activity based clustering analysis
构建的ARACNe-network通过RegProcess()函数将每一个cluster的ARACNe-network文件（pisces-1-net-final.tsv）与其相应的表达矩阵基因表达矩阵（MetaCells函数生成的表达矩阵）进行综合并生成一个调节子对象文件（pisces-1-net-pruned.rds）。

RegProcess <- function(a.file, exp.mat, out.dir, out.name = '.') {
  require(viper)
  processed.reg <- aracne2regulon(afile = a.file, eset = exp.mat, format = '3col')
  saveRDS(processed.reg, file = paste(out.dir, out.name, 'unPruned.rds', sep = ''))
  pruned.reg <- pruneRegulon(processed.reg, 50, adaptive = FALSE, eliminate = TRUE)
  saveRDS(pruned.reg, file = paste(out.dir, out.name, 'pruned.rds', sep = ''))
}

#针对每个cluster（3个为例）
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 <- AddProteinAssay(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。也可以运用蛋白矩阵进行下游的再分析以及更多的个性化分析。

左：差异基因热图；右：差异蛋白表达热图

生活很好，有你更好。

上一篇下一篇

猜你喜欢

热点阅读