--- a
+++ b/RNA-seq/AnalysisPipeline/6.3.cellphoneDB.analysis.R
@@ -0,0 +1,533 @@
+#' @description: cellphoneDB analysis
+
+library(psych)
+library(qgraph)
+library(igraph)
+library(tidyverse)
+library(ComplexHeatmap)
+library(circlize)
+library(ggpubr)
+library(dplyr)
+library(vegan)
+library(Seurat)
+
+source(file = "/home/longzhilin/Analysis_Code/code/cellphoneDB.dotplot.R")
+source(file = "/home/longzhilin/Analysis_Code/code/Plot.CircularChart.R")
+
+# TCGA data
+source(file = "/home/longzhilin/Analysis_Code/code/analysis.diff.survival.TCGA.R")
+source(file = "/home/longzhilin/Analysis_Code/code/survival.combined.R")
+DESeq2.normalized_counts <- readRDS("/data/active_data/lzl/RenalTumor-20200713/Data/TCGA/KIRC/Result/DESeq2.normalized_counts.rds")
+DESeq2.normalized_counts <- log2(DESeq2.normalized_counts+1)
+DESeq2.result <- readRDS("/data/active_data/lzl/RenalTumor-20200713/Data/TCGA/KIRC/Result/DESeq2.result.rds")
+clin.data <- readRDS("/data/active_data/lzl/RenalTumor-20200713/Data/TCGA/KIRC/Result/clin.data.rds")
+
+# ICB data
+source(file = "/home/longzhilin/Analysis_Code/SurvivalAnalysis/Cox.function.R")
+source(file = "/home/longzhilin/Analysis_Code/code/RCC.ICB.analysis.R")
+normalized_expression <- readRDS("/data/active_data/lzl/RenalTumor-20200713/Data/ICB.therapy/normalized_expression.log2.rds")
+patient.info.RNA <- readRDS("/data/active_data/lzl/RenalTumor-20200713/Data/ICB.therapy/patient.info.RNA.rds")
+
+patient.info.RNA$Sex <- gsub("^F|FEMALE|Female$", 0, patient.info.RNA$Sex)
+patient.info.RNA$Sex <- gsub("^M|Male|MALE$", 1, patient.info.RNA$Sex)
+patient.info.RNA$Tumor_Sample_Primary_or_Metastasis <- gsub("PRIMARY", 0, patient.info.RNA$Tumor_Sample_Primary_or_Metastasis)
+patient.info.RNA$Tumor_Sample_Primary_or_Metastasis <- gsub("METASTASIS", 1, patient.info.RNA$Tumor_Sample_Primary_or_Metastasis)
+
+data.merge <- readRDS("/data/active_data/lzl/RenalTumor-20200713/DataAnalysis-20210803/scRNA/data.merge.pro.rds")
+DefaultAssay(data.merge) <- "RNA"
+Idents(data.merge) <- data.merge$cellType_low
+Macrophage <- readRDS("/data/active_data/lzl/RenalTumor-20200713/DataAnalysis-20210803/scRNA/5.Immune/Macrophage/sub.scRNA.harmony.pro.rds") # 3 cluster
+CD8.T <- readRDS("/data/active_data/lzl/RenalTumor-20200713/DataAnalysis-20210803/scRNA/5.Immune/CD8T/sub.scRNA.harmony.pro.rds")
+idx1 <- match(colnames(Macrophage), colnames(data.merge)) # 4319 cells
+idx2 <- match(colnames(CD8.T), colnames(data.merge)) # 3728 cells
+cell.Type <- data.merge$cellType_low
+cell.Type[idx1] <- as.character(Macrophage$cellType3)
+cell.Type[idx2] <- paste0("CD8+ T-", CD8.T$cellType3)
+data.merge$cellType2 <- cell.Type
+cell.Type.new <- setdiff(unique(cell.Type), c("CD8+ T cell", "Macrophage"))
+data.merge.new <- subset(data.merge, subset = cellType2 %in% cell.Type.new)
+DefaultAssay(data.merge.new) <- "RNA"
+
+setwd("/data/active_data/lzl/RenalTumor-20200713/DataAnalysis-20210803/scRNA/6.CrossTalk/CellPhoneDB")
+mean.threshold <- 1
+pval.threshold <- 0.05
+colors <- c("#00A087", "#4DBBD5", "#E64B35", "#3C5488")
+
+## Filter custom ligand receptor, partner_a is ligand, partner_b is receptor
+# 3511*7
+ligand.receptor.data <- read.table("/data/ExtraDisk/sdd/longzhilin/Data/signatureGeneSet/Immune/ligand_receptor/ligand_receptor.data.csv", header = T, sep = ",", stringsAsFactors = F)
+ligand_receptor <- paste0(ligand.receptor.data$partner_a, "_", ligand.receptor.data$partner_b)
+
+All.interaction.number <- read.table("interaction_count.txt", header = T, sep = "\t", stringsAsFactors = F)
+cell.types <- All.interaction.number$X
+order.value <- c("CD4+ T cell", "Treg", paste0("CD8+ T-", levels(CD8.T$cellType3)), "Proliferative T cell", "NK/NKT cell", "B cell",
+                 levels(Macrophage$cellType3), "Monocyte", "Dendritic cell", "Neutrophil", "Mast cell",
+                 "Endothelium VCAM1+", "Endothelium VCAM1-", "Mesangial cell", "Tumor")
+data.merge.new$cellType2 <- factor(data.merge.new$cellType2, levels = order.value)
+Idents(data.merge.new) <- data.merge.new$cellType2
+order.value <- gsub("CD8\\+ T-Exhausted IEG", "CD8\\+ T-IEG", order.value)
+####################################Overall Cell Interaction Mode######## ##########################
+#Only keep user_curated, remove the interaction pair in cellphoneDB (including complex)
+mypvals <- read.delim("Renal/pvalues.txt", check.names = FALSE) # 3479 *411
+mymeans <- read.delim("Renal/means.txt", check.names = FALSE)
+sig.means <- read.delim("Renal/significant_means.txt", check.names = FALSE)
+
+mypvals.usr <- mypvals[which(mypvals$annotation_strategy == "user_curated"),] # 2930 * 411
+mymeans.usr <- mymeans[which(mymeans$annotation_strategy == "user_curated"),]
+sig.means.usr <- sig.means[which(sig.means$annotation_strategy == "user_curated"),]
+
+# Identification ligand and receptor, used to identify the first gene
+mypvals.usr$ligand <- rep("Ligand", nrow(mypvals.usr))
+mymeans.usr$ligand <- rep("Ligand", nrow(mymeans.usr))
+sig.means.usr$ligand <- rep("Ligand", nrow(sig.means.usr))
+
+res <- for(i in 1:nrow(mypvals.usr)){
+    index1 <- match(mypvals.usr$interacting_pair[i], ligand_receptor)
+    if(is.na(index1)){
+        mypvals.usr$ligand[i] <- "Receptor"
+        mymeans.usr$ligand[i] <- "Receptor"
+    }
+
+    index2 <- match(sig.means.usr$interacting_pair[i], ligand_receptor)
+    if(is.na(index2)){
+        sig.means.usr$ligand[i] <- "Receptor"
+    }
+}
+# Move the position of the ligand column, position the 11th column
+mypvals.usr <- mypvals.usr[,c(1:10, 412, 11:411)]
+mymeans.usr <- mymeans.usr[,c(1:10, 412, 11:411)]
+sig.means.usr <- sig.means.usr[,c(1:10, 413, 11:412)]
+
+# 2930 ligand-receptors
+saveRDS(mypvals.usr, file = "mypvals.usr.rds")
+saveRDS(mymeans.usr, file = "mymeans.usr.rds")
+saveRDS(sig.means.usr, file = "sig.means.usr.rds")
+write.table(mypvals.usr, file = "Renal/mypvals.usr.txt", row.names = F, col.names = T, sep = "\t", quote = F)
+write.table(mymeans.usr, file = "Renal/mymeans.usr.txt", row.names = F, col.names = T, sep = "\t", quote = F)
+write.table(sig.means.usr, file = "Renal/sig.means.usr.txt", row.names = F, col.names = T, sep = "\t", quote = F)
+
+####Statistics of the ligand-receptor interactions of each cell type
+#ligand
+sig.interaction.cellType <- lapply(cell.types, function(x){
+    mymeans.usr %>% dplyr::select("interacting_pair", "ligand", starts_with(x), ends_with(x)) %>% reshape2::melt() -> meansdf
+    colnames(meansdf)<- c("interacting_pair", "ligand", "CC","means")
+    mypvals.usr %>% dplyr::select("interacting_pair", "ligand", starts_with(x), ends_with(x)) %>% reshape2::melt()-> pvalsdf
+    colnames(pvalsdf)<- c("interacting_pair", "ligand", "CC","pvals")
+    pvalsdf$joinlab<- paste0(pvalsdf$interacting_pair, "_", pvalsdf$CC)
+    meansdf$joinlab<- paste0(meansdf$interacting_pair, "_", meansdf$CC)
+    pldf <- merge(pvalsdf, meansdf, by = "joinlab")
+    sig.interaction <- pldf %>% filter(means >= mean.threshold & pvals < pval.threshold)
+    sig.interaction <- sig.interaction[, -c(6,7,8)]
+    colnames(sig.interaction) <- c("joinlab", "interacting_pair", "ligand", "CC", "pvals", "means")
+
+    ####Split the interaction between ligand and receptor
+    a <- apply(sig.interaction, 1, function(y){
+        CC.interaction <- unlist(strsplit(as.character(y[4]), "\\|"))
+        if((y[3] == "Ligand" & CC.interaction[1] == x) | (y[3] == "Receptor" & CC.interaction[2] == x)){
+            return("ok")
+        }else{
+            return("no")
+        }
+    })
+    ligand.interactions <- sig.interaction[which(a=="ok"),]
+    receptor.interactions <- sig.interaction[which(a=="no"),]
+
+    sig.interaction <- sig.interaction[order(sig.interaction$means, decreasing = T),]
+    ligand.interactions <- ligand.interactions[order(ligand.interactions$means, decreasing = T),]
+    receptor.interactions <- receptor.interactions[order(receptor.interactions$means, decreasing = T),]      
+
+    return(list(All = sig.interaction, Ligand = ligand.interactions, Receptor = receptor.interactions))
+})
+names(sig.interaction.cellType) <- cell.types
+saveRDS(sig.interaction.cellType, file = "sig.interaction.cellType.rds")
+
+a <- lapply(sig.interaction.cellType, function(x){
+    return(x$All)
+})
+library(openxlsx)
+write.xlsx(a, file = "sig.interaction.cellType.xlsx", sheetName = names(sig.interaction.cellType), rowNames = F)
+
+####----Drawing a heat map of interaction
+
+interaction.counts <- sapply(names(sig.interaction.cellType), function(x){
+    interaction.pairs <- sig.interaction.cellType[[x]]$All
+    cellType.counts <- sapply(names(sig.interaction.cellType), function(y){
+        a <- t(as.data.frame(strsplit(x = as.character(interaction.pairs$CC), split = "\\|")))
+        idx1 <- which(a[,1]==y) 
+        idx2 <- which(a[,2]==y)
+        return(length(unique(c(idx1, idx2))))
+    })
+    idx <- which(cellType.counts == max(cellType.counts))
+
+    self.pairs <- length(which(interaction.pairs$CC == paste0(x, "|", x)))
+    cellType.counts[idx] <- self.pairs
+    return(cellType.counts)
+})
+saveRDS(interaction.counts, file = "sig.interaction.counts.rds")
+pdf("sig.interaction.number.heatmap.pdf")
+col_fun <- c("#124F8B", "#FED9B8", "#8B0B50")
+Heatmap(as.matrix(interaction.counts), name = "number of interactions", col = col_fun, row_names_gp = gpar(fontsize = 10), column_names_gp = gpar(fontsize = 10), width = unit(8, "cm"), height = unit(8, "cm"))
+Heatmap(log2(as.matrix(interaction.counts)), name = "number of interactions", col = col_fun, row_names_gp = gpar(fontsize = 10), column_names_gp = gpar(fontsize = 10), width = unit(8, "cm"), height = unit(8, "cm"))
+#只提取immune cell
+idx <- which(rownames(interaction.counts) %in% c("Endothelium VCAM1-", "Endothelium VCAM1+", "Mesangial cell"))
+Heatmap(as.matrix(interaction.counts[-idx, -idx]), name = "number of interactions", col = col_fun, row_names_gp = gpar(fontsize = 10), column_names_gp = gpar(fontsize = 10), width = unit(8, "cm"), height = unit(8, "cm"))
+Heatmap(log2(as.matrix(interaction.counts[-idx, -idx])), name = "number of interactions", col = col_fun, row_names_gp = gpar(fontsize = 10), column_names_gp = gpar(fontsize = 10), width = unit(8, "cm"), height = unit(8, "cm"))
+dev.off()
+
+####-----count the interaction number
+interaction.type.number <- sapply(sig.interaction.cellType, function(x){
+    return(c(nrow(x$Ligand), nrow(x$Receptor)))
+})
+rownames(interaction.type.number) <- c("Ligand", "Receptor")
+
+##1.all cell types
+interaction.type.number <- data.frame(interaction.type.number, check.names = F)
+interaction.type.number$Type <- c("Ligand", "Receptor")
+library(tidyr)
+res <- gather(interaction.type.number, cellType, count, -Type)
+saveRDS(res, file = "Ligand.Receptor.interaction.rds")
+
+res$cellType <- factor(res$cellType, levels = order.value)
+res.ligand <- res[which(res$Type == "Ligand"),]
+res.receptor <- res[which(res$Type == "Receptor"),]
+type <- rep("Lymphoid", nrow(res.ligand))
+type[c(1, 8, 13:16, 19)] <- "Myeloid"
+type[10] <- "Tumor"
+type[c(4, 6, 12)] <- "Other"
+res.ligand$type <- factor(type, levels = c("Lymphoid", "Myeloid", "Tumor", "Other"))
+res.receptor$type <- factor(type, levels = c("Lymphoid", "Myeloid", "Tumor", "Other"))
+pdf("Ligand.Receptor.interaction.pdf")
+ggbarplot(res, x="cellType", y="count", color = "Type", xlab = "", ylab = "Number of interactions", fill = "Type", palette = c("#E18727", "#0072B5"), x.text.angle = 60, position = position_dodge())
+ggbarplot(res.ligand, x="cellType", y="count", color = "type", xlab = "", ylab = "Number of interactions", fill = "type", palette = colors, sort.val = "desc", sort.by.groups = FALSE, x.text.angle = 60)
+ggbarplot(res.receptor, x="cellType", y="count", color = "type", xlab = "", ylab = "Number of interactions", fill = "type", palette = colors, sort.val = "desc", sort.by.groups = FALSE, x.text.angle = 60)
+
+#2.focus on the immune cell
+other <- c("Endothelium VCAM1-", "Endothelium VCAM1+", "Mesangial cell", "Tumor")
+index1 <- which(res$cellType %in% other)
+
+res.immune <- res[-index1,]
+ggbarplot(res.immune, x="cellType", y="count", color = "Type", xlab = "", ylab = "Number of interactions", fill = "Type", palette = c("#E18727", "#0072B5"), x.text.angle = 60, position = position_dodge())
+index1 <- which(res.ligand$cellType %in% other)
+ggbarplot(res.ligand[-index1,], x="cellType", y="count", color = "type", xlab = "", ylab = "Number of interactions", fill = "type", palette = colors, sort.val = "desc", sort.by.groups = FALSE, x.text.angle = 60) 
+ggbarplot(res.receptor[-index1,], x="cellType", y="count", color = "type", xlab = "", ylab = "Number of interactions", fill = "type", palette = colors, sort.val = "desc", sort.by.groups = FALSE, x.text.angle = 60)
+
+# heatmap展示ligand和receptor情况
+ligand.receptor.data <- data.frame(Ligand = res.ligand[,3], Receptor = res.receptor[,3])
+rownames(ligand.receptor.data) <- res.ligand$cellType
+row_split <- res.ligand$type
+Heatmap(ligand.receptor.data, name = "Number of Interaction", col = colorRamp2(c(0,1200), c("white", "red")), show_row_dend = F, show_column_dend = F, row_split = row_split, 
+        width = unit(2, "cm"), height = unit(10, "cm"), cell_fun = function(j, i, x, y, width, height, fill) {
+        grid.text(sprintf("%.0f", ligand.receptor.data[i, j]), x, y, gp = gpar(fontsize = 7))}
+        )
+dev.off()
+
+# heatmap shows the conditions of ligand and receptor
+All.interaction.number <- sapply(sig.interaction.cellType, function(x){
+    return(nrow(x$All))
+})
+All.interaction.number <- data.frame(cellType = names(All.interaction.number), interaction.num = All.interaction.number)
+
+class.type <- rep("Lymphoid", nrow(All.interaction.number))
+class.type[c(1, 8, 13:16, 19)] <- "Myeloid"
+class.type[which(All.interaction.number$cellType %in% c("Endothelium VCAM1+", "Endothelium VCAM1-", "Mesangial cell"))] <- "Other"
+class.type[which(All.interaction.number$cellType %in% c("Tumor"))] <- "Tumor"
+All.interaction.number$class <- factor(class.type, levels = c("Lymphoid", "Myeloid", "Tumor", "Other"))
+pdf("sig.interaction.number.pdf")
+ggbarplot(All.interaction.number, x="cellType", y="interaction.num", xlab = "", ylab = "Number of interactions", fill = "class", color = "white", palette = colors,sort.val = "desc",sort.by.groups = FALSE,x.text.angle = 60)
+dev.off()
+saveRDS(All.interaction.number, "All.interaction.number.rds")
+
+#### Prepare circular plots for all cell types
+track.matrix <- matrix(c(rep(0, nrow(All.interaction.number)), All.interaction.number$interaction.num), ncol = 2)
+rownames(track.matrix) <- All.interaction.number[,1]
+track.matrix <- track.matrix[order.value,]
+
+#function --- plot circular
+plot.circular <- function(track.matrix, gene.pair, mypvals, mymeans, data.merge = data.merge, idx = c(1, 3:10, 12), means.threshold = 1, pval.threshold = 0.01, risk.table = T, High.low = F, colors = c("#D95319", "#F39B7F", "#3B6793","#4285F4")){
+    gene.a <- gene.pair[1]
+    gene.b <- gene.pair[2]
+    mymeans %>% dplyr::filter(gene_a == gene.a & gene_b == gene.b) %>% reshape2::melt() -> meansdf
+    meansdf <- meansdf[, -idx]
+    colnames(meansdf)<- c("interacting_pair", "ligand", "CC","means")
+    mypvals %>% dplyr::filter(gene_a == gene.a & gene_b == gene.b) %>% reshape2::melt() -> pvalsdf
+    pvalsdf <- pvalsdf[, -idx]
+    colnames(pvalsdf)<- c("interacting_pair", "ligand", "CC","means")
+    pvalsdf$joinlab<- paste0(pvalsdf$interacting_pair, "_", pvalsdf$CC)
+    meansdf$joinlab<- paste0(meansdf$interacting_pair, "_", meansdf$CC)
+    pldf <- merge(pvalsdf, meansdf, by = "joinlab")
+    pldf <- pldf[,-c(1,6:8)]
+    colnames(pldf) <- c("interacting_pair", "ligand", "CC", "pvals", "means")
+    sig.interaction <- pldf%>% filter(means >= means.threshold & pvals < pval.threshold)
+
+    interaction.matrix <- apply(sig.interaction, 1, function(x){
+        if(x[2] == "Receptor"){
+            a <- unlist(strsplit(x[3], "\\|"))
+            b <- c(a[2], a[1], x[5])
+        }else{
+            a <- unlist(strsplit(x[3], "\\|"))
+            b <- c(a[1], a[2], x[5])
+        }
+        return(b)
+    })
+    interaction.matrix <- as.data.frame(t(interaction.matrix))
+    interaction.matrix$means <- as.numeric(interaction.matrix$means)
+    p <- Plot.CircularChart(track.matrix = track.matrix, interaction.matrix = interaction.matrix, 
+                            bg.col = NULL, cex = 0.8,
+                            link.col = NULL, arr.lwd = 0.3)
+    print(p)
+
+    survival.res <- analysis.diff.survival.TCGA(interest.gene = c(gene.a, gene.b), diff.gene.pro = DESeq2.result, exp.data.process = DESeq2.normalized_counts, clin.data = clin.data, EnhancedVolcano.plot = F, Box.plot = T, main = paste0(gene.a, "-", gene.b), meta.signature = F)
+    combined.res <- survival.combined(geneA = gene.a, geneB = gene.b, clin.data = clin.data, DESeq2.result = DESeq2.result, DESeq2.normalized_counts = DESeq2.normalized_counts, risk.table = risk.table, High.low = High.low, colors = colors)
+
+    p <- FeaturePlot(data.merge, features = gene.a, cols = c("lightgrey", "red"))
+    print(p)
+    p <- FeaturePlot(data.merge, features = gene.b, cols = c("lightgrey", "red"))
+    print(p)
+    p <- VlnPlot(data.merge, features = gene.a, pt.size = 0.5) + NoLegend()
+    print(p)
+    p <- VlnPlot(data.merge, features = gene.b, pt.size = 0.5) + NoLegend()
+    print(p)
+    return(interaction.matrix)
+}
+
+#########################################################analysis tumor cell################################################
+tumor <- sig.interaction.cellType[["Tumor"]]
+#Statistics and tumor interaction cell type distribution
+#ligand
+ligand.interact.tumor <- apply(tumor$Ligand, 1, function(x){
+    a <- unlist(strsplit(x[4], "\\|"))
+    return(a)
+})
+ligand.interact.tumor.freq <- data.frame(table(as.character(ligand.interact.tumor)))
+
+idx <- which(ligand.interact.tumor.freq$Var1 == "Tumor")
+ligand.interact.tumor.freq[idx, 2] <- ligand.interact.tumor.freq[idx, 2]-nrow(tumor$Ligand)
+
+receptor.interact.tumor <- apply(tumor$Receptor, 1, function(x){
+    a <- unlist(strsplit(x[4], "\\|"))
+    return(a)
+})
+receptor.interact.tumor.freq <- data.frame(table(as.character(receptor.interact.tumor)))
+idx <- which(receptor.interact.tumor.freq$Var1 == "Tumor")
+receptor.interact.tumor.freq[idx, 2] <- receptor.interact.tumor.freq[idx, 2]-nrow(tumor$Receptor)
+
+type <- rep("Lymphoid", nrow(receptor.interact.tumor.freq))
+type[c(7,10,12:13,16:18)] <- "Myeloid"
+type[20] <- "Tumor"
+type[c(8:9, 11)] <- "Other"
+receptor.interact.tumor.freq$Type <- factor(type, levels = c("Lymphoid", "Myeloid", "Tumor", "Other"))
+ligand.interact.tumor.freq$Type <- factor(type, levels = c("Lymphoid", "Myeloid", "Tumor", "Other"))
+
+ligand.interact.tumor.freq$Var1 <- factor(ligand.interact.tumor.freq$Var1, levels = order.value)
+receptor.interact.tumor.freq$Var1 <- factor(receptor.interact.tumor.freq$Var1, levels = order.value)
+pdf("Tumor/tumor.ligand.receptor.pdf")
+ggbarplot(ligand.interact.tumor.freq, x="Var1", y="Freq", color = "Type", xlab = "", palette = colors, ylab = "Number of Ligand", fill = "Type", sort.val = "desc", sort.by.groups = FALSE, x.text.angle = 60) 
+ggbarplot(receptor.interact.tumor.freq, x="Var1", y="Freq", color = "Type", xlab = "", palette = colors, ylab = "Number of Receptor", fill = "Type", sort.val = "desc", sort.by.groups = FALSE, x.text.angle = 60)
+ggbarplot(ligand.interact.tumor.freq[-c(8,9,11),], x="Var1", y="Freq", color = "Type", xlab = "", palette = colors[1:3], ylab = "Number of Ligand", fill = "Type", sort.val = "desc", sort.by.groups = FALSE, x.text.angle = 60) 
+ggbarplot(receptor.interact.tumor.freq[-c(8,9,11),], x="Var1", y="Freq", color = "Type", xlab = "", palette = colors[1:3], ylab = "Number of Receptor", fill = "Type", sort.val = "desc", sort.by.groups = FALSE, x.text.angle = 60)
+
+ligand.receptor.data <- data.frame(Ligand = ligand.interact.tumor.freq[,2], Receptor = receptor.interact.tumor.freq[,2])
+rownames(ligand.receptor.data) <- ligand.interact.tumor.freq$Var1
+row_split <- factor(ligand.interact.tumor.freq$Type, levels = c("Lymphoid", "Myeloid", "Tumor", "Other"))
+Heatmap(ligand.receptor.data, name = "Number of Interaction", col = colorRamp2(c(0,50), c("white", "red")), show_row_dend = F, show_column_dend = F, row_split = row_split, 
+        width = unit(2, "cm"), height = unit(8, "cm"), cell_fun = function(j, i, x, y, width, height, fill) {
+        grid.text(sprintf("%.0f", ligand.receptor.data[i, j]), x, y, gp = gpar(fontsize = 7))}
+        )
+
+# total number
+total.interaction <- data.frame(cellType = rownames(ligand.receptor.data), number = ligand.receptor.data[,1]+ligand.receptor.data[,2], type = ligand.interact.tumor.freq$Type)
+total.interaction$type <- factor(total.interaction$type, levels = c("Lymphoid", "Myeloid", "Tumor", "Other"))
+ggbarplot(total.interaction, x="cellType", y="number", color = "type", xlab = "", palette = colors, ylab = "Number of interaction", fill = "type", sort.val = "desc", sort.by.groups = FALSE, x.text.angle = 60) 
+dev.off()
+saveRDS(total.interaction, file = "Tumor/tumor.interaction.number.rds")
+
+#### Interaction between cancer cells and other cells
+#Screening significant interaction pairs-ligand
+selected_rows <- unique(tumor$All$interacting_pair[which(tumor$All$means>=mean.threshold & tumor$All$pvals<pval.threshold)])
+selected_columns <- sort(unique(as.character(tumor$All$CC)))
+pdf("Tumor/tumor.sig.interaction.pdf")
+p1 <- dot_plot(selected_rows = selected_rows, selected_columns = selected_columns, 
+               breaks = c(-2, -1, 0, 1, 2), limit = c(-2, 2), legend.key.size = unit(0.1, "inches"),
+               size = 8, max.dot.size = 3, means_path = "Renal/mymeans.usr.txt", pvalues_path = "Renal/mypvals.usr.txt")
+print(p1)
+dev.off()
+
+pdf("Tumor/VEGFA.expression.pdf", width = 3, height = 3)
+p <- DotPlot(data.merge.new, features = c("VEGFA", "NRP1", "KDR", "GPC1"), cols = c("#1e90ff", "#F15F30"), dot.scale = 3.5) + xlab("") + ylab("") + theme(axis.text.x = element_text(size = 8, angle = 90)) + theme(axis.text.y = element_text(size = 8)) + NoLegend()
+print(p)
+dev.off()
+pdf("Tumor/VEGFA.expression-test.pdf", width = 3, height = 3)
+p <- DotPlot(data.merge.new, features = c("VEGFA", "NRP1", "KDR", "GPC1"), cols = c("#1e90ff", "#F15F30"), dot.scale = 3.5) + xlab("") + ylab("") + theme(axis.text.x = element_text(size = 8, angle = 90)) + theme(axis.text.y = element_text(size = 8))
+print(p)
+dev.off()
+pdf("Tumor/interested.interaction.pairs.pdf")
+p <- VlnPlot(data.merge.new, features = c("MIF", "CD74", "CXCR4"), pt.size = 0, ncol = 2) + NoLegend()
+print(p)
+p <- VlnPlot(data.merge.new, features = c("VEGFA", "NRP1", "KDR", "GPC1"), pt.size = 0, ncol = 2) + NoLegend()
+print(p)
+p <- FeaturePlot(data.merge, features = c("VEGFA", "NRP1", "KDR", "GPC1"), cols = c("lightgrey", "red"))
+print(p)
+p <- VlnPlot(data.merge.new, features = c("CXCL2", "DPP4", "CD24", "SIGLEC10"), pt.size = 0, ncol = 2)
+print(p)
+p <- VlnPlot(data.merge.new, features = c("APOE", "APOC1", "TREM2", "VLDLR"), pt.size = 0, ncol = 2)
+print(p)
+p <- VlnPlot(data.merge.new, features = c("CD70", "CD27", "CALM1", "HMMR"), pt.size = 0, ncol = 2)
+print(p)
+dev.off()
+
+features <- c("MIF",  "VEGFA", "GPC1", "CD24",
+              "DPP4", "CD70", "CALM1", "APOE", "VLDLR") 
+pdf("Tumor/tumor.sig.ligand.receptor.survival.pdf")
+res <- analysis.diff.survival.TCGA(interest.gene = features, diff.gene.pro = DESeq2.result, exp.data.process = DESeq2.normalized_counts, clin.data = clin.data, EnhancedVolcano.plot = F, Box.plot = T, main = "cancer.sig.ligand.receptor", meta.signature = F)
+features <- c("MIF",  "VEGFA", "GPC1", 
+              "DPP4", "CD70", "CALM1", "APOE", "VLDLR")
+features.list <- as.list(features)
+names(features.list) <- features
+cox.res <- RCC.icb.analysis(signature.list = features.list, expresionMatrix = normalized_expression, clincal.info = patient.info.RNA)
+dev.off()
+
+
+#### Prepare circular plots 
+track.matrix <- matrix(c(rep(0, nrow(total.interaction)), total.interaction$number), ncol = 2)
+rownames(track.matrix) <- total.interaction[,1]
+track.matrix <- track.matrix[order.value,]
+
+###########################Tumor VS TAM
+tumor_TAM <- tumor$All[grep("TAM", tumor$All$CC),]
+selected_rows <- unique(tumor_TAM$interacting_pair[which(tumor_TAM$means>=mean.threshold & tumor_TAM$pvals<pval.threshold)])
+
+# features <- as.character(unlist(strsplit(selected_rows, "_")))
+# pdf("Tumor/Myeloid/tumor.sig.ligand.receptor.survival.pdf")
+# res <- analysis.diff.survival.TCGA(interest.gene = features, diff.gene.pro = DESeq2.result, exp.data.process = DESeq2.normalized_counts, clin.data = clin.data, EnhancedVolcano.plot = F, Box.plot = T, main = "Tumor.TAM", meta.signature = F)
+# dev.off()
+# write.csv(res,file = "Tumor/Myeloid/tumor.sig.ligand.receptor.survival.csv", row.names = F)
+
+TAM.DEGs <- readRDS("/data/active_data/lzl/RenalTumor-20200713/DataAnalysis-20210803/scRNA/5.Immune/Macrophage/cellType.DE.pro.rds")
+Tumor.DEGs <- readRDS("/data/active_data/lzl/RenalTumor-20200713/DataAnalysis-20210803/scRNA/2.Cluster/AnnotateCellType/cellType.sig.pos.DEGs.rds")
+Tumor.DEGs <- Tumor.DEGs %>% filter(cluster == "Tumor") %>% filter(p_val_adj<0.05 & avg_log2FC>0.25)
+TAM.DEGs <- TAM.DEGs %>% filter(p_val_adj<0.05 & avg_log2FC>0.25)
+DEGs <- unique(c(Tumor.DEGs$gene, TAM.DEGs$gene))
+interaction.genes <- unlist(str_split(selected_rows, pattern="_"))
+sig.DEGs <- intersect(interaction.genes, DEGs)
+
+interaction.info <- sapply(selected_rows, function(x){
+    a <- unlist(str_split(x, "_"))
+    first <- match(a[1], sig.DEGs)
+    second <- match(a[2], sig.DEGs)
+    if(!is.na(first)&!is.na(second)){
+        return("ok")
+    }else{
+        return("no")
+    }
+})
+
+idx <- which(interaction.info == "ok")
+selected_rows <- selected_rows[idx]
+a <- apply(tumor_TAM, 1, function(x){
+    b <- unlist(str_split(as.character(x[4]), "\\|"))
+    if((x[3]=="Ligand" & b[1]=="Tumor")|(x[3]=="Receptor" & b[2]=="Tumor")){
+        return("ok")
+    }else{
+        return("error")
+    }
+})
+selected_rows1 <- selected_rows[which(selected_rows %in% tumor_TAM$interacting_pair[which(a=="error")])]
+selected_columns1 <- as.character(tumor_TAM$CC[which(tumor_TAM$interacting_pair %in% selected_rows1)])
+selected_rows2 <- selected_rows[which(selected_rows %in% tumor_TAM$interacting_pair[which(a=="ok")])]
+selected_columns2 <- as.character(tumor_TAM$CC[which(tumor_TAM$interacting_pair %in% selected_rows2)])
+
+pdf("Tumor/Myeloid/tumor.TAM.interaction.receptor.pdf", width = 3.8, height = 4.2)
+p1 <- dot_plot(selected_rows = selected_rows1, selected_columns = unique(selected_columns1), 
+               breaks = c(-2, -1, 0, 1, 2), limit = c(-2, 2), legend.key.size = unit(0.1, "inches"), order.column = c("TAM-C1QB|Tumor", "TAM-RGCC|Tumor", "TAM-LGALS3|Tumor"),
+               size = 8, max.dot.size = 4, means_path = "Renal/mymeans.usr.txt", pvalues_path = "Renal/mypvals.usr.txt")
+print(p1)
+dev.off()
+pdf("Tumor/Myeloid/tumor.TAM.interaction.ligand.pdf", width = 3.8, height = 4.2)
+p1 <- dot_plot(selected_rows = selected_rows2, selected_columns = unique(selected_columns2), 
+               breaks = c(-2, -1, 0, 1, 2), limit = c(-2, 2), legend.key.size = unit(0.1, "inches"), order.column = c("Tumor|TAM-C1QB", "Tumor|TAM-RGCC", "Tumor|TAM-LGALS3"),
+               size = 8, max.dot.size = 4, means_path = "Renal/mymeans.usr.txt", pvalues_path = "Renal/mypvals.usr.txt")
+print(p1)
+dev.off()
+
+pdf("Tumor/Myeloid/SPP1.expression.pdf", width = 3, height = 3)
+p <- DotPlot(data.merge.new, features = c("SPP1", "PTGER4", "ITGA9", "ITGA4"), cols = c("#1e90ff", "#F15F30"), dot.scale = 3.5) + xlab("") + ylab("") + theme(axis.text.x = element_text(size = 8, angle = 90)) + theme(axis.text.y = element_text(size = 8)) + NoLegend()
+print(p)
+dev.off()
+pdf("Tumor/Myeloid/SPP1.expression-test.pdf", width = 3, height = 3)
+p <- DotPlot(data.merge.new, features = c("SPP1", "PTGER4", "ITGA9", "ITGA4"), cols = c("#1e90ff", "#F15F30"), dot.scale = 3.5) + xlab("") + ylab("") + theme(axis.text.x = element_text(size = 8, angle = 90)) + theme(axis.text.y = element_text(size = 8))
+print(p)
+dev.off()
+pdf("Tumor/Myeloid/interested.interaction.pairs.pdf")
+p <- VlnPlot(data.merge.new, features = c("MIF", "CD74", "CD44", "SPP1"), pt.size = 0, ncol = 2) + NoLegend()
+print(p)
+p <- VlnPlot(data.merge.new, features = c("CD24", "SIGLEC10", "APOE", "TREM2"), pt.size = 0, ncol = 2)
+print(p)
+p <- VlnPlot(data.merge.new, features = c("VIM", "CD44", "RPS19", "C5AR1"), pt.size = 0, ncol = 2)
+print(p)
+p <- VlnPlot(data.merge.new, features = c("PKM", "CD44", "CXCL8", "SDC2"), pt.size = 0, ncol = 2)
+print(p)
+p <- VlnPlot(data.merge.new, features = c("FN1", "PLAUR", "CXCL8", "SDC2"), pt.size = 0, ncol = 2)
+print(p)
+dev.off()
+
+features <- c("MIF",  "VEGFA", "GPC1", 
+              "DPP4", "CD70", "CALM1", "APOE", "VLDLR") 
+features.list <- as.list(features)
+names(features.list) <- features
+cox.res <- RCC.icb.analysis(signature.list = features.list, expresionMatrix = normalized_expression, clincal.info = patient.info.RNA)
+dev.off()
+
+###########################Tumor VS T cell
+tumor_lymphoid <- tumor$All[grep("Proliferative T cell|Treg|Exhausted", tumor$All$CC),]
+selected_rows <- unique(tumor_lymphoid$interacting_pair[which(tumor_lymphoid$means>=mean.threshold & tumor_lymphoid$pvals<pval.threshold)])
+
+DEGs <- readRDS("/data/active_data/lzl/RenalTumor-20200713/DataAnalysis-20210803/scRNA/5.Immune/CD8T/cellType.DE.rds")
+T.DEGs <- DEGs[["Exhausted"]]  %>% filter(p_val_adj<0.05 & avg_log2FC>0.25)
+all.DEGs <- readRDS("/data/active_data/lzl/RenalTumor-20200713/DataAnalysis-20210803/scRNA/2.Cluster/AnnotateCellType/cellType.sig.pos.DEGs.rds")
+other.DEGs <- all.DEGs %>% filter(cluster %in% c("Proliferative T cell", "Treg", "Tumor")) %>% filter(p_val_adj<0.05 & avg_log2FC>0.25)
+DEGs <- unique(c(T.DEGs$gene, other.DEGs$gene))
+interaction.genes <- unlist(str_split(selected_rows, pattern="_"))
+sig.DEGs <- intersect(interaction.genes, DEGs)
+
+interaction.info <- sapply(selected_rows, function(x){
+    a <- unlist(str_split(x, "_"))
+    first <- match(a[1], sig.DEGs)
+    second <- match(a[2], sig.DEGs)
+    if(!is.na(first)&!is.na(second)){
+        return("ok")
+    }else{
+        return("no")
+    }
+})
+idx <- which(interaction.info == "ok")
+selected_rows <- selected_rows[idx]
+a <- apply(tumor_lymphoid, 1, function(x){
+    b <- unlist(str_split(as.character(x[4]), "\\|"))
+    if((x[3]=="Ligand" & b[1]=="Tumor")|(x[3]=="Receptor" & b[2]=="Tumor")){
+        return("ok")
+    }else{
+        return("error")
+    }
+})
+selected_rows1 <- selected_rows[which(selected_rows %in% tumor_lymphoid$interacting_pair[which(a=="error")])]
+selected_columns1 <- as.character(tumor_lymphoid$CC[which(tumor_lymphoid$interacting_pair %in% selected_rows1)])
+selected_rows2 <- selected_rows[which(selected_rows %in% tumor_lymphoid$interacting_pair[which(a=="ok")])]
+selected_columns2 <- as.character(tumor_lymphoid$CC[which(tumor_lymphoid$interacting_pair %in% selected_rows2)])
+
+pdf("Tumor/Lymphoid/tumor.Lymphoid.interaction.receptor.pdf", width = 3.8, height = 4)
+p1 <- dot_plot(selected_rows = selected_rows1, selected_columns = unique(selected_columns1), 
+               breaks = c(-2, -1, 0, 1, 2), limit = c(-2, 2), legend.key.size = unit(0.1, "inches"), order.column = c("CD8+ T-Exhausted|Tumor", "Proliferative T cell|Tumor", "Treg|Tumor"),
+               size = 8, max.dot.size = 4, means_path = "Renal/mymeans.usr.txt", pvalues_path = "Renal/mypvals.usr.txt")
+print(p1)
+dev.off()
+pdf("Tumor/Lymphoid/tumor.Lymphoid.interaction.ligand.pdf", width = 3.65, height = 4)
+p1 <- dot_plot(selected_rows = selected_rows2, selected_columns = unique(selected_columns2), 
+               breaks = c(-2, -1, 0, 1, 2), limit = c(-2, 2), legend.key.size = unit(0.1, "inches"), order.column = c("Tumor|CD8+ T-Exhausted", "Tumor|Proliferative T cell", "Tumor|Treg"),
+               size = 8, max.dot.size = 4, means_path = "Renal/mymeans.usr.txt", pvalues_path = "Renal/mypvals.usr.txt")
+print(p1)
+dev.off()
+
+pdf("Tumor/Lymphoid/interested.interaction.pairs.pdf")
+p <- VlnPlot(data.merge.new, features = c("LTB", "LTBR", "CD27", "CD70"), pt.size = 0, ncol = 2) + NoLegend()
+print(p)
+dev.off()