talklr_workflow

library(talklr)
library(dplyr)
#> Warning: package 'dplyr' was built under R version 3.5.2
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

Single condition mode

The code below will read in a tab-separated expression table (first column is gene symbol,and remaining columns are mean expression values for each cell type). The expression values must be normalized for sequencing depth, but no need to normalize for gene length. The data MUST NOT be log2 transformed.

library(talklr)
glom_normal<-read.table(system.file("extdata", "glom_normal_data.txt", package = "talklr"),header=T,sep="\t",stringsAsFactors =F)

lr_glom_normal<-make_expressed_net(glom_normal,4,receptor_ligand,'product',1)
lr_glom_normal<-arrange(lr_glom_normal,desc(KL))
head(lr_glom_normal)
#>     Pair.Name Ligand.ApprovedSymbol
#> 1  NPNT_ITGA8                  NPNT
#> 2     F8_LRP1                    F8
#> 3   TNC_ITGA8                   TNC
#> 4 COL4A3_CD93                COL4A3
#> 5   VEGFA_KDR                 VEGFA
#> 6 COL4A5_CD93                COL4A5
#>                                        Ligand.Name Receptor.ApprovedSymbol
#> 1                                     nephronectin                   ITGA8
#> 2  coagulation factor VIII, procoagulant component                    LRP1
#> 3                                       tenascin C                   ITGA8
#> 4 collagen, type IV, alpha 3 (Goodpasture antigen)                    CD93
#> 5             vascular endothelial growth factor A                     KDR
#> 6                       collagen, type IV, alpha 5                    CD93
#>                                                         Receptor.Name DLRP HPMR
#> 1                                                   integrin, alpha 8 <NA> HPMR
#> 2                  low density lipoprotein receptor-related protein 1 <NA> HPMR
#> 3                                                   integrin, alpha 8 <NA> <NA>
#> 4                                                       CD93 molecule <NA> <NA>
#> 5 kinase insert domain receptor (a type III receptor tyrosine kinase) DLRP HPMR
#> 6                                                       CD93 molecule <NA> <NA>
#>   IUPHAR HPRD STRING.binding STRING.experiment HPMR.Ligand HPMR.Receptor
#> 1   <NA> HPRD STRING.binding              <NA>        NPNT         ITGA8
#> 2   <NA> HPRD           <NA>              <NA>          F8          LRP1
#> 3   <NA> <NA> STRING.binding              <NA>        <NA>         ITGA8
#> 4   <NA> HPRD           <NA>              <NA>        <NA>          CD93
#> 5 IUPHAR HPRD STRING.binding STRING.experiment        VEGF           KDR
#> 6   <NA> HPRD           <NA>              <NA>        <NA>          CD93
#>   PMID.Manual Pair.Source        Pair.Evidence ligand_podo ligand_mesa
#> 1        <NA>       known literature supported 1377.510999   18.220567
#> 2        <NA>       known literature supported    5.271237    8.635061
#> 3        <NA>       novel             putative    2.791624   80.301552
#> 4        <NA>       novel literature supported 1660.315851   11.243287
#> 5        <NA>       known literature supported 1806.634648    7.129777
#> 6        <NA>       novel literature supported  118.027705    1.829393
#>   ligand_endo receptor_podo receptor_mesa receptor_endo       KL
#> 1    3.091031     20.823977   1924.149513     12.957356 2.904744
#> 2  598.168237      4.855549    185.055508      3.553831 2.691040
#> 3    1.544840     20.823977   1924.149513     12.957356 2.687865
#> 4   19.130464      2.734476      4.199248    129.519739 2.683975
#> 5   75.417251     10.706085     22.096879   1145.777272 2.683641
#> 6    1.167388      2.734476      4.199248    129.519739 2.639632

We can then visualize the ligand-receptor wiring diagram for any of the top ligand-receptor pairs identified by talklr.

plot_lr_wiring(as.numeric(lr_glom_normal[1,17:19]),as.numeric(lr_glom_normal[1,20:22]),c("podo","mesa","endo"))
#> Loading required package: igraph
#> 
#> Attaching package: 'igraph'
#> The following objects are masked from 'package:dplyr':
#> 
#>     as_data_frame, groups, union
#> The following objects are masked from 'package:stats':
#> 
#>     decompose, spectrum
#> The following object is masked from 'package:base':
#> 
#>     union


plot_lr_wiring(as.numeric(lr_glom_normal[2,17:19]),as.numeric(lr_glom_normal[2,20:22]),c("podo","mesa","endo"))

Optional: comparison with differentially expressed gene-based approach

Typically, studies restrict ligand-receptor pairs to genes expressed 3 fold higher in one cell type compared to each of the other cell types, in other words, cell type marker genes.
In normal glomerulus, there are 198 such ligand-receptor pairs.

lr_deg_glom_normal<-make_deg_net(glom_normal,lr_glom_normal,3,1)
nrow(lr_deg_glom_normal)
#> [1] 198
head(lr_deg_glom_normal)
#>     Pair.Name Ligand.ApprovedSymbol
#> 1  NPNT_ITGA8                  NPNT
#> 2     F8_LRP1                    F8
#> 3   TNC_ITGA8                   TNC
#> 4 COL4A3_CD93                COL4A3
#> 5   VEGFA_KDR                 VEGFA
#> 6 COL4A5_CD93                COL4A5
#>                                        Ligand.Name Receptor.ApprovedSymbol
#> 1                                     nephronectin                   ITGA8
#> 2  coagulation factor VIII, procoagulant component                    LRP1
#> 3                                       tenascin C                   ITGA8
#> 4 collagen, type IV, alpha 3 (Goodpasture antigen)                    CD93
#> 5             vascular endothelial growth factor A                     KDR
#> 6                       collagen, type IV, alpha 5                    CD93
#>                                                         Receptor.Name DLRP HPMR
#> 1                                                   integrin, alpha 8 <NA> HPMR
#> 2                  low density lipoprotein receptor-related protein 1 <NA> HPMR
#> 3                                                   integrin, alpha 8 <NA> <NA>
#> 4                                                       CD93 molecule <NA> <NA>
#> 5 kinase insert domain receptor (a type III receptor tyrosine kinase) DLRP HPMR
#> 6                                                       CD93 molecule <NA> <NA>
#>   IUPHAR HPRD STRING.binding STRING.experiment HPMR.Ligand HPMR.Receptor
#> 1   <NA> HPRD STRING.binding              <NA>        NPNT         ITGA8
#> 2   <NA> HPRD           <NA>              <NA>          F8          LRP1
#> 3   <NA> <NA> STRING.binding              <NA>        <NA>         ITGA8
#> 4   <NA> HPRD           <NA>              <NA>        <NA>          CD93
#> 5 IUPHAR HPRD STRING.binding STRING.experiment        VEGF           KDR
#> 6   <NA> HPRD           <NA>              <NA>        <NA>          CD93
#>   PMID.Manual Pair.Source        Pair.Evidence ligand_podo ligand_mesa
#> 1        <NA>       known literature supported 1377.510999   18.220567
#> 2        <NA>       known literature supported    5.271237    8.635061
#> 3        <NA>       novel             putative    2.791624   80.301552
#> 4        <NA>       novel literature supported 1660.315851   11.243287
#> 5        <NA>       known literature supported 1806.634648    7.129777
#> 6        <NA>       novel literature supported  118.027705    1.829393
#>   ligand_endo receptor_podo receptor_mesa receptor_endo       KL
#> 1    3.091031     20.823977   1924.149513     12.957356 2.904744
#> 2  598.168237      4.855549    185.055508      3.553831 2.691040
#> 3    1.544840     20.823977   1924.149513     12.957356 2.687865
#> 4   19.130464      2.734476      4.199248    129.519739 2.683975
#> 5   75.417251     10.706085     22.096879   1145.777272 2.683641
#> 6    1.167388      2.734476      4.199248    129.519739 2.639632

Two condition mode

In the two condition mode, we will have one expression table for each condition. We will first identify a set of genes expressed in at least one cell type in either condition, and restrict ligand-receptor pairs to those expressed genes.

glom_normal<-read.table(system.file("extdata", "glom_normal_data.txt", package = "talklr"),header=T,sep="\t",stringsAsFactors =F)
glom_fsgs<-read.table(system.file("extdata", "glom_fsgs_data.txt", package = "talklr"),header=T,sep="\t",stringsAsFactors =F)
expressed_normal<-rowSums(glom_normal[,2:ncol(glom_normal)]>4)>=1
expressed_fsgs<-rowSums(glom_fsgs[,2:ncol(glom_fsgs)]>4)>=1
expressed_genes<-glom_normal$genes[expressed_normal|expressed_fsgs]

We will then identify interesting ligand receptor pairs for each condition, and then identify ligand-receptor pairs that changed between two conditions.

normal_net<-make_expressed_net_specify_expressed_genes(glom_normal,expressed_genes,receptor_ligand,'product')
fsgs_net<-make_expressed_net_specify_expressed_genes(glom_fsgs,expressed_genes,receptor_ligand,'product')
normal_net$fsgs_vs_normal_KL<-disease_vs_normal_KL(fsgs_net[,17:19],fsgs_net[,20:22],normal_net[,17:19],normal_net[,20:22],1,'produt')
perturbed_net<-arrange(normal_net,desc(fsgs_vs_normal_KL))
head(perturbed_net)
#>      Pair.Name Ligand.ApprovedSymbol                     Ligand.Name
#> 1   CCL2_ACKR2                  CCL2  chemokine (C-C motif) ligand 2
#> 2  S100A9_TLR4                S100A9 S100 calcium binding protein A9
#> 3 COL8A1_ITGA2                COL8A1    collagen, type VIII, alpha 1
#> 4    TNC_ITGA2                   TNC                      tenascin C
#> 5      C3_CD81                    C3          complement component 3
#> 6 COL4A1_ITGB8                COL4A1      collagen, type IV, alpha 1
#>   Receptor.ApprovedSymbol
#> 1                   ACKR2
#> 2                    TLR4
#> 3                   ITGA2
#> 4                   ITGA2
#> 5                    CD81
#> 6                   ITGB8
#>                                                  Receptor.Name DLRP HPMR IUPHAR
#> 1                                atypical chemokine receptor 2 DLRP <NA>   <NA>
#> 2                                         toll-like receptor 4 <NA> <NA>   <NA>
#> 3 integrin, alpha 2 (CD49B, alpha 2 subunit of VLA-2 receptor) <NA> <NA>   <NA>
#> 4 integrin, alpha 2 (CD49B, alpha 2 subunit of VLA-2 receptor) <NA> <NA>   <NA>
#> 5                                                CD81 molecule <NA> <NA>   <NA>
#> 6                                             integrin, beta 8 <NA> HPMR   <NA>
#>   HPRD STRING.binding STRING.experiment HPMR.Ligand HPMR.Receptor PMID.Manual
#> 1 HPRD           <NA>              <NA>        CCL2         CCBP2        <NA>
#> 2 <NA>           <NA>              <NA>        <NA>          <NA>    17767165
#> 3 HPRD           <NA>              <NA>        <NA>         ITGA2        <NA>
#> 4 <NA> STRING.binding              <NA>        <NA>         ITGA2        <NA>
#> 5 <NA> STRING.binding              <NA>          C3          CD81        <NA>
#> 6 <NA>           <NA>              <NA>      COL4A1         ITGB8        <NA>
#>   Pair.Source        Pair.Evidence ligand_podo ligand_mesa ligand_endo
#> 1       known literature supported   0.6655229   39.096000   0.9128985
#> 2       novel literature supported   0.0000000    4.030385   0.2326046
#> 3       novel literature supported   6.3640955  213.143616   2.5748886
#> 4       novel             putative   1.7916236   79.301552   0.5448398
#> 5       novel             putative   1.0782504    3.063415   0.5117477
#> 6       known literature supported   7.7050899  463.994434 118.0979251
#>   receptor_podo receptor_mesa receptor_endo        KL fsgs_vs_normal_KL
#> 1      0.346268     0.5806588     45.714711 2.7353417         1.2699917
#> 2      2.315512   133.4317970     47.590678 1.9460942         0.7622556
#> 3     80.654362     4.3548455      2.036168 2.4467662         0.6982411
#> 4     80.654362     4.3548455      2.036168 2.5152248         0.6815130
#> 5   1501.702615  1066.4289633   1201.783928 0.3637203         0.6505875
#> 6     97.189415     0.9618120      0.139797 2.2564404         0.6082486

Finally, we can then visualize rewiring of ligand-receptor pairs between two conditions.

pair2plot<-'CCL2_ACKR2' 
par(mfrow=c(1,2))
par(mar = c(0.3, 0.3, 0.3, 0.3))
plot_lr_wiring(as.numeric(normal_net[normal_net$Pair.Name==pair2plot,17:19]),as.numeric(normal_net[normal_net$Pair.Name==pair2plot,20:22]),c("podo","mesa","endo"))

plot_lr_wiring(as.numeric(fsgs_net[fsgs_net$Pair.Name==pair2plot,17:19]),as.numeric(fsgs_net[fsgs_net$Pair.Name==pair2plot,20:22]),c("podo","mesa","endo"))