STRINGdb分析单细胞亚群蛋白相互作用网络

研究蛋白之间的相互作用网络，有助于挖掘核心的调控基因。String是一个蛋白互作网络数据库，包含蛋白直接物理作用的互作关系的同时，还包含了蛋白之间以间接作用的互作关系。在线的界面很友好，可根据需求输入不同的protein name or protein sequence来做预测。在细胞数据分析中，我们经常需要那一个基因list（一般是差异基因）看看他们有什么特点，在STRINGdb中就可以仅输入基因名来获得可能的蛋白互做关系。

StringDB有多种可用物种，如人类、小鼠、斑马鱼、C.elengas and D.melanogaster等，可以对网络进行不同的修剪。为了达到我们的目的，我们使用了置信度最高的边。

我们可以用这个R包做：

从STRINGdb获取人网络/图。
对网络进行修剪，只得到高置信度边
创建邻接矩阵
将网络中的蛋白质id映射到邻接矩阵中的基因id

首先我们先载入需要的R包和我们熟悉的pbmc3k.final数据，并计算一个差异基因。

require(STRINGdb)
require(igraph)
require(biomaRt)
library(Seurat)
library(SeuratData)
library(tidyverse)

mk <- pbmc3k.final %>% FindAllMarkers(only.pos = T,logfc.threshold = .05)

然后，我们探究一下这个R包的功能：

 STRINGdb$methods()
 [1] ".objectPackage"                      ".objectParent"                      
 [3] "add_diff_exp_color"                  "add_proteins_description"           
 [5] "benchmark_ppi"                       "benchmark_ppi_pathway_view"         
 [7] "callSuper"                           "copy"                               
 [9] "enrichment_heatmap"                  "export"                             
[11] "field"                               "get_aliases"                        
[13] "get_annotations"                     "get_annotations_desc"               
[15] "get_bioc_graph"                      "get_clusters"                       
[17] "get_enrichment"                      "get_graph"                          
[19] "get_homologs"                        "get_homologs_besthits"              
[21] "get_homology_graph"                  "get_interactions"                   
[23] "get_link"                            "get_neighbors"                      
[25] "get_pathways_benchmarking_blackList" "get_png"                            
[27] "get_ppi_enrichment"                  "get_ppi_enrichment_full"            
[29] "get_proteins"                        "get_pubmed"                         
[31] "get_pubmed_interaction"              "get_subnetwork"                     
[33] "get_summary"                         "get_term_proteins"                  
[35] "getClass"                            "getRefClass"                        
[37] "import"                              "initFields"                         
[39] "initialize"                          "load"                               
[41] "load_all"                            "map"                                
[43] "mp"                                  "plot_network"                       
[45] "plot_ppi_enrichment"                 "post_payload"                       
[47] "remove_homologous_interactions"      "set_background"                     
[49] "show"                                "show#envRefClass"                   
[51] "trace"                               "untrace"                            
[53] "usingMethods"

想用?来获取帮助文档，居然是失败的，但是可用STRINGdb$help("map")来查看map函数的介绍及说明。下面我们来看看如何用STRINGdb来绘制蛋白相互作用网络。我们可以选中某个感兴趣的亚群的差异基因的top多少来做，这样做的意义是从蛋白互作上解释该一群的特异性。

首先我们画出这个亚群：


DimPlot(object = pbmc3k.final, reduction = "umap",order =T ,  sizes.highlight = .005,
        cells.highlight = WhichCells (pbmc3k.final,idents = "Naive CD4 T"), )+
  theme(legend.position = "none") + coord_fixed(1:1) + ggtitle(label = "Naive CD4 T")

然后取出它的差异基因绘制蛋白相互作用网络图。

NCD4T <- mk %>% filter( cluster == "Naive CD4 T") %>% top_n(20, avg_logFC )
NCD4T

           p_val avg_logFC pct.1 pct.2     p_val_adj     cluster      gene
1  5.612007e-124 0.5203411 0.997 0.975 7.696306e-120 Naive CD4 T     RPS25
2  2.513370e-118 0.5211234 0.994 0.972 3.446835e-114 Naive CD4 T      RPL9
3  1.075148e-108 0.5375120 0.991 0.977 1.474457e-104 Naive CD4 T     RPS3A
4  1.963031e-107 0.7300635 0.901 0.594 2.692101e-103 Naive CD4 T      LDHB
5   4.351761e-98 0.5102529 0.994 0.967  5.968004e-94 Naive CD4 T    RPS27A
6   1.284058e-88 0.5156347 0.990 0.904  1.760958e-84 Naive CD4 T     RPS29
7   1.606796e-82 0.9219135 0.436 0.110  2.203560e-78 Naive CD4 T      CCR7
8   4.198081e-77 0.6598608 0.838 0.406  5.757249e-73 Naive CD4 T      CD3D
9   2.316054e-54 0.5999356 0.726 0.399  3.176237e-50 Naive CD4 T      CD3E
10  3.191555e-50 0.6926087 0.628 0.358  4.376899e-46 Naive CD4 T     NOSIP
11  3.324866e-49 0.7296372 0.336 0.104  4.559722e-45 Naive CD4 T      LEF1
12  2.498572e-44 0.7119736 0.331 0.110  3.426542e-40 Naive CD4 T PRKCQ-AS1
13  7.614450e-43 0.6500625 0.438 0.185  1.044246e-38 Naive CD4 T   PIK3IP1
14  3.318236e-39 0.6020963 0.197 0.041  4.550628e-35 Naive CD4 T      FHIT
15  6.485255e-33 0.6420019 0.263 0.088  8.893879e-29 Naive CD4 T       MAL
16  1.055905e-31 0.5311329 0.385 0.179  1.448068e-27 Naive CD4 T      TCF7
17  7.374411e-30 0.7636255 0.245 0.084  1.011327e-25 Naive CD4 T   LDLRAP1
18  1.487326e-29 0.5330400 0.689 0.494  2.039719e-25 Naive CD4 T   C6orf48
19  5.321844e-18 0.5211056 0.418 0.260  7.298377e-14 Naive CD4 T  C12orf57
20  1.037256e-16 0.5693087 0.214 0.096  1.422493e-12 Naive CD4 T      CD8B

看到一些RP基因，这些我们就不再纳入蛋白互作中了，过滤掉。

NCD4T <- mk %>% filter( cluster == "Naive CD4 T") %>% top_n(30, avg_logFC )  %>%
  filter(gene  %in% setdiff(rownames(pbmc3k.final),c(grep("RP",rownames(pbmc3k.final),value = T))  ))
NCD4T

           p_val avg_logFC pct.1 pct.2     p_val_adj     cluster      gene
1  1.963031e-107 0.7300635 0.901 0.594 2.692101e-103 Naive CD4 T      LDHB
2   1.606796e-82 0.9219135 0.436 0.110  2.203560e-78 Naive CD4 T      CCR7
3   4.198081e-77 0.6598608 0.838 0.406  5.757249e-73 Naive CD4 T      CD3D
4   1.858081e-55 0.5047172 0.910 0.801  2.548173e-51 Naive CD4 T      NPM1
5   2.316054e-54 0.5999356 0.726 0.399  3.176237e-50 Naive CD4 T      CD3E
6   3.191555e-50 0.6926087 0.628 0.358  4.376899e-46 Naive CD4 T     NOSIP
7   3.324866e-49 0.7296372 0.336 0.104  4.559722e-45 Naive CD4 T      LEF1
8   2.498572e-44 0.7119736 0.331 0.110  3.426542e-40 Naive CD4 T PRKCQ-AS1
9   7.614450e-43 0.6500625 0.438 0.185  1.044246e-38 Naive CD4 T   PIK3IP1
10  1.403899e-41 0.5000026 0.841 0.701  1.925307e-37 Naive CD4 T    TMEM66
11  3.318236e-39 0.6020963 0.197 0.041  4.550628e-35 Naive CD4 T      FHIT
12  3.894244e-35 0.5029688 0.597 0.333  5.340566e-31 Naive CD4 T      IL7R
13  6.485255e-33 0.6420019 0.263 0.088  8.893879e-29 Naive CD4 T       MAL
14  5.199576e-32 0.4756840 0.152 0.029  7.130698e-28 Naive CD4 T LINC00176
15  1.055905e-31 0.5311329 0.385 0.179  1.448068e-27 Naive CD4 T      TCF7
16  7.374411e-30 0.7636255 0.245 0.084  1.011327e-25 Naive CD4 T   LDLRAP1
17  1.487326e-29 0.5330400 0.689 0.494  2.039719e-25 Naive CD4 T   C6orf48
18  4.623459e-25 0.4911351 0.397 0.210  6.340612e-21 Naive CD4 T  LEPROTL1
19  3.179825e-18 0.4884449 0.438 0.284  4.360812e-14 Naive CD4 T    NDFIP1
20  5.321844e-18 0.5211056 0.418 0.260  7.298377e-14 Naive CD4 T  C12orf57
21  1.037256e-16 0.5693087 0.214 0.096  1.422493e-12 Naive CD4 T      CD8B
22  3.785959e-16 0.4750692 0.396 0.255  5.192064e-12 Naive CD4 T    DNAJB1


example1_mapped <- string_db$map( NCD4T, "gene", removeUnmappedRows = TRUE ) # 等一会。

下面就可以直接绘制了。

hits <- example1_mapped$STRING_id[1:20]
head(hits)

[1] "9606.ENSP00000229319" "9606.ENSP00000246657" "9606.ENSP00000300692" "9606.ENSP00000296930"
[5] "9606.ENSP00000354566" "9606.ENSP00000375726"

string_db$plot_network( hits )

选择特定的cluster的差异基因做个string_network看看他们之间的关系也是文章信息的一个补充啊。用stringdb作图后，stringdb还提供了一个小功能STRING payload mechanism，就是能根据基因表达的上下调信息对互作网络图上的点进行标注颜色，比如上调标注为红色，下调标注为绿色。

pct1map <- string_db$add_diff_exp_color( subset(NCD4T, pct.1 >0.7),logFcColStr="avg_logFC" )
payload_id <- string_db$post_payload( pct1map$gene,colors=pct1map$color )
string_db$plot_network( hits, payload_id=payload_id )

当然，想要画出好看的网络图cytoscape还是比较好用的。我们可以导出作图文件：

info <- string_db$get_interactions(hits)
head(info)

                  from                   to neighborhood neighborhood_transferred fusion
1 9606.ENSP00000246657 9606.ENSP00000265165            0                        0      0
2 9606.ENSP00000254322 9606.ENSP00000296930            0                        0      0
3 9606.ENSP00000265165 9606.ENSP00000300692            0                        0      0
4 9606.ENSP00000246657 9606.ENSP00000306157            0                        0      0
5 9606.ENSP00000300692 9606.ENSP00000306157            0                        0      0
6 9606.ENSP00000300692 9606.ENSP00000331172            0                        0      0
  cooccurence homology coexpression coexpression_transferred experiments experiments_transferred
1           0        0            0                        0           0                       0
2           0        0            0                        0           0                       0
3           0        0            0                        0           0                       0
4           0        0            0                        0           0                       0
5           0        0            0                        0           0                       0
6           0        0            0                        0         576                       0
  database database_transferred textmining textmining_transferred combined_score
1        0                    0        165                      0            165
2        0                    0        176                      0            175
3        0                    0        163                      0            163
4        0                    0        700                      0            700
5        0                    0        191                      0            191
6      900                    0        176                    164            966

StringDB提供了一种计算GO(过程、功能和成分)、KEGG和反应体途径、PubMed、UniProt和PFAM/INTERPRO/SMARTdomains的富集的方法，所有这些都是在一个简单的调用。富集本身是用hy- per几何检验计算的，而FDR是用Benjamini-Hochberg程序计算的。这些我们就不再一一演示了。感兴趣可以看看他的文档：https://bioconductor.org/packages/release/bioc/vignettes/STRINGdb/inst/doc/STRINGdb.pdf

本文的出发点在表达谱，经由莫一群的差异基因list，通过StringDB 落脚在蛋白相互作用上。也可以说为我们认识细胞亚群提供了一个侧面，特别是我们还可以在string_network图上显著地标出特定的基因，以使我们的故事更加圆满。

看到人家可以把基因list利用到什么程度了吗：

STRING:蛋白质相互作用(PPI网络)数据库简介
 STRING网站+Cytoscape软件制作精美蛋白互作网络图(PPI)
https://bioconductor.riken.jp/packages/3.0/bioc/html/STRINGdb.html

STRINGDB包的简单使用