作者输瓜，Evil Genius

世事無常

315打假日

2月26日獲知家里電信詐騙芹缔，到今日過去了17天，從一開始的震驚拌阴，到冷靜也僅用了3天， 我特別感謝那些幫助了我的人，很多人無償捐助了我欠痴。

很多人都還是學(xué)生迄靠，將來都會走向社會，進(jìn)入崗位喇辽，其中有一些人也會遇到很大的挫折掌挚，我希望大家遇到挫折的時候可以想起我，我這么倒霉的情況下菩咨，依然要相信吠式，生活還是很美好的，大多數(shù)人的挫折比起我來旦委，也就不再是挫折了奇徒。

今天我們來分享一個關(guān)于蛋白活性推斷的內(nèi)容，最近一段時間因為一篇文章的發(fā)表缨硝，運(yùn)用基因表達(dá)來推斷蛋白活性，文章在Single-cell protein activity analysis identifies recurrence-associated renal tumor macrophages罢低，雜志 Cell查辩，頂刊，其中就用到了單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)來推斷蛋白活性网持，其中用到的軟件是viper,2021年5月的一個軟件宜岛，值得關(guān)注。

推斷原理

VIPER(Virtual Inference of Protein-activity by Enriched Regulon analysis)算法允許在單個樣本的基礎(chǔ)上功舀，從基因表達(dá)譜數(shù)據(jù)計算蛋白質(zhì)活性推斷萍倡。它利用最直接受特定蛋白質(zhì)調(diào)控的基因表達(dá)，如轉(zhuǎn)錄因子(TF)的靶標(biāo)辟汰，作為其活性的準(zhǔn)確推斷手段列敲。

viper實現(xiàn)了一種專門用于估計調(diào)控活性的算法，該算法考慮了調(diào)節(jié)子的作用模式帖汞， regulator-target gene相互作用的可信度和每個靶基因調(diào)控的多效性戴而。

VIPER在這個包中提供了兩種推斷方法:多樣本版本(msVIPER)設(shè)計用于基于多個樣本或表達(dá)譜的基因表達(dá)特征，以及單樣本版本(VIPER)翩蘸，它在逐個樣本的基礎(chǔ)上估計相對蛋白質(zhì)活性所意，從而允許將典型的基因表達(dá)矩陣(即多個樣本中的多個mRNA)轉(zhuǎn)換為蛋白質(zhì)活性矩陣，表示每個樣本中每個蛋白質(zhì)的相對活性。

看一下實例代碼

安裝扶踊，其中bcellViper提供了示例數(shù)據(jù)和需要的調(diào)控網(wǎng)絡(luò)作為參考

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

需要輸入兩個文件

• gene expression signature
• an appropriate cell context-specific regulatory network. This regulatory network is provided in the format of a class object

調(diào)控文件通常是由 ARACNe的輸出文件產(chǎn)生的泄鹏，我們來看一下分析過程：

• 由于ARACNe的資源消耗問題，所有對于單細(xì)胞數(shù)據(jù)針對每個cluster進(jìn)行計算

• 為了生成準(zhǔn)確的秧耗、魯棒性好的ARACNe network命满，ARACNe需要輸入表達(dá)矩陣中細(xì)胞的大部分轉(zhuǎn)錄結(jié)構(gòu)相同的數(shù)據(jù)。對于單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)而言绣版，這需要在生成ARACNe network之前將數(shù)據(jù)中的細(xì)胞進(jìn)行clustering胶台。這個cluster可以通過多種方式獲取：任何一種用于單細(xì)胞聚類分群的方法都可以杂抽，也可以是簡單的通過前幾個主成分進(jìn)行的簡單聚類分群诈唬。PISCES包中的Clustering方法有：Partition Around Medioids (PAM), Multi-Way K-Means, and Louvain with Resolution Optimization。PISCES軟件使用的是基于Seurat與PISCES R Package 對數(shù)據(jù)進(jìn)行的two step optimize resolution cluster：所有的clustering step均分兩步完成缩麸。Seurat中FindNeighbors與FindClusters函數(shù)使用的是Louvain算法铸磅，這種算法的缺陷是會導(dǎo)致過度分群。因此杭朱，在0.01 ~ 1.0的分辨率（resolution）范圍內(nèi)進(jìn)行聚類以0.01為間隔阅仔，并在每個分辨率值上評估聚類質(zhì)量，以在此范圍內(nèi)選擇最佳的聚類方式弧械。對于每個分辨率值（resolution）八酒，將cluster中的細(xì)胞數(shù)再次取樣至1000，并計算這1000個細(xì)胞及其cluster標(biāo)簽的silhouette score刃唐。對于基因表達(dá)數(shù)據(jù)羞迷，Pearson correlation被用于細(xì)胞距離矩陣（就是用CorDist函數(shù)算出來的dist.mat）被用于計算silhouette score。對于VIPER推斷出來的蛋白活性數(shù)據(jù)画饥，VIPER包中的ViperSimilarity函數(shù)計算出的distance metric被用于計算silhouette score衔瓮。這個過程會針對這1000個細(xì)胞隨機(jī)進(jìn)行100次，然后得出一個針對一個resolution的mean and standard deviation of average silhouette score抖甘。選擇使平均silhouette score最大化的最高resolution值作為對數(shù)據(jù)進(jìn)行聚類而不過度聚類的最佳resolution热鞍。
  Clustering完成后就可以產(chǎn)生meta-cells用于輸入ARACNe：將cluster中距離最近的10個細(xì)胞的reads相加后進(jìn)一步re-normalizing,生成一個具有250個sample的矩陣用于后續(xù)ARACNe（這個地方操作起來有點復(fù)雜，資源可以的話將所有細(xì)胞輸入進(jìn)行計算）衔彻。
  如果在數(shù)據(jù)集中的不同細(xì)胞類型已經(jīng)進(jìn)行了定義和注釋薇宠，那么cell-type specific networks可以基于細(xì)胞注釋得出。然而米奸，由于無監(jiān)督下的（無細(xì)胞定義與注釋）PISCES 計算可以進(jìn)一步確認(rèn)實驗的設(shè)計是否有問題并可能進(jìn)一步得出新的生物學(xué)發(fā)現(xiàn)昼接，因此推薦無監(jiān)督下的PISCES 計算。這里大家就用定義好的Seurat分群就可以了悴晰。

# 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的矩陣進(jìn)行下游推斷即可逐工，如果cluster的細(xì)胞數(shù)過多，則需要下采用或者合并漂辐。

ARACNe-network generation

從這里開始ARACNe-AP的教程泪喊，這部分內(nèi)容可以基于windows平臺的Terminal或者linux平臺完成，這里使用的是linux平臺髓涯，方法大同小異袒啼。

• 先下載基于ARACNe-AP作者放在Github上的ARACNe-AP壓縮包，將解壓后的文件上傳至linux服務(wù)器并放在后續(xù)ARACNe工作目錄下并解壓纬纪；
• 分別下載并安裝JAVA和ANT蚓再，并設(shè)置相應(yīng)環(huán)境變量
  下載JAVA-JDK：https://www.oracle.com/java/technologies/downloads/#java8
  下載ANTs:https://ant.apache.org/bindownload.cgi
• Xshell或MAC-Terminal登錄linux服務(wù)器

mkdir JAVA_JDK#創(chuàng)建JAVA安裝位置
mkdir ANTs    #創(chuàng)建ANTs安裝位置

• 將下載好的JDK與ANTs安裝文件分別傳送至服務(wù)器的“/JAVA_JDK”與“/ANTs”文件夾中
• 安裝軟件:分別進(jìn)入到/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

• 配置環(huán)境文件包各，進(jìn)入到個人目錄下
  vi .profile #使用文本編輯器打開并編輯“.profile”文件

• 分別為JAVA與ANT增加環(huán)境變量,編輯一下內(nèi)容加入到“.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

• 進(jìn)入ARACNe工作目錄摘仅，檢查JAVA與ANT是否安裝成功, 如果安裝成功輸入以下命令會彈出詳細(xì)版本號

source .profile
java -version
ant -version

• 構(gòu)建JAVA需要的aracne.jar文件

ant main

• 構(gòu)建ARACNe-network需要兩個文件：表達(dá)矩陣（Matrix）與基因列表(list)，注意表達(dá)矩陣中使用ensembl ID问畅，那么基因列表也需要是ensembl ID, 如果矩陣中使用gene symble,那么基因列表也需要是gene symble娃属。這兩個文件都放在ARACNe-AP-master.zip解壓后的test文件夾中，運(yùn)行參數(shù)中的文件名一定要與你文件夾中的文件名一致护姆，否則java在調(diào)用這個文件的時候會報錯
• 表達(dá)矩陣：單細(xì)胞基因表達(dá)矩陣矾端，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的表達(dá)矩陣都需要跑三遍ARACNe流程，再將最后的三個文件合并成一個柒凉，即得到這個cluster的ARACNe-network 。當(dāng)然基因集可以合并

#Step1 Calculate threshold with a fixed seed
java -Xmx5G -jar dist/aracne.jar -e 單細(xì)胞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 單細(xì)胞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 單細(xì)胞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 單細(xì)胞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

• 最后在當(dāng)前工作目錄下進(jìn)入“outputFolder", 里面的”network.txt“就是要的結(jié)果篓跛，我們將一個矩陣針對tfs-hugo.txt膝捞、cotfs-hugo.txt、sig-hugol.txt的三次運(yùn)行結(jié)果合并在一起愧沟，并使用文本打開network.txt文件蔬咬，將其表頭“Regulator Target MI pvalue”刪除，保存后修改文件后綴為tsv沐寺，就可以將其輸入RegProcess()函數(shù)了林艘。最終得到的文件就是這個矩陣的ARACNe-network，幾個矩陣就有幾個ARACNe-network混坞。

• 然后將輸出結(jié)果中的 .adj file和基因表達(dá)矩陣狐援，轉(zhuǎn)換成 object.

Protein activity based clustering analysis
構(gòu)建的ARACNe-network通過RegProcess()函數(shù)將每一個cluster的ARACNe-network文件（pisces-1-net-final.tsv）與其相應(yīng)的表達(dá)矩陣基因表達(dá)矩陣（MetaCells函數(shù)生成的表達(dá)矩陣）進(jìn)行綜合并生成一個調(diào)節(jié)子對象文件（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()函數(shù)生成的文件保存為rds格式，先使用readRDS()函數(shù)讀進(jìn)來啥酱。
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()函數(shù)構(gòu)建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對象里去，并設(shè)置為當(dāng)前默認(rèn)的active.assay
sce.combined.sct <- CPMTransform(sce.combined.sct)#針對PISCESassay進(jìn)行CPMTransform()镶殷、GESTransform()標(biāo)準(zhǔn)化禾酱。
sce.combined.sct <- GESTransform(sce.combined.sct)
sce.combined.sct <- PISCESViper(sce.combined.sct, net.list)#蛋白活性推斷

推斷出最終的VIPER矩陣（蛋白活性矩陣），就可以對細(xì)胞重新進(jìn)行clustering绘趋。VIPER結(jié)果通常允許分辨出更小的颤陶、轉(zhuǎn)錄上更不同的cluster。然后這些cluster可以被用于鑒定master regulator protein陷遮。也可以運(yùn)用蛋白矩陣進(jìn)行下游的再分析以及更多的個性化分析滓走。

左：差異基因熱圖；右：差異蛋白表達(dá)熱圖

生活很好拷呆，有你更好闲坎。