Cell-cell communication analysis
This report summarizes the proof-of-concept for LIANA cell-cell communication analysis
As a pilot data Oelen v3 dataset was selected
Example of LIANA output for aggregated Seurat object is presented below:
all_cell_all_donors_LIANA_scores <- read.table('/Users/korshe/Documents/Data_Groningen/L_R_all_cells.tsv',sep='\t',header=T)
head(all_cell_all_donors_LIANA_scores)## source target ligand receptor aggregate_rank mean_rank
## 1 other monocyte PF4 LRP1 3.306481e-07 155.2
## 2 other DC PF4 THBD 1.433560e-06 243.6
## 3 monocyte monocyte S100A8 CD68 1.625094e-06 82.0
## 4 other monocyte PF4 LDLR 1.728796e-06 152.0
## 5 monocyte monocyte S100A9 CD68 9.298329e-06 82.4
## 6 monocyte monocyte S100A9 CD36 1.153991e-05 133.6
## natmi.edge_specificity natmi.rank connectome.weight_sc connectome.rank
## 1 0.4513651 16 2.514876 3
## 2 0.3682824 30 2.395713 5
## 3 0.2198468 140 1.397114 105
## 4 0.3513807 33 2.367788 6
## 5 0.1688186 251 1.422410 97
## 6 0.1664300 259 1.283525 140
## logfc.logfc_comb logfc.rank sca.LRscore sca.rank cellphonedb.pvalue
## 1 1.663739 19 0.9185359 737 0
## 2 1.565113 31 0.9021337 1151 0
## 3 2.011162 2 0.9475565 162 0
## 4 1.637895 22 0.9202780 698 0
## 5 2.166145 1 0.9545581 62 0
## 6 1.983350 3 0.9419176 265 0
## cellphonedb.rank BH_aggregated_rank Bonf_aggregated_rank cell_types
## 1 1 0.003343845 0.003343845 other_monocyte
## 2 1 0.004370827 0.014497597 other_DC
## 3 1 0.004370827 0.016434576 monocyte_monocyte
## 4 1 0.004370827 0.017483309 other_monocyte
## 5 1 0.018806800 0.094034002 monocyte_monocyte
## 6 1 0.019450523 0.116703139 monocyte_monocyte
## id cells
## 1 other_monocyte_PF4_LRP1 PF4_LRP1
## 2 other_DC_PF4_THBD PF4_THBD
## 3 monocyte_monocyte_S100A8_CD68 S100A8_CD68
## 4 other_monocyte_PF4_LDLR PF4_LDLR
## 5 monocyte_monocyte_S100A9_CD68 S100A9_CD68
## 6 monocyte_monocyte_S100A9_CD36 S100A9_CD36Check what are the most common cell-cell interactions (CCI)
frequency_of_communication <- as.data.frame(table(all_cell_all_donors_LIANA_scores$cell_types))
frequency_of_communication <- frequency_of_communication[order(frequency_of_communication$Freq, decreasing = T),]
frequency_of_communication## Var1 Freq
## 28 NK_monocyte 38
## 35 other_monocyte 35
## 21 monocyte_monocyte 25
## 26 NK_CD8T 21
## 29 NK_NK 12
## 37 other_other 10
## 9 CD4T_monocyte 9
## 36 other_NK 9
## 16 DC_DC 7
## 33 other_CD8T 7
## 15 DC_CD8T 6
## 12 CD8T_CD8T 5
## 17 DC_monocyte 5
## 22 monocyte_NK 5
## 23 monocyte_other 5
## 34 other_DC 5
## 10 CD4T_NK 4
## 30 NK_other 4
## 2 B_CD8T 3
## 11 CD4T_other 3
## 18 monocyte_B 3
## 20 monocyte_DC 3
## 27 NK_DC 3
## 1 B_B 2
## 7 CD4T_CD8T 2
## 3 B_monocyte 1
## 4 B_other 1
## 5 CD4T_B 1
## 6 CD4T_CD4T 1
## 8 CD4T_DC 1
## 13 CD8T_DC 1
## 14 DC_B 1
## 19 monocyte_CD8T 1
## 24 NK_B 1
## 25 NK_CD4T 1
## 31 other_B 1
## 32 other_CD4T 1library(ggplot2)
ggplot(data=frequency_of_communication, aes(reorder(Var1,Freq), y=Freq)) + geom_bar(stat="identity") + coord_flip() + theme_light() + ylab('Number of interactions') + xlab('Cell to cell communication')
To calculate values per donor particular donor was excluded from the matrix and afterwards cell-cell communication scores were calculated. A donor specific value was represented as a result of substracting results form the object, where donor was removed from the aggregated (whole) matrix
Example of the output of [Object without DONOR X]
all_cell_per_donors_LIANA_scores_substract_from_aggregated <- read.table('/Users/korshe/Documents/Data_Groningen/L_R_ALL_CELLSrig_values_with_mean.tsv',sep='\t',header=T)
head(all_cell_per_donors_LIANA_scores_substract_from_aggregated)## id mean_rank LLDeep_0117 LLDeep_1300 LLDeep_0615
## 1 B_B_BTLA_CD79A 1052.0 1039.0 1055.2 1054.4
## 2 B_B_CD70_TNFRSF13B 1406.2 1409.8 1406.6 1432.8
## 3 B_CD8T_HLA-DPB1_LAG3 653.8 654.6 657.4 659.8
## 4 B_CD8T_HLA-DQA1_LAG3 747.6 749.0 752.0 753.2
## 5 B_CD8T_HLA-DRA_LAG3 531.8 529.8 533.4 538.6
## 6 B_monocyte_RPS19_C5AR1 331.0 332.4 331.2 329.2
## LLDeep_0923 LLDeep_0705 LLDeep_1313 LLDeep_1198 LLDeep_0022 LLDeep_0062
## 1 1046.6 1051.4 1052.6 1064.8 1055.0 1045.8
## 2 1419.2 1404.4 1405.0 1391.0 1405.8 1412.2
## 3 652.2 657.8 649.8 648.6 658.2 659.2
## 4 745.8 752.4 742.4 738.2 755.2 747.6
## 5 533.0 537.2 528.2 527.4 535.8 532.4
## 6 332.8 332.8 330.0 331.2 330.2 330.8
## LLDeep_0773 LLDeep_1171 LLDeep_0064 LLDeep_1081 LLDeep_0492 LLDeep_0032
## 1 1036.6 1049.8 1050.6 1048.0 1041.4 1054.4
## 2 1396.6 1414.0 1410.6 1383.4 1401.4 1402.4
## 3 651.4 655.6 655.8 653.2 650.4 657.0
## 4 740.2 754.2 750.4 742.0 747.8 749.8
## 5 528.2 531.6 533.2 525.4 531.0 534.6
## 6 331.6 328.8 332.6 319.6 328.6 328.6
## LLDeep_0458 LLDeep_0526 LLDeep_1191 LLDeep_0717 LLDeep_1318 LLDeep_1045
## 1 1051.0 1067.8 1083.6 1056.6 1044.6 1049.2
## 2 1402.4 1419.0 1435.4 1401.8 1408.2 1412.0
## 3 652.8 649.0 652.2 647.8 651.8 646.4
## 4 748.4 743.4 740.4 736.2 744.2 742.6
## 5 529.6 528.8 524.2 526.0 529.4 523.4
## 6 328.6 332.6 331.8 330.4 335.6 328.4
## LLDeep_0853 LLDeep_0768 LLDeep_1133 LLDeep_1189 LLDeep_0854 LLDeep_1160
## 1 1052.4 1046.0 1052.6 1034.2 1052.8 1073.2
## 2 1402.2 1383.0 1406.4 1379.4 1396.8 1452.2
## 3 655.0 654.4 654.4 652.8 659.4 650.6
## 4 748.6 741.8 751.2 749.0 757.0 743.4
## 5 531.4 531.6 536.8 533.0 538.2 529.8
## 6 328.0 334.8 328.0 336.6 330.6 332.4
## LLDeep_1289 LLDeep_0986 LLDeep_1197 LLDeep_0244 LLDeep_0810
## 1 1042.8 1055.6 1043.8 1047.8 1051.0
## 2 1392.8 1389.8 1394.4 1398.2 1439.8
## 3 648.0 652.4 651.2 657.0 649.2
## 4 742.2 743.0 745.0 751.8 747.8
## 5 525.4 529.0 530.8 533.6 529.6
## 6 327.6 332.8 328.0 330.2 333.0Overall 243 of significant interactions was found
dim(all_cell_per_donors_LIANA_scores_substract_from_aggregated)## [1] 243 34Example of the output of [Whole Object] - [Object without DONOR X]
all_cell_per_donors_LIANA_scores <- read.table('/Users/korshe/Documents/Data_Groningen/L_R_ALL_CELLS_substract_values.tsv',sep='\t',header=T)
head(all_cell_per_donors_LIANA_scores)## id mean_rank LLDeep_0117 LLDeep_1300 LLDeep_0615
## 1 B_B_BTLA_CD79A 1052.0 13.0 -3.2 -2.4
## 2 B_B_CD70_TNFRSF13B 1406.2 -3.6 -0.4 -26.6
## 3 B_CD8T_HLA-DPB1_LAG3 653.8 -0.8 -3.6 -6.0
## 4 B_CD8T_HLA-DQA1_LAG3 747.6 -1.4 -4.4 -5.6
## 5 B_CD8T_HLA-DRA_LAG3 531.8 2.0 -1.6 -6.8
## 6 B_monocyte_RPS19_C5AR1 331.0 -1.4 -0.2 1.8
## LLDeep_0923 LLDeep_0705 LLDeep_1313 LLDeep_1198 LLDeep_0022 LLDeep_0062
## 1 5.4 0.6 -0.6 -12.8 -3.0 6.2
## 2 -13.0 1.8 1.2 15.2 0.4 -6.0
## 3 1.6 -4.0 4.0 5.2 -4.4 -5.4
## 4 1.8 -4.8 5.2 9.4 -7.6 0.0
## 5 -1.2 -5.4 3.6 4.4 -4.0 -0.6
## 6 -1.8 -1.8 1.0 -0.2 0.8 0.2
## LLDeep_0773 LLDeep_1171 LLDeep_0064 LLDeep_1081 LLDeep_0492 LLDeep_0032
## 1 15.4 2.2 1.4 4.0 10.6 -2.4
## 2 9.6 -7.8 -4.4 22.8 4.8 3.8
## 3 2.4 -1.8 -2.0 0.6 3.4 -3.2
## 4 7.4 -6.6 -2.8 5.6 -0.2 -2.2
## 5 3.6 0.2 -1.4 6.4 0.8 -2.8
## 6 -0.6 2.2 -1.6 11.4 2.4 2.4
## LLDeep_0458 LLDeep_0526 LLDeep_1191 LLDeep_0717 LLDeep_1318 LLDeep_1045
## 1 1.0 -15.8 -31.6 -4.6 7.4 2.8
## 2 3.8 -12.8 -29.2 4.4 -2.0 -5.8
## 3 1.0 4.8 1.6 6.0 2.0 7.4
## 4 -0.8 4.2 7.2 11.4 3.4 5.0
## 5 2.2 3.0 7.6 5.8 2.4 8.4
## 6 2.4 -1.6 -0.8 0.6 -4.6 2.6
## LLDeep_0853 LLDeep_0768 LLDeep_1133 LLDeep_1189 LLDeep_0854 LLDeep_1160
## 1 -0.4 6.0 -0.6 17.8 -0.8 -21.2
## 2 4.0 23.2 -0.2 26.8 9.4 -46.0
## 3 -1.2 -0.6 -0.6 1.0 -5.6 3.2
## 4 -1.0 5.8 -3.6 -1.4 -9.4 4.2
## 5 0.4 0.2 -5.0 -1.2 -6.4 2.0
## 6 3.0 -3.8 3.0 -5.6 0.4 -1.4
## LLDeep_1289 LLDeep_0986 LLDeep_1197 LLDeep_0244 LLDeep_0810
## 1 9.2 -3.6 8.2 4.2 1.0
## 2 13.4 16.4 11.8 8.0 -33.6
## 3 5.8 1.4 2.6 -3.2 4.6
## 4 5.4 4.6 2.6 -4.2 -0.2
## 5 6.4 2.8 1.0 -1.8 2.2
## 6 3.4 -1.8 3.0 0.8 -2.0dim(all_cell_per_donors_LIANA_scores)## [1] 243 34Checking variance of cell-communication scores
library(Matrix)
library(matrixStats)
all_cell_per_donors_LIANA_scores_matrix <- (all_cell_per_donors_LIANA_scores[,-2])
rownames(all_cell_per_donors_LIANA_scores_matrix) <- all_cell_per_donors_LIANA_scores_matrix[,1]
all_cell_per_donors_LIANA_scores_matrix <- all_cell_per_donors_LIANA_scores_matrix[,-1]
mean <- rowMeans(all_cell_per_donors_LIANA_scores_matrix)
median <- rowMedians(as.matrix(all_cell_per_donors_LIANA_scores_matrix))
sd_st <- rowSds(as.matrix(all_cell_per_donors_LIANA_scores_matrix))
variance_st <- rowVars(as.matrix(all_cell_per_donors_LIANA_scores_matrix))
cv_st <- sd_st/mean
liana_stats <- as.data.frame(rownames(all_cell_per_donors_LIANA_scores_matrix))
colnames(liana_stats) <-'cell_to_cell_interaction'
liana_stats$mean <- mean
liana_stats$median <- median
liana_stats$sd <- sd_st
liana_stats$variance <- variance_st
liana_stats$cv <- cv_st <- sd_st/mean
liana_stats<- liana_stats[order(liana_stats$sd, decreasing = T),]
# all_cell_per_donors_LIANA_scores_matrix<- lapply(all_cell_per_donors_LIANA_scores_matrix, function(x) as.numeric(as.character(x)))
#all_cell_per_donors_LIANA_scores_matrix <- as.numeric(all_cell_per_donors_LIANA_scores_matrix)
head(liana_stats)## cell_to_cell_interaction mean median sd variance cv
## 9 CD4T_CD4T_B2M_CD3D 6.13125 1.6 37.10173 1376.5383 6.051251
## 240 other_other_PF4_THBD -1.31875 4.7 35.17069 1236.9777 -26.669720
## 28 CD4T_other_IL32_ITGB3 1.05000 2.3 20.44846 418.1394 19.474720
## 237 other_other_LTBP1_ITGB5 1.48750 2.7 19.68474 387.4889 13.233437
## 186 other_DC_HMGB1_THBD 0.81250 -2.3 18.76867 352.2631 23.099905
## 242 other_other_TLN1_ITGB5 0.12500 -0.2 18.53414 343.5142 148.273087Visualizing the most variable features
library(tidyr)##
## Attaching package: 'tidyr'## The following objects are masked from 'package:Matrix':
##
## expand, pack, unpackall_cell_per_donors_LIANA_scores_matrix_2 <- all_cell_per_donors_LIANA_scores[,-2]
all_cell_per_donors_LIANA_scores_matrix_long <- gather(all_cell_per_donors_LIANA_scores_matrix_2, donor, CCI, 2:ncol(all_cell_per_donors_LIANA_scores_matrix_2))Distribution of CCI terms among the donors
library(forcats)
all_cell_per_donors_LIANA_scores_matrix_long %>%
ggplot( aes(x=reorder(donor, CCI), y=CCI, fill=donor)) +
geom_boxplot() +
xlab("class") +
theme(legend.position="none") +
xlab("") + coord_flip()
all_cell_per_donors_LIANA_scores_matrix_long_variable <- all_cell_per_donors_LIANA_scores_matrix_long[all_cell_per_donors_LIANA_scores_matrix_long$id %in% liana_stats$cell_to_cell_interaction[1:9],]The most variable CCI terms are:
all_cell_per_donors_LIANA_scores_matrix_long_variable %>%
ggplot( aes(x=reorder(donor, CCI), y=CCI, fill=donor)) +
geom_point() +
xlab("class") +
theme(legend.position="none") +
xlab("") + theme(axis.text.x = element_blank())+
xlab('donor') + facet_wrap(~id, scale="free") #+ theme_light()# + theme(legend.position = "none")#+ coord_flip() 
#
# ggplot(data, aes(x=variety, y=note, fill=treatment)) +
# geom_boxplot() +
# facet_wrap(~variety, scale="free")all_cell_per_donors_LIANA_scores_matrix_long_variable <- all_cell_per_donors_LIANA_scores_matrix_long[all_cell_per_donors_LIANA_scores_matrix_long$id %in% liana_stats$cell_to_cell_interaction[10:29],]The next set of variable terms:
all_cell_per_donors_LIANA_scores_matrix_long_variable %>%
ggplot( aes(x=reorder(donor, CCI), y=CCI, fill=donor)) +
geom_point() +
xlab("class") +
theme(legend.position="none") +
xlab("") + theme(axis.text.x = element_blank())+
xlab('donor') + facet_wrap(~id, scale="free") #+ theme_light()# + theme(legend.position = "none")#+ coord_flip() 
#
# ggplot(data, aes(x=variety, y=note, fill=treatment)) +
# geom_boxplot() +
# facet_wrap(~variety, scale="free")Let’s check if the donors or features cluster together
mtscaled <- as.matrix(scale(all_cell_per_donors_LIANA_scores_matrix))
# create heatmap and don't reorder columns
heatmap(mtscaled, Colv=F, scale='none')
As we see, there is no strong enrichment or correlations among the
donors and features, therefore covariate correction might not be needed
now.