R/LRC3_Alpha.R

In XiaojuanShen/LRC3: Ligand-Receptor based Communication

#' Prepare the data
#'
#'This function is to prepare the data.
#'
#' @param dd The input expresstion matrix. The format of the data should be a dataframe or a matrix,
#' with rownames as the gene names, and colomn names as the cell type names.
#' It could also be an SingleCellExperiment objects with thecell types names stored in the “cell_type1” slot of the colData.
#' @param SingleCellExperimentObject if \code{SingleCellExperimentObject = TRUE},the imput data format is SingleCellExperiment object,
#' if \code{SingleCellExperimentObject = FALSE},the imput data format is an expression matrix.
#' @param log_counts parameter for 'SingleCellExperimentObject', if \code{log_counts = FALSE}, use none log normalized data,
#' if \code{log_counts = TRUE},use log normalized data.
#' @param SparseMatrix if \code{SparseMatrix = TRUE}, the SingleCellExperimentObject input data is a sparse matrix.
#' @param Species choosing a species of the dataset
#' @return data_list an object that contain the expression matrix dd, the cell type name index, and cell type names.
#' @export
#'
PrepareData <- function(dd = dd,SingleCellExperimentObject = FALSE,log_counts = FALSE,SparseMatrix = FALSE, Species = c("Human","Mouse")){
  if (SingleCellExperimentObject == TRUE){

    if (is.null(dd$cell_type1)){stop("Error! there is no cell_type1 slot! cell_type1 slot in the SingleCellExperimentObject is required")}

    CellTypes <- dd$cell_type1
    # CellTypes <- dd$ClusterID#this is for the MCA data
    if (log_counts == TRUE){
      dd <- logcounts(dd)
      # dd <- exprs(dd)
    } else if (log_counts == FALSE){
      dd <- counts(dd)
      # dd <- dd[,-IndexNA]
    }
    if (length(which(duplicated(rownames(dd)))) != 0){
      dd <- dd[-duplicated(rownames(dd)),]
    }
    Species <- match.arg(Species)
    if (Species == "Human"){
      LR_pairs <- LRpairs
    } else if (Species == "Mouse"){
      LR_pairs <- Mouse_LRpairs
    }

    Ligand <- as.character(unique(LR_pairs$Ligand.ApprovedSymbol))
    Receptor <- as.character(unique(LR_pairs$Receptor.ApprovedSymbol))
    LigandRecptor <- c(Ligand,Receptor)
    dd <- dd[which(rownames(dd) %in% LigandRecptor),]

    if (SparseMatrix == TRUE){
      dd <- as.matrix(dd)
    }
    colnames(dd) <- CellTypes
  }

  if (SingleCellExperimentObject == FALSE){
    if (length(which(duplicated(rownames(dd)))) != 0){
      dd <- dd[-duplicated(rownames(dd)),]
    }
    Species <- match.arg(Species)
    if (Species == "Human"){
      LR_pairs <- LRpairs
    } else if (Species == "Mouse"){
      LR_pairs <- Mouse_LRpairs
    }

    Ligand <- as.character(unique(LR_pairs$Ligand.ApprovedSymbol))
    Receptor <- as.character(unique(LR_pairs$Receptor.ApprovedSymbol))
    LigandRecptor <- c(Ligand,Receptor)
    dd <- dd[which(rownames(dd) %in% LigandRecptor),]
  }

    if (length(unique(colnames(dd))) == 1){
      print("Warning! There is only one cell type in this data!")
      # break
    }
  CellTypeNames <- unique(colnames(dd))
  CellTypeNames <- gsub("-","_",CellTypeNames)
  CellTypeNames <- sort(CellTypeNames)
  if(length(which(rowSums(dd)== 0)) !=0){
    dd <- dd[-which(rowSums(dd) == 0),]

  }
  if (length(which(duplicated(rownames(dd)))) != 0){
    dd <- dd[-which(duplicated(rownames(dd))),]
  }
  Number_CellType <- c(1:length(CellTypeNames))
  PreColumnNames <- colnames(dd)
  ChangedColumnNames <- gsub("-","_",PreColumnNames)
  Colnames_Data <- as.list(ChangedColumnNames)
  Num_Colnames <- sapply(Colnames_Data,function(x){which(CellTypeNames %in% x)})
  colnames(dd) <- Num_Colnames
  data_list <- list()
  data_list[[1]] <- dd
  data_list[[2]] <- Num_Colnames
  data_list[[3]] <- CellTypeNames
  return(data_list)
}



#' Get a list of results.
#'
#' This function Performs calculations on the communication between different cell types and generate a list of results.The
#' result list contains 13 matrices or vectors with each gives one type of information on the input data.
#' The name of the 13 elements are contact_number_Remaining, scaled_contact_number_remaining,
#' Remaining_contactINF, Remaining_scaled_contactINF, rem_LRpairs,
#' LigandExpressionLevelINF, ReceptorExpressionLevelINF,
#' ScaledLigandExpressionLevelINF, ScaledReceptorExpressionLevelINF,
#' p, k, GroupedUniqueLRMatrix, GroupedNum_Colnames.
#' @param data_list the input data generated by function PrepareData.
#' @param PvalAsThreshold parameter used to set the Pval as the threshold.
#' if \code{PvalAsThreshold = TRUE},the algrithm will perform a random permutation process on the expression matrix,
#' and get a threshold from the permutation for each of the gene.
#' if \code{PvalAsThreshold = FALSE},the algrithm will use one fixed threshold for all the genes.The defalt is TRUE.
#' @param FixedThreshold threshold when using a fixed value instead of Pval.
#' Uses should adust this threshold as needed.The defalt is set as FixedThreshold = 1.
#' @param MultiThreshold Boolean value to set weather use multiple threshold or only one threshold
#' @param PercentForPval threshold of the Pval. After the permutation on the expression matrix,
#' the Pval for each of the gene is chosen by this threshold from the permutation procedure.
#'  The default is 0.95, which means the the value of the top 5 percent.
#' @param PercentForPtest threshold for the P value test, the default is 0.95.
#' @param DefaultLRlist if \code{DefaultLRlist = TRUE}, the a defult list containning 2570 LR pairs will be used.
#' if \code{DefaultLRlist = FALSE}, users need to define their own parameter for of 'LR_pairs'.
#' @param Species select the species that users are analyzing. LRC3 could analyzing both human and mouse datasets.
#' @param LR_pairs the imput LRpairs users have to provide when setting the \code{DefaultLRlist = FALSE}
#' The user can also choose the Membrane or NonMembrane LR pairs for both Mouse and Human that curated and included in this package,
#' by setting \code{LR_pairs = Human_Membrane} or \code{LR_pairs = Human_NonMembrane} for Human,
#' and  \code{LR_pairs = Mouse_Membrane}, \code{LR_pairs = Mouse_NonMembrane} for Mouse.
#' @param ncore the number of cores to be used
#' @importFrom parallel makeCluster stopCluster
#' @importFrom doParallel registerDoParallel
#' @importFrom parallel detectCores
#' @importFrom foreach foreach
#' @importFrom foreach %dopar%
#' @importFrom abind abind
#' @return LRC3_list a list that contains 13 elements.
#' 1st element is the contact number between different cell types.
#' 2nd element is the scaled contact number between different cell types.
#' 3rd element is an array containing the information of ligand-receptor pairs used between different cell types.
#' 4th element is an array containing the scaled information of ligand-receptor pairs used between different cell types.
#' 5th element is the ligand-receptor names that are used in the input dataset.
#' 6th element is the expression level of ligand that are used in each cell type.
#' 7th element is the expression level of receptor that are used in each cell type.
#' 8th element is the scaled expression level of ligand that are used in each cell type.
#' 9th element is the scaled the expression level of receptor that are used in each cell type.
#' 10th element is the threshold for ligand.
#' 11th element is the threshold for receptor.
#' 12th element is the expression matrix ordered by cell type.
#' 13th element is the index of cell types.
#' 14th element is the contact number between each cell type after testing the significance of the total contact number
#' @export

LRC3_INF <- function(data_list,PvalAsThreshold = TRUE, FixedThreshold = 1,MultiThreshold = FALSE,PercentForPval = 0.95,PercentForPtest = 0.95,DefaultLRlist = TRUE,Species = c("Human","Mouse"),LR_pairs = LR_pairs,ncore = 2){

  LRmatrix <- function(data_list,DefaultLRlist = TRUE,Species = c("Human","Mouse")){
    dd <- data_list[[1]]
    rnames <- rownames(dd)

    if (DefaultLRlist == TRUE){
      Species <- match.arg(Species)
      if (Species == "Human"){
        LR_pairs <- LRpairs
      }else if (Species == "Mouse") {
        LR_pairs <- Mouse_LRpairs
      }

    } else if (DefaultLRlist == FALSE){
      LR_pairs <- as.data.frame(LR_pairs)
      names(LR_pairs) <- c("Ligand.ApprovedSymbol","Receptor.ApprovedSymbol")
    }

    # LR_pairs <- read.csv("/Users/xs2/R/Genome_biology/Ligand-receptor-pairs.csv",header = T)
    # LR_pairs <- read.csv("/Users/xs2/R/Genome_biology/LRpairs_Nature.csv")

    # LR_pairs <- LRpairs[c1 <- c(895,1773,1785,2531, 682,2201,2494,2515,2528,2533,180,181,182,
    #                             183,184,185,186,187,1036,1037,1038,1039,1922,1923,1924,1925,1926,1927,
    #                             656,1421,1422,1423,1424,1450,432,2250,2456,1446,1447,1450,1451,1453,1454,
    #                             1455,844,845,846,847,848,849,850,851,2356,2357,2358,2359,2360,
    #                             2504,2505,2506,2507,2508,2509,2521,2522,2523,2524,2525,2526,2527,2528,2529,2530,2531,2532,2533,2534,2535,2536,
    #                             2537,2538,2539,2540,2541,2542,2543,2544,1223,2356,2357,2358,2359,2360,2361,2362,872,873,874,875,876,877),]
    # LR_pairs <- LR_pairs[-which(duplicated(LR_pairs)),]
    # LR_pairs <- read.csv("/Users/xs2/R/Genome_biology/DEC_LRpairs.csv")
    Num_Colnames <- data_list[[2]]


    ligand_list <- as.character(unique(LR_pairs$Ligand.ApprovedSymbol))
    receptor_list <- as.character(unique(LR_pairs$Receptor.ApprovedSymbol))
    rem_LRpairs <- LR_pairs[which(LR_pairs$Ligand.ApprovedSymbol %in% rnames),]#[1] 682   3
    rem_LRpairs <- rem_LRpairs[which(rem_LRpairs$Receptor.ApprovedSymbol %in% rnames),]#[1] 301   3

    rem_L <- rem_LRpairs$Ligand.ApprovedSymbol
    rem_R <- rem_LRpairs$Receptor.ApprovedSymbol

    if (length(rem_L) <= 1 | length(rem_R) <=1){
      L_dd <- 0
      R_dd <- 0
      rem_LRpairs <- 0
      UniqueLRMatrix <- 0
      GroupedUniqueLRMatrix <- 0
      GroupedNum_Colnames <- 0

    } else if (length(rem_L) > 1 & length(rem_R) >1){
      rem_unique_LR <- unique(c(as.character(rem_L),as.character(rem_R)))
      UniqueLRMatrix <- dd[which(rownames(dd) %in% rem_unique_LR),]
      GroupedUniqueLRMatrix <- dd[,which(Num_Colnames == 1)]
      GroupedNum_Colnames <- Num_Colnames[which(Num_Colnames == 1)]
      for (i in 2:length(unique(colnames(dd)))){
        TEM_Matrix <- dd[,which(Num_Colnames == i)]
        GroupedNum_Colnames <- c(GroupedNum_Colnames,Num_Colnames[which(Num_Colnames == i)])
        GroupedUniqueLRMatrix <- cbind(GroupedUniqueLRMatrix,TEM_Matrix)
      }

      L_dd <- data.frame(matrix(nrow = length(rem_L),ncol = length(colnames(dd))))
      R_dd <- data.frame(matrix(nrow = length(rem_L),ncol = length(colnames(dd))))
      Index_A <- rep(0,length(rem_L))
      Index_B <- rep(0,length(rem_L))
      for (i in 1:length(rem_L)){
        Index_A[i] <- which(rownames(dd) == rem_L[i])
        Index_B[i] <- which(rownames(dd) == rem_R[i])
      }
      L_dd <- dd[Index_A,]
      R_dd <- dd[Index_B,]
    }

    # which(rowSums(L_dd)>0)

    LRmatrix_list <- vector("list")
    LRmatrix_list[[1]] <- L_dd
    LRmatrix_list[[2]] <- R_dd
    LRmatrix_list[[3]] <- rem_LRpairs
    LRmatrix_list[[4]] <- UniqueLRMatrix
    LRmatrix_list[[5]] <- GroupedUniqueLRMatrix
    LRmatrix_list[[6]] <- GroupedNum_Colnames
    names(LRmatrix_list) <- c("L_dd","R_dd","rem_LRpairs","UniqueLRMatrix","GroupedUniqueLRMatrix","GroupedNum_Colnames")
    return(LRmatrix_list)
  }

  LR_mean <- function(LRmatrix_list,data_list){
    L_dd <- LRmatrix_list[[1]]
    R_dd <- LRmatrix_list[[2]]
    UniqueLRMatrix <- LRmatrix_list[[4]]
    GroupedUniqueLRMatrix <- LRmatrix_list[[5]]
    dd <- data_list[[1]]
    GroupedNum_Colnames <- LRmatrix_list[[6]]
    Num_cell_type <- length(unique(colnames(dd)))

    L_mean <- NULL
    tem_mean <- NULL
    UniqueLR_mean <- NULL

    for (i in 1:Num_cell_type){
      if (length(L_dd[1, colnames(L_dd) == i]) != 1) {
        tem_mean <- rowMeans(L_dd[, colnames(L_dd) == i])
        L_mean <- cbind(L_mean, tem_mean)
        tem_UniqueLRmean <- rowMeans(UniqueLRMatrix[, colnames(L_dd) == i])
        UniqueLR_mean <- cbind(UniqueLR_mean,tem_UniqueLRmean)
      } else if (length(L_dd[1, colnames(L_dd) == i]) == 1) {
        tem_mean <- L_dd[, colnames(L_dd) == i]
        L_mean <- cbind(L_mean, tem_mean)
        tem_UniqueLRmean <- UniqueLRMatrix[, colnames(L_dd) == i]
        UniqueLR_mean <- cbind(UniqueLR_mean,tem_UniqueLRmean)
      }
    }
    colnames(L_mean) <- 1:Num_cell_type
    colnames(UniqueLR_mean) <- 1:Num_cell_type


    R_mean <- NULL
    tem_mean <- NULL
    for (i in 1:Num_cell_type){
      if (length(L_dd[1,colnames(R_dd) == i]) != 1){
        tem_mean <- rowMeans(R_dd[, colnames(R_dd) == i])
        R_mean <- cbind(R_mean, tem_mean)
      }else {
        tem_mean <- R_dd[,colnames(R_dd) == i]
        R_mean <- cbind(R_mean, tem_mean)
      }
    }
    colnames(R_mean) <- 1:Num_cell_type

    LRmean_list <- list()
    LRmean_list[[1]] <- L_mean
    LRmean_list[[2]] <- R_mean
    LRmean_list[[3]] <- UniqueLR_mean
    LRmean_list[[4]] <- GroupedUniqueLRMatrix
    LRmean_list[[5]] <- GroupedNum_Colnames
    names(LRmean_list) <- c("L_mean","R_mean","UniqueLR_mean","GroupedUniqueLRMatrix","GroupedNum_Colnames")
    return(LRmean_list)
  }

  RandomMean <- function(LRmatrix_list,data_list,PercentForPval = 0.95,ncore = 2, MultiThreshold = FALSE){
    L_dd <- LRmatrix_list[[1]]
    R_dd <- LRmatrix_list[[2]]
    dd <- data_list[[1]]
    Num_cell_type <- length(unique(colnames(dd)))
    random_aa_mean <- data.frame(matrix(ncol = 1000, nrow =  nrow(L_dd)))#generate a dataframe to deposite the mean values for random generated new matrix
    random_bb_mean <- data.frame(matrix(ncol = 1000, nrow =  nrow(R_dd)))
    LR_dd <- rbind(L_dd,R_dd)


    RandomMeanMatrix <- function(LR_dd,Num_cell_type){

      random_LR <- LR_dd[,sample(1:ncol(LR_dd))]
      colnames(random_LR) <- colnames(LR_dd)

      LR_Random_Mean <- NULL
      for (i in 1:Num_cell_type){
        if (length(random_LR[1, colnames(random_LR) == i]) != 1) { #considering the effect of only one cell type.remove the error in mean function when there is only one cell in this cell type.
          j <- rowMeans(random_LR[, colnames(random_LR) == i])
          LR_Random_Mean <- cbind(LR_Random_Mean, j)
        } else {
          j <- random_LR[, colnames(random_LR) == i]
          LR_Random_Mean <- cbind(LR_Random_Mean, j)
        }
      }
      colnames(LR_Random_Mean) <- 1:Num_cell_type

      RandomMeanList <- list()
      RandomMeanList[[1]] <- LR_Random_Mean[1:(nrow(LR_Random_Mean)/2),]
      RandomMeanList[[2]] <- LR_Random_Mean[((nrow(LR_Random_Mean)/2)+1):nrow(LR_Random_Mean),]
      return(RandomMeanList)
    }
    # if (nrow(dd) == nrow(SampleData)){
    #   ncore = 3
    # } else if (nrow(dd) != nrow(SampleData)){
    #   ncore <- parallel::detectCores()
    # }
    # ncore <- parallel::detectCores()
    # set.seed(123)
    # ncore = 9
    set.seed(12345)
    seeds.random <- sample(1:1000)
    cl <- parallel::makeCluster(ncore)
    doParallel::registerDoParallel(cl)
    randomMeanList <- as.list(rep(0,1000))
    ptm <- proc.time()
    randomMeanList <- foreach::foreach(k = 1:1000) %dopar%
    {
      set.seed(seeds.random[k])
      RandomMeanMatrix(LR_dd,Num_cell_type)
    }
    parallel::stopCluster(cl)
    proc.time() - ptm
    random_aa_mean <- lapply(randomMeanList,function(x){x[[1]]})
    random_bb_mean <- lapply(randomMeanList,function(x){x[[2]]})
    aa_matrix <- do.call(cbind,random_aa_mean)
    bb_matrix <- do.call(cbind,random_bb_mean)
    p <- apply(aa_matrix,1,function(x){sort(x)[ncol(aa_matrix)*PercentForPval]})
    k <- apply(bb_matrix,1,function(x){sort(x)[ncol(aa_matrix)*PercentForPval]})
    p[which(p == 0)] <- 0.01
    k[which(k==0)] <- 0.01

    if (MultiThreshold == TRUE){
      p[which(p < FixedThreshold)] <-  FixedThreshold
      k[which(k < FixedThreshold)] <- FixedThreshold
    }

    RandomMean_list <- list()
    RandomMean_list[[1]] <- random_aa_mean
    RandomMean_list[[2]] <- random_bb_mean
    RandomMean_list[[3]] <- p
    RandomMean_list[[4]] <- k
    names(RandomMean_list) <- c("random_aa_mean","random_bb_mean","p","k")
    return(RandomMean_list)
  }

  LRC3_Original <- function(LRmean_list,RandomMean_list,PvalAsThreshold = TRUE, FixedThreshold = 1){

    L_mean <- LRmean_list[[1]]
    R_mean <- LRmean_list[[2]]
    p <- RandomMean_list[[3]]
    k <- RandomMean_list[[4]]
    contactINF <- array(0,dim = c(ncol(L_mean),ncol(L_mean),nrow(L_mean)))
    scaled_contactINF <- array(0,dim = c(ncol(L_mean),ncol(L_mean),nrow(L_mean)))
    LigandExpressionLevelINF <- array(0,dim = c(ncol(L_mean),ncol(L_mean),nrow(L_mean)))
    ReceptorExpressionLevelINF <- array(0,dim = c(ncol(L_mean),ncol(L_mean),nrow(L_mean)))
    ScaledLigandExpressionLevelINF <- array(0,dim = c(ncol(L_mean),ncol(L_mean),nrow(L_mean)))
    ScaledReceptorExpressionLevelINF <- array(0,dim = c(ncol(L_mean),ncol(L_mean),nrow(L_mean)))
    #
    for (i in 1:ncol(L_mean)){
      for (j in 1:ncol(L_mean)){
        c1 <- L_mean[,i]
        c2 <- R_mean[,j]
        LigandExpressionLevelINF[i,j,] <- c1
        ReceptorExpressionLevelINF[i,j,] <- c2
        # ScaledLigandExpressionLevelINF[i,j,] <- c1/p
        # ScaledReceptorExpressionLevelINF[i,j,] <- c2/k
        if (PvalAsThreshold == TRUE){
          ScaledLigandExpressionLevelINF[i,j,] <- c1/p
          ScaledReceptorExpressionLevelINF[i,j,] <- c2/k

          vec <- as.numeric(c1 > p & c2 > k)
          contactINF[i,j,] <- vec
          # ScaledLigandExpressionLevelINF[i,j,] <- ScaledLigandExpressionLevelINF[i,j,]*vec
          # ScaledReceptorExpressionLevelINF[i,j,] <- ScaledReceptorExpressionLevelINF[i,j,]*vec
          vec[which(vec > 0)] <- (c1[which(vec > 0)]/p[which(vec > 0)]+c2[which(vec > 0)]/k[which(vec > 0)])/2

        } else if (PvalAsThreshold == FALSE){
          ScaledLigandExpressionLevelINF[i,j,] <- c1/FixedThreshold
          ScaledReceptorExpressionLevelINF[i,j,] <- c2/FixedThreshold

          vec <- as.numeric(c1 > FixedThreshold & c2 > FixedThreshold)
          contactINF[i,j,] <- vec
          # ScaledLigandExpressionLevelINF[i,j,] <- ScaledLigandExpressionLevelINF[i,j,]*vec
          # ScaledReceptorExpressionLevelINF[i,j,] <- ScaledReceptorExpressionLevelINF[i,j,]*vec
          vec[which(vec > 0)] <- (c1[which(vec > 0)]/FixedThreshold + c2[which(vec > 0)]/FixedThreshold)/2
        }
        if (length(which(vec > 10))){
          vec[which(vec>10)] <- 10
        }
        scaled_contactINF[i,j,] <- vec
      }
    }

    contact_number_original <- matrix(0,ncol = ncol(L_mean), nrow = ncol(L_mean))
    for (i in 1:ncol(L_mean)){
      contact_number_original[i,] <- apply(contactINF[i,,],1,sum)
    }
    scaled_contact_number_original <- matrix(0,ncol = ncol(L_mean), nrow = ncol(L_mean))
    for (i in 1:ncol(L_mean)){
      scaled_contact_number_original[i,] <- apply(scaled_contactINF[i,,],1,sum)

    }

    LRC3_Original_list <- list()
    LRC3_Original_list[[1]] <- contact_number_original
    LRC3_Original_list[[2]] <- scaled_contact_number_original
    LRC3_Original_list[[3]] <- contactINF
    LRC3_Original_list[[4]] <- scaled_contactINF
    LRC3_Original_list[[5]] <- LigandExpressionLevelINF
    LRC3_Original_list[[6]] <- ReceptorExpressionLevelINF
    LRC3_Original_list[[7]] <- ScaledLigandExpressionLevelINF
    LRC3_Original_list[[8]] <- ScaledReceptorExpressionLevelINF
    names(LRC3_Original_list) <- c("contact_number_original","scaled_contact_number_original","contactINF","scaled_contactINF","LigandExpressionLevelINF","ReceptorExpressionLevelINF","ScaledLigandExpressionLevelINF","ScaledReceptorExpressionLevelINF")
    return(LRC3_Original_list)
  }


  LRPairsOriginal_CellTypeSearch <- function(CellType1,CellType2,LRC3_Original_list,LRmatrix_list){
    contactINF <- LRC3_Original_list[[3]]
    rem_LRpairs <- LRmatrix_list[[3]]
    LRPairs_Vec <- contactINF[CellType1,CellType2,]
    LRPairs_InUse <- rem_LRpairs[which(LRPairs_Vec == 1),]
    rownames(LRPairs_InUse) <- c(1:nrow(LRPairs_InUse))
    colnames(LRPairs_InUse) <- c("Ligand","Receptor")
    LRPairs_InUse_Transpose <- t(LRPairs_InUse)
    Table_LRpairs_InUse <- as.table(LRPairs_InUse_Transpose)
    return(Table_LRpairs_InUse)
  }

  LRC3_Random <- function(data_list,RandomMean_list,LRmean_list,PvalAsThreshold = TRUE, FixedThreshold = 1,PercentForPtest = 0.95,ncore = 2){

    random_aa_mean <- RandomMean_list[[1]]
    random_bb_mean <- RandomMean_list[[2]]
    L_mean <- LRmean_list[[1]]
    R_mean <- LRmean_list[[2]]
    dd <- data_list[[1]]

    TotalSize <- nrow(random_bb_mean[[1]])*ncol(random_aa_mean[[1]])*ncol(random_aa_mean[[1]])
      if (TotalSize <= 1.25*10^6) {
      # if (nrow(dd) == nrow(SampleData)){
      #   ncore = 3
      # } else if (nrow(dd) != nrow(SampleData)){
      #   ncore <- parallel::detectCores()
      # }
      # ncore = 9
      # ncore <- detectCores()
      cl <- makeCluster(ncore)
      registerDoParallel(cl)
      LRC3_Random_list <- as.list(rep(0,1000))
      LRC3_Random_list <- foreach::foreach(i = 1:1000) %dopar%
      {
        LRC3_Contacts <- function(Random_aa_Mean,Random_bb_Mean){

          L_mean <- as.matrix(Random_aa_Mean)
          R_mean <- as.matrix(Random_bb_Mean)

          p <- RandomMean_list[[3]]
          k <- RandomMean_list[[4]]


          contactINF <- array(0,dim = c(ncol(L_mean),ncol(L_mean),nrow(L_mean)))
          scaled_contactINF <- array(0,dim = c(ncol(L_mean),ncol(L_mean),nrow(L_mean)))

          for (i in 1:ncol(L_mean)){
            for (j in 1:ncol(L_mean)){
              c1 <- L_mean[,i]
              c2 <- R_mean[,j]
              if (PvalAsThreshold == TRUE){
                vec <- as.numeric(c1 > p & c2 > k)
                contactINF[i,j,] <- vec
                vec[which(vec > 0)] <- (c1[which(vec > 0)]/p[which(vec > 0)]+c2[which(vec > 0)]/k[which(vec > 0)])/2

              } else if (PvalAsThreshold == FALSE){
                vec <- as.numeric(c1 > FixedThreshold & c2 > FixedThreshold)
                contactINF[i,j,] <- vec
                vec[which(vec > 0)] <- (c1[which(vec > 0)]/FixedThreshold + c2[which(vec > 0)]/FixedThreshold)/2

              }

              if (length(which(vec > 10))){
                vec[which(vec>10)] <- 10
              }
              scaled_contactINF[i,j,] <- vec
            }
          }
          LRC3_Contact_list <- list()
          LRC3_Contact_list[[1]] <- contactINF
          LRC3_Contact_list[[2]] <- scaled_contactINF
          names(LRC3_Contact_list) <- c("contactINF","scaled_contactINF")
          return(LRC3_Contact_list)
        }
        LRC3_Contacts(Random_aa_Mean = random_aa_mean[[i]],Random_bb_Mean = random_bb_mean[[i]])
      }
      stopCluster(cl)

      contactINF_list <- lapply(LRC3_Random_list,function(x){x[[1]]})
      scaled_contactINF_list <- lapply(LRC3_Random_list,function(x){x[[2]]})

      Ncolumns <- ncol(random_aa_mean[[1]])
      P_matrix <- data.frame(matrix(ncol = Ncolumns,nrow = Ncolumns))
      highly_occurence_positions <- array(0,dim = c(Ncolumns,Ncolumns,2000))
      Frequency_highly_occurence_positions <- array(0,dim = c(Ncolumns,Ncolumns,2000))

      for (i in 1: Ncolumns){
        for (j in 1: Ncolumns){
          VEC <- t(sapply(contactINF_list,function(x){x[i,j,]}))
          OccurenceEachLRpairs <- apply(VEC,2,function(x){sum(x)})
          PositionTenMore <- which(OccurenceEachLRpairs > (1000*(1-PercentForPtest)))
          if (length(PositionTenMore) > 0){
            OccurenceTenMore <- OccurenceEachLRpairs[PositionTenMore]
            OrderedPositionTenMore <- PositionTenMore[order(OccurenceTenMore,decreasing = T)]
            OrderedOccurenceTenMore <- OccurenceTenMore[order(OccurenceTenMore,decreasing = T)]
            highly_occurence_positions[i,j,1:length(PositionTenMore)] <- OrderedPositionTenMore
            Frequency_highly_occurence_positions[i,j,1:length(PositionTenMore)] <- OrderedOccurenceTenMore
          }
          P_vec <- sapply(contactINF_list,function(x){sum(x[i,j,])})
          P_matrix[i,j] <- sort(P_vec)[length(P_vec)*0.95]#PtestForContactsNumber
        }
      }



    } else if (TotalSize > 1.25*10^6){
      if (ncol(random_aa_mean[[1]])/10 <= 2){
        Num_Steps <- 1
      } else if (ncol(random_aa_mean[[1]])/10 > 2 & ncol(random_aa_mean[[1]])/10 <= 3){
        Num_Steps <- 2
      } else if (ncol(random_aa_mean[[1]])/10 > 3 & ncol(random_aa_mean[[1]])/10 <= 4){
        Num_Steps <- 4
      } else if (ncol(random_aa_mean[[1]])/10 > 4 & ncol(random_aa_mean[[1]])/10 <= 5){
        Num_Steps <- 5
      } else if (ncol(random_aa_mean[[1]])/10 > 5 & ncol(random_aa_mean[[1]])/10 <= 6){
        Num_Steps <- 8
      }  else if (ncol(random_aa_mean[[1]])/10 > 6 & ncol(random_aa_mean[[1]])/10 <=10){
        Num_Steps <- 25
      } else if (ncol(random_aa_mean[[1]])/10 > 10) {
        print("Warning! The cell type number exceed maximum!")
      }
      if (ncol(random_aa_mean[[1]])/10 > 10){
        # break
        stop("Warning! The cell type number exceed maximum!")
      }
      # Num_Steps <- 10
      List_P_matrix <- list()
      List_highly_occurence_positions <- list()
      List_Frequency_highly_occurence_positions <- list()
      for (ns in 1:Num_Steps) {

        # if (nrow(dd) == nrow(SampleData)){
        #   ncore = 3
        # } else if (nrow(dd) != nrow(SampleData)){
        #   ncore <- parallel::detectCores()
        # }
        # ncore = 9
        # ncore <- detectCores()
        cl <- makeCluster(ncore)
        registerDoParallel(cl)
        LRC3_Random_list <- as.list(rep(0,1000))
        LRC3_Random_list <- foreach::foreach(l = 1:1000) %dopar%
        {
          LRC3_Contacts <- function(Random_aa_Mean,Random_bb_Mean,index_ns){

            L_mean <- as.matrix(Random_aa_Mean)
            R_mean <- as.matrix(Random_bb_Mean)

            p <- RandomMean_list[[3]][index_ns]
            k <- RandomMean_list[[4]][index_ns]


            contactINF <- array(0,dim = c(ncol(L_mean),ncol(L_mean),nrow(L_mean)))
            scaled_contactINF <- array(0,dim = c(ncol(L_mean),ncol(L_mean),nrow(L_mean)))

            for (i in 1:ncol(L_mean)){
              for (j in 1:ncol(L_mean)){

                c1 <- L_mean[,i]
                c2 <- R_mean[,j]
                if (PvalAsThreshold == TRUE){
                  vec <- as.numeric(c1 > p & c2 > k)
                  contactINF[i,j,] <- vec
                  vec[which(vec > 0)] <- (c1[which(vec > 0)]/p[which(vec > 0)]+c2[which(vec > 0)]/k[which(vec > 0)])/2

                } else if (PvalAsThreshold == FALSE){
                  vec <- as.numeric(c1 > FixedThreshold & c2 > FixedThreshold)
                  contactINF[i,j,] <- vec
                  vec[which(vec > 0)] <- (c1[which(vec > 0)]/FixedThreshold + c2[which(vec > 0)]/FixedThreshold)/2

                }

                if (length(which(vec > 10))){
                  vec[which(vec>10)] <- 10
                }
                scaled_contactINF[i,j,] <- vec
              }
            }
            LRC3_Contact_list <- list()
            LRC3_Contact_list[[1]] <- contactINF
            LRC3_Contact_list[[2]] <- scaled_contactINF
            names(LRC3_Contact_list) <- c("contactINF","scaled_contactINF")
            return(LRC3_Contact_list)
          }

          if (ns != Num_Steps){
            index_ns <- (((ns-1)*(round(nrow(random_aa_mean[[l]])/Num_Steps))) +1 ):(ns*(round(nrow(random_aa_mean[[l]])/Num_Steps)))
            Random_aa_Mean <- random_aa_mean[[l]][(((ns-1)*(round(nrow(random_aa_mean[[l]])/Num_Steps))) +1 ):(ns*(round(nrow(random_aa_mean[[l]])/Num_Steps))),]
            Random_bb_Mean <- random_bb_mean[[l]][(((ns-1)*(round(nrow(random_bb_mean[[l]])/Num_Steps))) +1 ):(ns*(round(nrow(random_bb_mean[[l]])/Num_Steps))),]
          } else if (ns == Num_Steps){
            index_ns <- (((ns-1)*(round(nrow(random_bb_mean[[l]])/Num_Steps))) +1 ):nrow(random_bb_mean[[l]])
            Random_aa_Mean <- random_aa_mean[[l]][(((ns-1)*(round(nrow(random_aa_mean[[l]])/Num_Steps))) +1 ):nrow(random_aa_mean[[l]]),]
            Random_bb_Mean <- random_bb_mean[[l]][(((ns-1)*(round(nrow(random_bb_mean[[l]])/Num_Steps))) +1 ):nrow(random_bb_mean[[l]]),]
          }

          LRC3_Contacts(Random_aa_Mean = Random_aa_Mean,Random_bb_Mean = Random_bb_Mean,index_ns =index_ns)
          # LRC3_Contacts(Random_aa_Mean = random_aa_mean[[i]][1:(round(nrow(random_aa_mean[[i]])/Num_Steps)),],Random_bb_Mean = random_bb_mean[[i]][1:(round(nrow(random_bb_mean[[i]])/Num_Steps)),])
        }
        stopCluster(cl)

        contactINF_list <- lapply(LRC3_Random_list,function(x){x[[1]]})
        scaled_contactINF_list <- lapply(LRC3_Random_list,function(x){x[[2]]})

        Ncolumns <- ncol(random_aa_mean[[1]])
        # P_matrix <- data.frame(matrix(ncol = Ncolumns,nrow = Ncolumns))
        P_matrix <- array(0,dim = c(Ncolumns,Ncolumns,ncol(contactINF_list[[1]][1,,])))
        highly_occurence_positions <- array(0,dim = c(Ncolumns,Ncolumns,ncol(contactINF_list[[1]][1,,])))
        Frequency_highly_occurence_positions <- array(0,dim = c(Ncolumns,Ncolumns,ncol(contactINF_list[[1]][1,,])))

        for (i in 1: Ncolumns){
          for (j in 1: Ncolumns){
            VEC <- t(sapply(contactINF_list,function(x){x[i,j,]}))
            OccurenceEachLRpairs <- apply(VEC,2,function(x){sum(x)})
            PositionTenMore <- which(OccurenceEachLRpairs > (1000*(1-PercentForPtest)))
            if (length(PositionTenMore) > 0){
              OccurenceTenMore <- OccurenceEachLRpairs[PositionTenMore]
              OrderedPositionTenMore <- PositionTenMore[order(OccurenceTenMore,decreasing = T)]
              OrderedOccurenceTenMore <- OccurenceTenMore[order(OccurenceTenMore,decreasing = T)]
              highly_occurence_positions[i,j,1:length(PositionTenMore)] <- OrderedPositionTenMore
              Frequency_highly_occurence_positions[i,j,1:length(PositionTenMore)] <- OrderedOccurenceTenMore
            }
            P_vec <- sapply(contactINF_list,function(x){sum(x[i,j,])})
            # P_matrix[i,j] <- sort(P_vec)[length(P_vec)*PtestForContactsNumber]
            P_matrix[i,j,] <- P_vec
          }
        }
        List_P_matrix[[ns]] <- P_matrix
        List_highly_occurence_positions[[ns]] <- highly_occurence_positions
        List_Frequency_highly_occurence_positions[[ns]] <- Frequency_highly_occurence_positions
      }
      P_matrix <- do.call(abind::abind,List_P_matrix)
      highly_occurence_positions <- do.call(abind::abind,List_highly_occurence_positions)
      Frequency_highly_occurence_positions <- do.call(abind::abind,List_Frequency_highly_occurence_positions)
      OrderedOccurence <- array(dim = dim(highly_occurence_positions))
      P_matrix2 <- matrix(ncol = Ncolumns,nrow = Ncolumns)
      for (i in 1:Ncolumns){
        for (j in 1:Ncolumns){
          OrderedOccurence[i,j,] <- order(Frequency_highly_occurence_positions[i,j,],decreasing = T)
          highly_occurence_positions[i,j,] <- highly_occurence_positions[i,j,][OrderedOccurence[i,j,]]
          Frequency_highly_occurence_positions[i,j,] <- Frequency_highly_occurence_positions[i,j,][OrderedOccurence[i,j,]]
          P_matrix2[i,j] <- sort(P_matrix[i,j,])[length(P_matrix[i,j,])*0.95]#PtestForContactsNumber
        }
        # OrderedOccurence[i,,] <- t(apply(Frequency_highly_occurence_positions[i,,],1,function(x){order(x,decreasing = T)}))
      }
      P_matrix <- P_matrix2
    }

    ########################################################################################################################################################################################################
    LRC3_Random_RMList <- list()
    LRC3_Random_RMList[[1]] <- P_matrix
    LRC3_Random_RMList[[2]] <- highly_occurence_positions
    LRC3_Random_RMList[[3]] <- Frequency_highly_occurence_positions
    return(LRC3_Random_RMList)
  }


  LRC3_Remaining <- function(LRC3_Original_list,LRmatrix_list,LRC3_Random_RMList,LRmean_list,data_list,RandomMean_list,PvalAsThreshold = TRUE, FixedThreshold = 1){
    contact_number_original <- LRC3_Original_list[[1]]
    scaled_contact_number_original <- LRC3_Original_list[[2]]
    contactINF <- LRC3_Original_list[[3]]
    scaled_contactINF <- LRC3_Original_list[[4]]
    LigandExpressionLevelINF <- LRC3_Original_list[[5]]
    ReceptorExpressionLevelINF <- LRC3_Original_list[[6]]
    ScaledLigandExpressionLevelINF <- LRC3_Original_list$ScaledLigandExpressionLevelINF
    ScaledReceptorExpressionLevelINF <- LRC3_Original_list$ScaledReceptorExpressionLevelINF
    rem_LRpairs <- LRmatrix_list[[3]]
    P_matrix <- LRC3_Random_RMList[[1]]
    highly_occurence_positions <- LRC3_Random_RMList[[2]]
    Frequency_highly_occurence_positions <- LRC3_Random_RMList[[3]]
    L_mean <- LRmean_list[[1]]
    R_mean <- LRmean_list[[2]]
    CellTypeNames <- data_list[[3]]
    p <- RandomMean_list[[3]]
    k <- RandomMean_list[[4]]
    if (PvalAsThreshold == FALSE) {
      p[1:length(p)] <- FixedThreshold
      k[1:length(k)] <- FixedThreshold
    }


    Remaining_contactINF <- contactINF
    Remaining_scaled_contactINF <- scaled_contactINF


    for (i in 1:ncol(highly_occurence_positions)){
      for (j in 1:ncol(highly_occurence_positions)){

        Original_vec <- contactINF[i,j,]
        LR_Original <- which(Original_vec > 0)
        random_OccurenceLR <- as.numeric(highly_occurence_positions[i,j,])
        random_OccurenceLR <- random_OccurenceLR[which(random_OccurenceLR != 0)]
        Remove_LR <- LR_Original[which(LR_Original %in% random_OccurenceLR)]
        Remaining_contactINF[i,j,][Remove_LR] <- 0
        Remaining_scaled_contactINF[i,j,][Remove_LR] <- 0
      }
    }
    LigandExpressionLevelINF <- LigandExpressionLevelINF*Remaining_contactINF
    ReceptorExpressionLevelINF <- ReceptorExpressionLevelINF*Remaining_contactINF
    ScaledLigandExpressionLevelINF <- ScaledLigandExpressionLevelINF*Remaining_contactINF
    ScaledLigandExpressionLevelINF[which(ScaledLigandExpressionLevelINF>=7)] <- 7
    ScaledReceptorExpressionLevelINF <- ScaledReceptorExpressionLevelINF*Remaining_contactINF
    ScaledReceptorExpressionLevelINF[which(ScaledReceptorExpressionLevelINF>=7)] <- 7

    contact_number_Remaining <- matrix(0,ncol = ncol(P_matrix), nrow = ncol(P_matrix))
    for (i in 1:ncol(P_matrix)){
      contact_number_Remaining[i,] <- apply(Remaining_contactINF[i,,],1,sum)
      biological_colnames <- CellTypeNames
      colnames(contact_number_Remaining) <- biological_colnames
      rownames(contact_number_Remaining) <- biological_colnames
    }
    scaled_contact_number_remaining <- matrix(0,ncol = ncol(P_matrix), nrow = ncol(P_matrix))
    for (i in 1:ncol(P_matrix)){
      scaled_contact_number_remaining[i,] <- apply(Remaining_scaled_contactINF[i,,],1,sum)

      biological_colnames <- CellTypeNames
      colnames(scaled_contact_number_remaining) <- biological_colnames
      rownames(scaled_contact_number_remaining) <- biological_colnames
    }
    PtestContactNumber_Index <- matrix(as.numeric(contact_number_Remaining > P_matrix),ncol = ncol(contact_number_Remaining))
    PtestContactNumber <- contact_number_Remaining*PtestContactNumber_Index

    LRC3_Remaining_list <- list()
    LRC3_Remaining_list[[1]] <- contact_number_Remaining
    LRC3_Remaining_list[[2]] <- scaled_contact_number_remaining
    LRC3_Remaining_list[[3]] <- Remaining_contactINF
    LRC3_Remaining_list[[4]] <- Remaining_scaled_contactINF
    LRC3_Remaining_list[[5]] <- rem_LRpairs
    LRC3_Remaining_list[[6]] <- LigandExpressionLevelINF
    LRC3_Remaining_list[[7]] <- ReceptorExpressionLevelINF
    LRC3_Remaining_list[[8]] <- ScaledLigandExpressionLevelINF
    LRC3_Remaining_list[[9]] <- ScaledReceptorExpressionLevelINF
    LRC3_Remaining_list[[10]] <- p
    LRC3_Remaining_list[[11]] <- k
    LRC3_Remaining_list[[12]] <- LRmean_list[[4]]
    LRC3_Remaining_list[[13]] <- LRmean_list[[5]]
    LRC3_Remaining_list[[14]] <- PtestContactNumber
    names(LRC3_Remaining_list) <- c("contact_number_Remaining","scaled_contact_number_remaining","Remaining_contactINF","Remaining_scaled_contactINF","rem_LRpairs","LigandExpressionLevelINF","ReceptorExpressionLevelINF","ScaledLigandExpressionLevelINF","ScaledReceptorExpressionLevelINF","p","k","GroupedUniqueLRMatrix","GroupedNum_Colnames","PtestContactNumber")
    return(LRC3_Remaining_list)
  }

  LRmatrix_list <- LRmatrix(data_list,DefaultLRlist,Species)#get the L_dd and R_dd

  if (sum(LRmatrix_list[[1]]) == 0){
    print("Error In the data, no Ligands or Receptors found in the data!")
  }else if (sum(LRmatrix_list[[1]]) != 0){
    LRmean_list <- LR_mean(LRmatrix_list,data_list)
    RandomMean_list <- RandomMean(LRmatrix_list,data_list,PercentForPval,ncore,MultiThreshold)
    LRC3_Original_list <- LRC3_Original(LRmean_list,RandomMean_list,PvalAsThreshold, FixedThreshold)
    LRC3_Random_RMList <- LRC3_Random(data_list,RandomMean_list,LRmean_list,PvalAsThreshold, FixedThreshold,PercentForPtest,ncore)
    LRC3_list <- LRC3_Remaining(LRC3_Original_list,LRmatrix_list,LRC3_Random_RMList,LRmean_list,data_list,RandomMean_list,PvalAsThreshold, FixedThreshold)
    return(LRC3_list)
  }

}

#'
#'The Contact Matrix Plot
#'
#' This function performs the Plot on the contacts between different cell types.
#' It gives 4 different plots in total.
#' IndexNumber 1 is the matrix showing the contact number between different cell types in two direction: ligand to receptor(OUT), receptor to ligand(IN).
#' IndexNumber 2 is the scaled matrix showing the scaled contact number between different cell types. The scaled number is the ratio of the expression level of that gene to the threshold.
#' IndexNumber 3 is the undirectional contact matrix,which add the In and OUT signals between two cell types.
#' IndexNumber 4 is the undirectional scaled contact matrix.
#' Both one column and one row in the plot represent one cell type,respectively.
#' The color in the plot represents the contact number.
#'
#' @param LRC3_list  the input data generated by function LRC3_INF.
#' @param IndexNumber the Index of the plot.
#' @param Cluster if \code{Cluster =FALSE}, the visualization of the result will not be clustered.
#' @param SaveImage if \code{SaveImage =FALSE}, save image as pdf file.
#' @param FileName if \code{SaveImage =TRUE}, give the name to the pdf file.
#' @importFrom pheatmap pheatmap
#' @export
LRC3_ContactMatrixPlot <- function(LRC3_list,IndexNumber,Cluster =FALSE,SaveImage = FALSE,FileName = "Contactmatrix"){
  LRC3_PlotMatrix <- list()
  LRC3_PlotMatrix[[1]] <- LRC3_list[[1]]
  LRC3_PlotMatrix[[2]] <- LRC3_list[[2]]
  LRC3_PlotMatrix[[3]] <- LRC3_list[[1]] +t(LRC3_list[[1]])
  diag(LRC3_PlotMatrix[[3]]) <- diag(LRC3_PlotMatrix[[3]])/2
  LRC3_PlotMatrix[[4]] <- LRC3_list[[2]]+t(LRC3_list[[2]])
  diag(LRC3_PlotMatrix[[4]]) <- diag(LRC3_PlotMatrix[[4]])/2
  # LRC3_PlotMatrix[[5]] <- LRC3_list[[14]]

  if (sum(LRC3_PlotMatrix[[IndexNumber]]) > 0){
    if (SaveImage == FALSE){
      pheatmap::pheatmap(LRC3_PlotMatrix[[IndexNumber]], cluster_rows = Cluster, cluster_cols = Cluster)

    } else if (SaveImage == TRUE){
      pdf(file= paste0(FileName,".pdf"))
      pheatmap::pheatmap(LRC3_PlotMatrix[[IndexNumber]], cluster_rows = Cluster, cluster_cols = Cluster)
      dev.off()
    }
  } else if (sum(LRC3_PlotMatrix[[IndexNumber]]) == 0){
    print("There is no contact in this dataset!")
  }
}

#' Visualize the connections between different cell types.
#'
#' This function extracts and plots connections between cell types using different ligand-receptor pathways
#' based on the usage frequencies and expression level of the ligand-receptor pairs.
#' Ligand-receptor pairs that have been used for more than 1 times or only 1 time among cell types are shown seperately.
#' IndexNumber 1 shows connections of ligand-receptor pairs used more than 1 times. By choosing different parameters,
#' users can choose to view connections based on the mean expression level of a cell type,
#' or proportion of cells in a cell type that significantly express the specific pathways.
#' IndexNumber 2 shows the usage frequencies of each specific ligand or receptor by each cell type.
#' Due to significant test, some connections of certain pathways between cell types are removed during calculation.
#' IndexNumber 3 shows the connection of ligand-receptor pairs used only once among cell types.
#' By setting the parameter, users can view connections based on mean expression of a cell type,
#' or proportion of cells in a cell type that significantly express the specific ligands and receptors.
#'
#' @param LRC3_list the input data generated by function LRC3_INF.
#' @param IndexOfPlot the Index of the plot.
#' if \code{IndexOfPlot = 1},view connections with ligand-receptor pairs used more than once among cell types.
#' if \code{IndexOfPlot = 2},view the usage frequencies of each ligand or receptor by each cell type.
#' if \code{IndexOfPlot = 3},view the connections with ligand-receptor pairs used only once among cell types.
#' @param By_Proportion if \code{By_Proportion = FALSE}, view the connections based on the mean expression of cell types.
#' if \code{By_Proportion = TRUE}, view the connections based on the proportion of cells within a cell type that significantly express the specific ligand or receptor.
#' @param MultiFigure parameter for plotting multiple figures or ONE figure.
#' if \code{MultiFigure = FALSE}, show all the connections on ONE figure.
#' The figure would be very intense when there are a large number of connections.
#' if \code{MultiFigure = FALSE}, show connections on multiple figures. It gives a maximun number of 5 figures based on the total number of connections.
#' Apart from the last one figure, each figure shows a fixed number of 50 ligand-receptor patheays. The last figure shows all the remaining pairs.
#' The ligand-receptor pairs is first ranked by their usage frequency by a decreasing order.
#' Then pairs with same usage frequency are further ranked by their expression level by a decreasing order.
#' @param PlotResults if \code{PlotResults = TRUE}, plot out the figures.
#' @param SaveImage if \code{SaveImage = FALSE},plot without save the figure. if \code{SaveImage = TRUE},save the figure as pdf.
#' @param ReturnList if \code{ReturnList = TRUE},return the results as a list.
#' @importFrom pheatmap pheatmap
#' @importFrom RColorBrewer brewer.pal
#' @return LRC3_Connection_list a list that contains three elements.
#' 1st is the table for SharedLRpairs
#' 2nd is the table for UniqueLRpairs
#' 3rd is the Total LRpairs Ranked by mean, including both shared and unique LRpairs
#' @export

LRC3_Connection <- function(LRC3_list,IndexOfPlot = 1,By_Proportion = FALSE,MultiFigure = FALSE,
                           PlotResults =TRUE,ReturnList=FALSE,SaveImage = FALSE){
  LRpairsTotalUsage <- NULL
  ScaledLRpairsTotalUsage <- NULL
  LRpairsTotalUsage <- apply(LRC3_list$Remaining_contactINF,3,sum)
  ScaledLRpairsTotalUsage <- apply(LRC3_list$Remaining_scaled_contactINF,3,sum)
  GroupedUniqueLRMatrix <- LRC3_list$GroupedUniqueLRMatrix
  GroupedNum_Colnames <-LRC3_list$GroupedNum_Colnames
  p <- LRC3_list$p
  k <- LRC3_list$k

  # table(LRpairsTotalUsage)
  TopUsageLRpairs <- LRC3_list$rem_LRpairs[which(LRpairsTotalUsage >= 2),]
  if(nrow(TopUsageLRpairs) == 0 ){
    OrderedTopLRpairsCellTypeConnectionMatrix <- 0
    OrderedTopLRpairsCellTypeConnectionMatrix_frequency <- 0
  } else if (nrow(TopUsageLRpairs) == 1){
    OrderedTopLRpairsCellTypeConnectionMatrix <- 1
    OrderedTopLRpairsCellTypeConnectionMatrix_frequency <- 1
  } else if (nrow(TopUsageLRpairs) > 1){
    RownamesForTopUsageLRpairs <- paste0(TopUsageLRpairs$Ligand.ApprovedSymbol,"--",TopUsageLRpairs$Receptor.ApprovedSymbol)
    TopUsageScaledLigandExpressionLevelINF <- LRC3_list$ScaledLigandExpressionLevelINF[,,which(LRpairsTotalUsage >= 2)]
    TopUsageScaledReceptorExpressionLevelINF <- LRC3_list$ScaledReceptorExpressionLevelINF[,,which(LRpairsTotalUsage >= 2)]

    TopLigandMatrix_Signal <- matrix(0,ncol = nrow(TopUsageScaledLigandExpressionLevelINF[1,,]), nrow = ncol(TopUsageScaledLigandExpressionLevelINF[1,,]))
    for (i in 1:nrow(TopUsageScaledLigandExpressionLevelINF[1,,])){
      TopLigandMatrix_Signal[,i] <- apply(TopUsageScaledLigandExpressionLevelINF[i,,],2,function(x){x[which(x!=0)][1]})
      TopLigandMatrix_Signal[,i][is.na(TopLigandMatrix_Signal[,i])] <- 0
    }
    colnames(TopLigandMatrix_Signal) <- paste0(colnames(LRC3_list$contact_number_Remaining),"-L")

    TopReceptorMatrix_Signal <- matrix(0,ncol = nrow(TopUsageScaledLigandExpressionLevelINF[1,,]), nrow = ncol(TopUsageScaledLigandExpressionLevelINF[1,,]))
    for (i in 1:nrow(TopUsageScaledLigandExpressionLevelINF[1,,])){
      TopReceptorMatrix_Signal[,i] <- apply(TopUsageScaledReceptorExpressionLevelINF[,i,],2,function(x){x[which(x!=0)][1]})
      TopReceptorMatrix_Signal[,i][is.na(TopReceptorMatrix_Signal[,i])] <- 0
    }
    colnames(TopReceptorMatrix_Signal) <- paste0(colnames(LRC3_list$contact_number_Remaining),"-R")

    TopLigandMatrix_Frequency <- matrix(0,ncol = nrow(TopUsageScaledLigandExpressionLevelINF[1,,]), nrow = ncol(TopUsageScaledLigandExpressionLevelINF[1,,]))
    for (i in 1:nrow(TopUsageScaledLigandExpressionLevelINF[1,,])){
      TopLigandMatrix_Frequency[,i] <- apply(TopUsageScaledLigandExpressionLevelINF[i,,],2,function(x){length(which(x > 0))})
    }
    colnames(TopLigandMatrix_Frequency) <- paste0(colnames(LRC3_list$contact_number_Remaining),"-L")
    TopReceptorMatrix_Frequency <- matrix(0,ncol = nrow(TopUsageScaledLigandExpressionLevelINF[1,,]), nrow = ncol(TopUsageScaledLigandExpressionLevelINF[1,,]))
    for (i in 1:nrow(TopUsageScaledLigandExpressionLevelINF[1,,])){
      TopReceptorMatrix_Frequency[,i] <- apply(TopUsageScaledReceptorExpressionLevelINF[,i,],2,function(x){length(which(x > 0))})
    }
    colnames(TopReceptorMatrix_Frequency) <- paste0(colnames(LRC3_list$contact_number_Remaining),"-R")


    TopLRpairsCellTypeConnectionMatrix <- cbind(TopLigandMatrix_Signal,TopReceptorMatrix_Signal)
    TopLRpairsCellTypeConnectionMatrix_frequency <- cbind(TopLigandMatrix_Frequency,TopReceptorMatrix_Frequency)
    rownames(TopLRpairsCellTypeConnectionMatrix) <- paste0(TopUsageLRpairs$Ligand.ApprovedSymbol,"--",TopUsageLRpairs$Receptor.ApprovedSymbol)
    rownames(TopLRpairsCellTypeConnectionMatrix_frequency) <- paste0(TopUsageLRpairs$Ligand.ApprovedSymbol,"--",TopUsageLRpairs$Receptor.ApprovedSymbol)

    OrderedTopLRpairsCellTypeConnectionMatrix <- TopLRpairsCellTypeConnectionMatrix[order(rowSums(TopLRpairsCellTypeConnectionMatrix_frequency),decreasing = T),]
    OrderedTopLRpairsCellTypeConnectionMatrix_frequency <- TopLRpairsCellTypeConnectionMatrix_frequency[order(rowSums(TopLRpairsCellTypeConnectionMatrix_frequency),decreasing = T),]

    FrequencyTopLRs <- sort(rowSums(TopLRpairsCellTypeConnectionMatrix_frequency)/2,decreasing = T)

    Ordered_TopCommonLRpairsConnection <- NULL
    Ordered_TopCommonLRpairsConnection_Frequency <- NULL
    NameForOne <- NULL
    for (i in max(FrequencyTopLRs):min(FrequencyTopLRs)){
      sub_Matrix <- OrderedTopLRpairsCellTypeConnectionMatrix[which(FrequencyTopLRs == i),]
      sub_Matrix_Frequency <- OrderedTopLRpairsCellTypeConnectionMatrix_frequency[which(FrequencyTopLRs == i),]
      if (length(which(FrequencyTopLRs ==i)) != 1){
        vec_NonezeroNum <- apply(sub_Matrix,1,function(x){length(which(x != 0))})
        OrderOfMatrix <- order(rowSums(sub_Matrix)/vec_NonezeroNum,decreasing = T)
        sub_Matrix <- sub_Matrix[OrderOfMatrix,]
        sub_Matrix_Frequency <- sub_Matrix_Frequency[OrderOfMatrix,]
      }
      Ordered_TopCommonLRpairsConnection <- rbind(Ordered_TopCommonLRpairsConnection,sub_Matrix)
      Ordered_TopCommonLRpairsConnection_Frequency <- rbind(Ordered_TopCommonLRpairsConnection_Frequency,sub_Matrix_Frequency)
      if (length(which(FrequencyTopLRs ==i)) == 1){
        NameForOne <- cbind(NameForOne,names(which(FrequencyTopLRs == i)))
      }
    }

    rownames(Ordered_TopCommonLRpairsConnection)[which(rownames(Ordered_TopCommonLRpairsConnection) == "sub_Matrix")] <- c(NameForOne)
    rownames(Ordered_TopCommonLRpairsConnection_Frequency)[which(rownames(Ordered_TopCommonLRpairsConnection_Frequency)== "sub_Matrix_Frequency")] <- c(NameForOne)

    OrderedTopLRpairsCellTypeConnectionMatrix <- Ordered_TopCommonLRpairsConnection
    OrderedTopLRpairsCellTypeConnectionMatrix_frequency <- Ordered_TopCommonLRpairsConnection_Frequency
  }

 # CelltypeConnection_List[[1]] <- OrderedTopLRpairsCellTypeConnectionMatrix

  UniqueLRpairs <- LRC3_list$rem_LRpairs[which(LRpairsTotalUsage == 1),]
  if (nrow(UniqueLRpairs) == 0 ){
    TopUniqueCommunication_RankByMean <- 0

  } else if(nrow(UniqueLRpairs) ==1) {
    TopUniqueCommunication_RankByMean <- 1
  } else if (nrow(UniqueLRpairs) != 0 & nrow(UniqueLRpairs) !=1 ){

    UniqueScaledLigandExpressionLevelINF <- LRC3_list$ScaledLigandExpressionLevelINF[,,which(LRpairsTotalUsage == 1)]
    UniqueScaledReceptorExpressionLevelINF <- LRC3_list$ScaledReceptorExpressionLevelINF[,,which(LRpairsTotalUsage == 1)]

    UniqueLigandMatrix_Signal <- matrix(0,ncol = nrow(UniqueScaledLigandExpressionLevelINF[1,,]), nrow = ncol(UniqueScaledLigandExpressionLevelINF[1,,]))
    # for (i in 1:nrow(TopUsageScaledLigandExpressionLevelINF[1,,])){
    for (i in 1:nrow(LRC3_list[[1]])){
    UniqueLigandMatrix_Signal[,i] <- apply(UniqueScaledLigandExpressionLevelINF[i,,],2,function(x){x[which(x!=0)[1]]})
      UniqueLigandMatrix_Signal[,i][is.na(UniqueLigandMatrix_Signal[,i])] <- 0
    }
    colnames(UniqueLigandMatrix_Signal) <- paste0(colnames(LRC3_list$contact_number_Remaining),"-L")

    UniqueReceptorMatrix_Signal <- matrix(0,ncol = nrow(UniqueScaledReceptorExpressionLevelINF[1,,]), nrow = ncol(UniqueScaledReceptorExpressionLevelINF[1,,]))
    for (i in 1:nrow(LRC3_list[[1]])){
      UniqueReceptorMatrix_Signal[,i] <- apply(UniqueScaledReceptorExpressionLevelINF[,i,],2,function(x){x[which(x!=0)][1]})
      UniqueReceptorMatrix_Signal[,i][is.na(UniqueReceptorMatrix_Signal[,i])] <- 0
    }
    colnames(UniqueReceptorMatrix_Signal) <- paste0(colnames(LRC3_list$contact_number_Remaining),"-R")
    UniqueLRpairsCellTypeConnectionMatrix <- cbind(UniqueLigandMatrix_Signal,UniqueReceptorMatrix_Signal)
    rownames(UniqueLRpairsCellTypeConnectionMatrix) <- paste0(UniqueLRpairs$Ligand.ApprovedSymbol,"--",UniqueLRpairs$Receptor.ApprovedSymbol)

    OrderedUniqueLRpairsCellTypeConnectionMatrix_RankByMin <- UniqueLRpairsCellTypeConnectionMatrix[order(pmin(apply(UniqueLigandMatrix_Signal,1,max),apply(UniqueReceptorMatrix_Signal,1,max)),decreasing = T),]
    OrderedUniqueLRpairsCellTypeConnectionMatrix_RankByMean <- UniqueLRpairsCellTypeConnectionMatrix[order((apply(UniqueLigandMatrix_Signal,1,max) + apply(UniqueReceptorMatrix_Signal,1,max))/2,decreasing = T),]#order by the mean of max(L) and Max(R)


    TopUniqueCommunication_RankByMean <- OrderedUniqueLRpairsCellTypeConnectionMatrix_RankByMean

  }

  if (ReturnList == TRUE){
    CelltypeConnection_List <- list()

    Unique_frequency <- 0
    if (nrow(UniqueLRpairs) != 0 & !(is.null(nrow(OrderedTopLRpairsCellTypeConnectionMatrix)))){
      TotalConnection <- rbind(OrderedTopLRpairsCellTypeConnectionMatrix,TopUniqueCommunication_RankByMean)
      Unique_frequency <- matrix(as.numeric(TopUniqueCommunication_RankByMean != 0),nrow = nrow(TopUniqueCommunication_RankByMean))
      rownames(Unique_frequency) <- rownames(TopUniqueCommunication_RankByMean)
      TotalFrequencyMatrix <- rbind(OrderedTopLRpairsCellTypeConnectionMatrix_frequency,Unique_frequency)
    } else if (nrow(UniqueLRpairs) != 0 & is.null(nrow(OrderedTopLRpairsCellTypeConnectionMatrix))){
      TotalConnection <- TopUniqueCommunication_RankByMean
      Unique_frequency <- matrix(as.numeric(TopUniqueCommunication_RankByMean != 0),nrow = nrow(TopUniqueCommunication_RankByMean))
      rownames(Unique_frequency) <- rownames(TopUniqueCommunication_RankByMean)
      TotalFrequencyMatrix <- Unique_frequency
      # TotalConnection <- OrderedTopLRpairsCellTypeConnectionMatrix
    } else if (nrow(UniqueLRpairs) == 0 & !(is.null(nrow(OrderedTopLRpairsCellTypeConnectionMatrix)))){
      TotalConnection <- OrderedTopLRpairsCellTypeConnectionMatrix
      TotalFrequencyMatrix <- OrderedTopLRpairsCellTypeConnectionMatrix_frequency
    } else if (nrow(UniqueLRpairs) == 0 & is.null(nrow(OrderedTopLRpairsCellTypeConnectionMatrix))){
      TotalConnection <- 0
      print("There is no connection for the input data!")
    }
    # if(sum(Unique_frequency) != 0){
    #   rownames(Unique_frequency) <- rownames(TopUniqueCommunication_RankByMean)
    #   TotalFrequencyMatrix <- rbind(OrderedTopLRpairsCellTypeConnectionMatrix_frequency,Unique_frequency)
    #
    # } else if (sum(Unique_frequency) == 0){
    #   TotalFrequencyMatrix <- OrderedTopLRpairsCellTypeConnectionMatrix_frequency
    # }

    CelltypeConnection_List[[1]] <- TotalConnection
    CelltypeConnection_List[[2]] <- TopUniqueCommunication_RankByMean
    CelltypeConnection_List[[3]] <- TotalFrequencyMatrix

    if (nrow(UniqueLRpairs) != 0 & !(is.null(nrow(OrderedTopLRpairsCellTypeConnectionMatrix)))){
      TotalConnection_RankByMean <- rbind(OrderedTopLRpairsCellTypeConnectionMatrix,TopUniqueCommunication_RankByMean)
    } else if (nrow(UniqueLRpairs) != 0 & is.null(nrow(OrderedTopLRpairsCellTypeConnectionMatrix))){
      TotalConnection_RankByMean <- TopUniqueCommunication_RankByMean
      # TotalConnection_RankByMean <- OrderedTopLRpairsCellTypeConnectionMatrix
    } else if (nrow(UniqueLRpairs) == 0 & !(is.null(nrow(OrderedTopLRpairsCellTypeConnectionMatrix)))){
      TotalConnection_RankByMean <- OrderedTopLRpairsCellTypeConnectionMatrix

    } else if (nrow(UniqueLRpairs) == 0 & is.null(nrow(OrderedTopLRpairsCellTypeConnectionMatrix))){
      TotalConnection_RankByMean <- 0
    }
    TotalConnection_RankByMean_Ligand <- TotalConnection_RankByMean[,1:ncol(TotalConnection_RankByMean)/2]
    TotalConnection_RankByMean_Receptor <- TotalConnection_RankByMean[,((ncol(TotalConnection_RankByMean)/2)+1):ncol(TotalConnection_RankByMean)]
    TotalConnection_RankByMean <- TotalConnection_RankByMean[order((apply(TotalConnection_RankByMean_Ligand,1,max) + apply(TotalConnection_RankByMean_Receptor,1,max)),decreasing = TRUE),]
    CelltypeConnection_List[[4]] <- TotalConnection_RankByMean#this is ranking the LRpairs by mean expression, including shared and unique LRpairs
    names(CelltypeConnection_List) <- c("SharedLRpairs","UniqueLRpairs","FrequencyOfSharedLRpairs","TotalLRpairs_RankByMean")

  }




  ################################the following part is for proportion of cells of Commonly Used LRpairs


    if (By_Proportion == TRUE){
    if (IndexOfPlot == 1){
      TopMatrix <- OrderedTopLRpairsCellTypeConnectionMatrix
    } else if (IndexOfPlot == 3){
      TopMatrix <- TopUniqueCommunication_RankByMean
    }

    if (sum(TopMatrix) == 0){
      Com_CellTypeConnection_dd <- 0
    } else if(sum(TopMatrix) != 0){
      if (length(TopMatrix) ==1) {
        Com_CellTypeConnection_dd <- 1
      }else if(length(TopMatrix) !=1){
        rem_ComLRpairs <- rownames(TopMatrix)
        rem_ComL <- sapply(strsplit(rem_ComLRpairs,"--"),"[[",1)
        rem_ComR <- sapply(strsplit(rem_ComLRpairs,"--"),"[[",2)
        ComLigand_Matrix <- TopMatrix[,1:(ncol(TopMatrix)/2)]
        ComReceptor_Matrix <- TopMatrix[,(ncol(TopMatrix)/2+1):ncol(TopMatrix)]

        ComLigand_Binary <- ComLigand_Matrix
        ComLigand_Binary[which(ComLigand_Binary != 0)] <- 1
        ComReceptor_Binary <- ComReceptor_Matrix
        ComReceptor_Binary[which(ComReceptor_Binary !=0)] <- 1

        index_ComLigand <- sapply(as.list(rem_ComL),function(x){which(rownames(GroupedUniqueLRMatrix) == x)})
        index_ComReceptor <- sapply(as.list(rem_ComR),function(x){which(rownames(GroupedUniqueLRMatrix) == x)})

        ComLigand_dd <- GroupedUniqueLRMatrix[index_ComLigand,]# the single cells for each cell type
        ComReceptor_dd <- GroupedUniqueLRMatrix[index_ComReceptor,]

        ComLigand_colnames <- sapply(strsplit(colnames(ComLigand_Binary),"-"),"[[",1)
        ComReceptor_colnames <- sapply(strsplit(colnames(ComReceptor_Binary),"-"),"[[",1)

        Index_P <- sapply(as.list(rem_ComL),function(x){which(names(p) == x)}[1])#get the positions of remaining p
        Index_K <- sapply(as.list(rem_ComR),function(x){which(names(k) == x)}[1])

        rem_P <- p[Index_P]
        rem_K <- k[Index_K]

        ComLigand_dd_Connection <- NULL
        ComReceptor_dd_Connection <- NULL

        for (i in 1:length(ComLigand_colnames)){
          sub_ComLigand_dd <- ComLigand_dd[,which(GroupedNum_Colnames== i)]
          sub_ComLigand_Binary <- ComLigand_Binary[,which(ComLigand_colnames == ComLigand_colnames[i])]
          tem_ComLigand_dd <- sub_ComLigand_dd*sub_ComLigand_Binary

          tem_ComLigand_dd[which(tem_ComLigand_dd < rem_P)] <- 0
          tem_ComLigand_dd <- tem_ComLigand_dd/rem_P
          tem_ComLigand_dd[which(tem_ComLigand_dd > 7)] <- 7
          colnames(tem_ComLigand_dd) <- paste0(colnames(tem_ComLigand_dd),".",c(1:ncol(tem_ComLigand_dd)))
          ComLigand_dd_Connection <- cbind(ComLigand_dd_Connection,tem_ComLigand_dd)

          sub_ComReceptor_dd <- ComReceptor_dd[,which(colnames(ComReceptor_dd)== i)]
          sub_ComReceptor_Binary <- ComReceptor_Binary[,which(ComReceptor_colnames == ComReceptor_colnames[i])]
          tem_ComReceptor_dd <- sub_ComReceptor_dd*sub_ComReceptor_Binary
          tem_ComReceptor_dd[which(tem_ComReceptor_dd < rem_K)] <- 0
          tem_ComReceptor_dd <- tem_ComReceptor_dd/rem_K
          tem_ComReceptor_dd[which(tem_ComReceptor_dd > 7)] <- 7
          colnames(tem_ComReceptor_dd) <- paste0(colnames(tem_ComReceptor_dd),".",c(1:ncol(tem_ComReceptor_dd)))
          ComReceptor_dd_Connection <- cbind(ComReceptor_dd_Connection,tem_ComReceptor_dd)

        }


        Unique_CellType <- colnames(LRC3_list$contact_number_Remaining)

        rep_Num <- NULL
        rep_vec <- NULL
        for (i in 1:length(Unique_CellType)){
          rep_Num[i] <- length(which(colnames(ComLigand_dd) == i))
          rep_subvec <- rep(i-1,rep_Num[i])
          rep_vec <- c(rep_vec,rep_subvec)

        }
        rep_vec <- rep(rep_vec,2)

        Com_CellTypeConnection_dd <- cbind(ComLigand_dd_Connection,ComReceptor_dd_Connection)
        rownames(Com_CellTypeConnection_dd) <- paste0(rownames(ComLigand_dd_Connection),"--",rownames(ComReceptor_dd_Connection))
        colnames(Com_CellTypeConnection_dd)[1:(ncol(Com_CellTypeConnection_dd)/2)] <- paste0(colnames(Com_CellTypeConnection_dd)[1:(ncol(Com_CellTypeConnection_dd)/2)],"-L")
        colnames(Com_CellTypeConnection_dd)[(ncol(Com_CellTypeConnection_dd)/2+1):ncol(Com_CellTypeConnection_dd)] <- paste0(colnames(Com_CellTypeConnection_dd)[(ncol(Com_CellTypeConnection_dd)/2+1):ncol(Com_CellTypeConnection_dd)],"-R")
      }

    }
      }

  breaksList = seq(min(OrderedTopLRpairsCellTypeConnectionMatrix_frequency)-1,max(OrderedTopLRpairsCellTypeConnectionMatrix_frequency),by = 1)
  breaksList2 = seq(min(OrderedTopLRpairsCellTypeConnectionMatrix)-1,max(OrderedTopLRpairsCellTypeConnectionMatrix)+1,by = 1)

  if (PlotResults == TRUE){
    if (IndexOfPlot == 2) {

      if (nrow(TopUsageLRpairs) == 0){
        print("There is no ligand-receptor pairs used more than 1 time.")
      } else if (nrow(TopUsageLRpairs) == 1){
        print("There is only 1 ligand-receptor pair that used more than 1 time.")
      } else if (nrow(TopUsageLRpairs) > 1){
        # library(pheatmap)
        if (MultiFigure ==FALSE){
          annotation <- data.frame(LR = factor(c(rep(0,ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2),
                                                 rep(1,ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2)),
                                               labels = c("Ligand", "Receptor")))
          rownames(annotation) <- colnames(OrderedTopLRpairsCellTypeConnectionMatrix)
          LR <- c("navyblue", "red")
          names(LR) <- c("Ligand", "Receptor")
          anno_colors <- list(LR = LR)
          unique_labels <- sort(unique(FrequencyTopLRs))

          annotationrow <- data.frame(N = factor(c(FrequencyTopLRs),
                                                 labels = unique_labels))
          rownames(annotationrow) <- rownames(OrderedTopLRpairsCellTypeConnectionMatrix)

          if (SaveImage == TRUE) {
            pdf(file="Rplot%03d.pdf")
          }
          if (SaveImage == FALSE){
            #x11()
          }
          pheatmap::pheatmap(OrderedTopLRpairsCellTypeConnectionMatrix_frequency,
                             cluster_rows = F, cluster_cols = F,
                             fontsize=7,fontsize_col=9, fontsize_row = 7,

                             legend_breaks = breaksList+1,
                             legend_labels = breaksList+1,
                             gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
                             breaks = breaksList,main = "Usage Frequency",
                             color = colorRampPalette(rev(brewer.pal(n = 9, name = "RdYlBu")))(length(breaksList)-1), # Defines the vector of colors for the legend (it has to be of the same lenght of breaksList)
                             annotation_col = annotation,
                             annotation_row = annotationrow,
                             annotation_colors = anno_colors)
        } else if (MultiFigure == TRUE){

          annotation <- data.frame(LR = factor(c(rep(0,ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2),
                                                 rep(1,ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2)),
                                               labels = c("Ligand", "Receptor")))
          rownames(annotation) <- colnames(OrderedTopLRpairsCellTypeConnectionMatrix)
          LR <- c("navyblue", "red")
          names(LR) <- c("Ligand", "Receptor")
          anno_colors <- list(LR = LR)
          unique_labels <- sort(unique(FrequencyTopLRs))

          annotationrow <- data.frame(N = factor(c(FrequencyTopLRs),
                                                 labels = unique_labels))
          rownames(annotationrow) <- rownames(OrderedTopLRpairsCellTypeConnectionMatrix)

          if (nrow(OrderedTopLRpairsCellTypeConnectionMatrix_frequency) <= 51) {
            if (SaveImage == TRUE) {
              pdf(file="Rplot%03d.pdf")
            }
            if (SaveImage == FALSE){
              #x11()
            }
            pheatmap::pheatmap(OrderedTopLRpairsCellTypeConnectionMatrix_frequency,
                               cluster_rows = F, cluster_cols = F,
                               fontsize=7,fontsize_col=7, fontsize_row = 7,
                               legend_breaks = breaksList+1,
                               legend_labels = breaksList+1,
                               gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
                               breaks = breaksList,main = "Usage Frequency",
                               color = colorRampPalette(rev(brewer.pal(n = 9, name = "RdYlBu")))(length(breaksList)-1), # Defines the vector of colors for the legend (it has to be of the same lenght of breaksList)
                               annotation_col = annotation,
                               annotation_row = annotationrow,
                               annotation_colors = anno_colors)#gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
          } else if (nrow(OrderedTopLRpairsCellTypeConnectionMatrix_frequency) > 51 & nrow(OrderedTopLRpairsCellTypeConnectionMatrix_frequency) <= 101){
            if (SaveImage == TRUE) {
              pdf(file="Rplot%03d.pdf")
            }
            if (SaveImage == FALSE){
              #x11()
            }
            pheatmap::pheatmap(OrderedTopLRpairsCellTypeConnectionMatrix_frequency[1:50,],
                               cluster_rows = F, cluster_cols = F,
                               fontsize=7,fontsize_col=7, fontsize_row = 7,
                               legend_breaks = breaksList+1,
                               legend_labels = breaksList+1,
                               breaks = breaksList,main = "Usage Frequency",
                               color = colorRampPalette(rev(brewer.pal(n = 9, name = "RdYlBu")))(length(breaksList)-1), # Defines the vector of colors for the legend (it has to be of the same lenght of breaksList)
                               gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
                               annotation_col = annotation,
                               annotation_row = annotationrow,
                               annotation_colors = anno_colors)#gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
            # dev.copy(jpeg,filename="plot.png");
            # dev.off ();
            if (SaveImage == FALSE){
              #x11()
            }
            pheatmap::pheatmap(OrderedTopLRpairsCellTypeConnectionMatrix_frequency[51:nrow(OrderedTopLRpairsCellTypeConnectionMatrix_frequency),],
                               cluster_rows = F, cluster_cols = F,
                               fontsize=7,fontsize_col=7, fontsize_row = 7,
                               legend_breaks = breaksList+1,
                               legend_labels = breaksList+1,
                               breaks = breaksList,
                               color = colorRampPalette(rev(brewer.pal(n = 9, name = "RdYlBu")))(length(breaksList)-1), # Defines the vector of colors for the legend (it has to be of the same lenght of breaksList)
                               gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,main = "Usage Frequency",
                               annotation_col = annotation,
                               annotation_row = annotationrow,
                               annotation_colors = anno_colors)#gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
          }else if (nrow(OrderedTopLRpairsCellTypeConnectionMatrix_frequency) > 101 & nrow(OrderedTopLRpairsCellTypeConnectionMatrix_frequency) <= 151){
            if (SaveImage == TRUE) {
              pdf(file="Rplot%03d.pdf")
            }
            if (SaveImage == FALSE){
              #x11()
            }
            pheatmap::pheatmap(OrderedTopLRpairsCellTypeConnectionMatrix_frequency[1:50,],
                               cluster_rows = F, cluster_cols = F,
                               fontsize=7,fontsize_col=7, fontsize_row = 7,
                               legend_breaks = breaksList+1,
                               legend_labels = breaksList+1,
                               breaks = breaksList,
                               color = colorRampPalette(rev(brewer.pal(n = 9, name = "RdYlBu")))(length(breaksList)-1), # Defines the vector of colors for the legend (it has to be of the same lenght of breaksList)
                               main = "Usage Frequency",
                               gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
                               annotation_col = annotation,
                               annotation_row = annotationrow,
                               annotation_colors = anno_colors)#gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
            if (SaveImage == FALSE){
              #x11()
            }
            pheatmap::pheatmap(OrderedTopLRpairsCellTypeConnectionMatrix_frequency[51:100,],
                               cluster_rows = F, cluster_cols = F,
                               fontsize=7,fontsize_col=7, fontsize_row = 7,
                               legend_breaks = breaksList+1,
                               legend_labels = breaksList+1,
                               breaks = breaksList,
                               color = colorRampPalette(rev(brewer.pal(n = 9, name = "RdYlBu")))(length(breaksList)-1), # Defines the vector of colors for the legend (it has to be of the same lenght of breaksList)
                               main = "Usage Frequency",
                               gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
                               annotation_col = annotation,
                               annotation_row = annotationrow,
                               annotation_colors = anno_colors) #gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,

            if (SaveImage == FALSE){
              #x11()
            }
            pheatmap::pheatmap(OrderedTopLRpairsCellTypeConnectionMatrix_frequency[100:nrow(OrderedTopLRpairsCellTypeConnectionMatrix_frequency),],
                               cluster_rows = F, cluster_cols = F,
                               fontsize=7,fontsize_col=7, fontsize_row = 7,
                               legend_breaks = breaksList+1,
                               legend_labels = breaksList+1,
                               breaks = breaksList,main = "Usage Frequency",
                               color = colorRampPalette(rev(brewer.pal(n = 9, name = "RdYlBu")))(length(breaksList)-1), # Defines the vector of colors for the legend (it has to be of the same lenght of breaksList)
                               gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
                               annotation_col = annotation,
                               annotation_row = annotationrow,
                               annotation_colors = anno_colors)#gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
            # graphics.off()
          } else if (nrow(OrderedTopLRpairsCellTypeConnectionMatrix_frequency) > 151 & nrow(OrderedTopLRpairsCellTypeConnectionMatrix_frequency) <= 201){

            if (SaveImage == TRUE) {
              pdf(file="Rplot%03d.pdf")
            }
            if (SaveImage == FALSE){
              #x11()
            }
            pheatmap::pheatmap(OrderedTopLRpairsCellTypeConnectionMatrix_frequency[1:50,],
                               cluster_rows = F, cluster_cols = F,
                               fontsize=7,fontsize_col=7, fontsize_row = 7,

                               legend_breaks = breaksList+1,
                               legend_labels = breaksList+1,
                               breaks = breaksList,
                               color = colorRampPalette(rev(brewer.pal(n = 9, name = "RdYlBu")))(length(breaksList)-1), # Defines the vector of colors for the legend (it has to be of the same lenght of breaksList)
                               main = "Usage Frequency",
                               gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
                               annotation_col = annotation,
                               annotation_row = annotationrow,
                               annotation_colors = anno_colors)#gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
            if (SaveImage == FALSE){
              #x11()
            }
            pheatmap::pheatmap(OrderedTopLRpairsCellTypeConnectionMatrix_frequency[51:100,],
                               cluster_rows = F, cluster_cols = F,
                               fontsize=7,fontsize_col=7, fontsize_row = 7,

                               legend_breaks = breaksList+1,
                               legend_labels = breaksList+1,
                               breaks = breaksList,
                               color = colorRampPalette(rev(brewer.pal(n = 9, name = "RdYlBu")))(length(breaksList)-1), # Defines the vector of colors for the legend (it has to be of the same lenght of breaksList)
                               main = "Usage Frequency",
                               gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
                               annotation_col = annotation,
                               annotation_row = annotationrow,
                               annotation_colors = anno_colors) #gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
            if (SaveImage == FALSE){
              #x11()
            }
            pheatmap::pheatmap(OrderedTopLRpairsCellTypeConnectionMatrix_frequency[101:150,],
                               cluster_rows = F, cluster_cols = F,
                               fontsize=7,fontsize_col=7, fontsize_row = 7,
                               legend_breaks = breaksList+1,
                               legend_labels = breaksList+1,
                               breaks = breaksList,main = "Usage Frequency",
                               color = colorRampPalette(rev(brewer.pal(n = 9, name = "RdYlBu")))(length(breaksList)-1), # Defines the vector of colors for the legend (it has to be of the same lenght of breaksList)
                               gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
                               annotation_col = annotation,
                               annotation_row = annotationrow,
                               annotation_colors = anno_colors)#gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
            if (SaveImage == FALSE){
              #x11()
            }
            pheatmap::pheatmap(OrderedTopLRpairsCellTypeConnectionMatrix_frequency[151:nrow(OrderedTopLRpairsCellTypeConnectionMatrix_frequency),],
                               cluster_rows = F, cluster_cols = F,
                               fontsize=7,fontsize_col=7, fontsize_row = 7,
                               legend_breaks = breaksList+1,
                               legend_labels = breaksList+1,
                               breaks = breaksList,main = "Usage Frequency",
                               color = colorRampPalette(rev(brewer.pal(n = 9, name = "RdYlBu")))(length(breaksList)-1), # Defines the vector of colors for the legend (it has to be of the same lenght of breaksList)
                               gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
                               annotation_col = annotation,
                               annotation_row = annotationrow,
                               annotation_colors = anno_colors)#gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
            # graphics.off()
          }else if (nrow(OrderedTopLRpairsCellTypeConnectionMatrix_frequency) > 201){
            if (SaveImage == TRUE) {
              pdf(file="Rplot%03d.pdf")
            }
            if (SaveImage == FALSE){
              #x11()
            }
            pheatmap::pheatmap(OrderedTopLRpairsCellTypeConnectionMatrix_frequency[1:50,],
                               cluster_rows = F, cluster_cols = F,
                               fontsize=7,fontsize_col=7, fontsize_row = 7,
                               legend_breaks = breaksList+1,
                               legend_labels = breaksList+1,
                               breaks = breaksList,
                               color = colorRampPalette(rev(brewer.pal(n = 9, name = "RdYlBu")))(length(breaksList)-1), # Defines the vector of colors for the legend (it has to be of the same lenght of breaksList)
                               main = "Usage Frequency",
                               gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
                               annotation_col = annotation,
                               annotation_row = annotationrow,
                               annotation_colors = anno_colors)#gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
            if (SaveImage == FALSE){
              #x11()
            }
            pheatmap::pheatmap(OrderedTopLRpairsCellTypeConnectionMatrix_frequency[51:100,],
                               cluster_rows = F, cluster_cols = F,
                               fontsize=7,fontsize_col=7, fontsize_row = 7,
                               legend_breaks = breaksList+1,
                               legend_labels = breaksList+1,
                               breaks = breaksList,
                               color = colorRampPalette(rev(brewer.pal(n = 9, name = "RdYlBu")))(length(breaksList)-1), # Defines the vector of colors for the legend (it has to be of the same lenght of breaksList)
                               main = "Usage Frequency",
                               gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
                               annotation_col = annotation,
                               annotation_row = annotationrow,
                               annotation_colors = anno_colors) #gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,

            if (SaveImage == FALSE){
              #x11()
            }
            pheatmap::pheatmap(OrderedTopLRpairsCellTypeConnectionMatrix_frequency[101:150,],
                               cluster_rows = F, cluster_cols = F,
                               fontsize=7,fontsize_col=7, fontsize_row = 7,
                               legend_breaks = breaksList+1,
                               legend_labels = breaksList+1,
                               breaks = breaksList,main = "Usage Frequency",
                               color = colorRampPalette(rev(brewer.pal(n = 9, name = "RdYlBu")))(length(breaksList)-1), # Defines the vector of colors for the legend (it has to be of the same lenght of breaksList)
                               gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
                               annotation_col = annotation,
                               annotation_row = annotationrow,
                               annotation_colors = anno_colors)#gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
            if (SaveImage == FALSE){
              #x11()
            }
            pheatmap::pheatmap(OrderedTopLRpairsCellTypeConnectionMatrix_frequency[151:200,],
                               cluster_rows = F, cluster_cols = F,
                               fontsize=7,fontsize_col=7, fontsize_row = 7,
                               legend_breaks = breaksList+1,
                               legend_labels = breaksList+1,
                               breaks = breaksList,main = "Usage Frequency",
                               color = colorRampPalette(rev(brewer.pal(n = 9, name = "RdYlBu")))(length(breaksList)-1), # Defines the vector of colors for the legend (it has to be of the same lenght of breaksList)
                               gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
                               annotation_col = annotation,
                               annotation_row = annotationrow,
                               annotation_colors = anno_colors)#gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,

            if (SaveImage == FALSE){
              #x11()
            }
            pheatmap::pheatmap(OrderedTopLRpairsCellTypeConnectionMatrix_frequency[200:nrow(OrderedTopLRpairsCellTypeConnectionMatrix_frequency),],
                               cluster_rows = F, cluster_cols = F,
                               fontsize=7,fontsize_col=7, fontsize_row = 7,
                               legend_breaks = breaksList+1,
                               legend_labels = breaksList+1,
                               breaks = breaksList,main = "Usage Frequency",
                               color = colorRampPalette(rev(brewer.pal(n = 9, name = "RdYlBu")))(length(breaksList)-1), # Defines the vector of colors for the legend (it has to be of the same lenght of breaksList)
                               gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
                               annotation_col = annotation,
                               annotation_row = annotationrow,
                               annotation_colors = anno_colors)#gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
            # graphics.off()
          }
        }
      }


    } else if (IndexOfPlot == 1){

      OrderedTopLRpairsCellTypeConnectionMatrix[which(OrderedTopLRpairsCellTypeConnectionMatrix == 0)] <- NA
      # Com_CellTypeConnection_dd[which(Com_CellTypeConnection_dd == 0)] <- NA

      PlotMean_Shared <- function(input = OrderedTopLRpairsCellTypeConnectionMatrix){
        pheatmap::pheatmap(OrderedTopLRpairsCellTypeConnectionMatrix,
                           cluster_rows = F, cluster_cols = F,
                           color = colorRampPalette(c("#ffeda0","#fed976","#feb24c","#fd8d3c","#fc4e2a","#e31a1c","#bd0026"))(1000),#yellow to red
                           fontsize=7,fontsize_col=7, fontsize_row = 7,main = "Shared LRpairs",
                           gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
                           annotation_col = annotation,
                           annotation_row = annotationrow,
                           na_col = c("white"),
                           annotation_colors = anno_colors)#gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
      }

      PlotProportion_shared <-

      if (nrow(TopUsageLRpairs) == 0){
        print("There is no ligand-receptor pairs used more than 1 time.")
      } else if (nrow(TopUsageLRpairs) == 1){
        print("There is only 1 ligand-receptor pair that used more than 1 time.")
      } else if (nrow(TopUsageLRpairs) > 1){
        if (MultiFigure == FALSE){
          if (By_Proportion == FALSE){
            annotation <- data.frame(LR = factor(c(rep(0,ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2),
                                                   rep(1,ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2)),
                                                 labels = c("Ligand", "Receptor")))
            rownames(annotation) <- colnames(OrderedTopLRpairsCellTypeConnectionMatrix)

            LR <- c("navyblue", "red")
            names(LR) <- c("Ligand", "Receptor")
            anno_colors <- list(LR = LR)

            unique_labels <- sort(unique(FrequencyTopLRs))
            annotationrow <- data.frame(N = factor(c(FrequencyTopLRs),
                                                   labels = unique_labels))
            rownames(annotationrow) <- rownames(OrderedTopLRpairsCellTypeConnectionMatrix)

            if (SaveImage == TRUE) {
              pdf(file="Rplot%03d.pdf")
            }
            if (SaveImage == FALSE){
              #x11()
            }
            PlotMean_Shared(input = OrderedTopLRpairsCellTypeConnectionMatrix)
            # pheatmap::pheatmap(OrderedTopLRpairsCellTypeConnectionMatrix,
            #                    cluster_rows = F, cluster_cols = F,
            #                    color = colorRampPalette(c("#ffeda0","#fed976","#feb24c","#fd8d3c","#fc4e2a","#e31a1c","#bd0026"))(1000),#yellow to red
            #                    fontsize=7,fontsize_col=7, fontsize_row = 7,main = "Shared LRpairs",
            #                    gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
            #                    annotation_col = annotation,
            #                    annotation_row = annotationrow,
            #                    na_col = c("white"),
            #                    annotation_colors = anno_colors)#gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
          } else if (By_Proportion == TRUE){
            # Com_CellTypeConnection_dd[which(Com_CellTypeConnection_dd == 0)] <- NA
            annotation <- data.frame(CellType = factor(rep_vec, labels = c(Unique_CellType)),LR = factor(c(rep(0,ncol(Com_CellTypeConnection_dd)/2),
                                                                                                           rep(1,ncol(Com_CellTypeConnection_dd)/2)),
                                                                                                         labels = c("Ligand", "Receptor")))
            rownames(annotation) <- colnames(Com_CellTypeConnection_dd)


            LR <- c("navyblue", "red")
            names(LR) <- c("Ligand", "Receptor")
            anno_colors <- list(LR = LR)

            unique_labels <- sort(unique(FrequencyTopLRs))
            annotationrow <- data.frame(N = factor(c(FrequencyTopLRs),
                                                   labels = unique_labels))
            rownames(annotationrow) <- rownames(OrderedTopLRpairsCellTypeConnectionMatrix)
            if (SaveImage == TRUE) {
              pdf(file="CellProportion%03d.pdf")
            }
            if (SaveImage == FALSE){
              #x11()
            }
            pheatmap::pheatmap(Com_CellTypeConnection_dd,
                               cluster_rows = F, cluster_cols = F,
                               color = colorRampPalette(c("#ffffff","#000000","#007700","#00bb00","#33cc00","#669900","#996600","#cc3300","#ff0000"))(1000),

                               main = "Shared LRpairs",
                               fontsize=7,fontsize_col=7, fontsize_row = 7,
                               gaps_col = ncol(Com_CellTypeConnection_dd)/2,
                               show_colnames = F,
                               annotation_col = annotation,
                               annotation_row = annotationrow,
                               na_col = c("white"),
                               annotation_colors = anno_colors)
          }
        } else if (MultiFigure == TRUE){



          if (By_Proportion == FALSE){
            # library(pheatmap)
            annotation <- data.frame(LR = factor(c(rep(0,ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2),
                                                   rep(1,ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2)),
                                                 labels = c("Ligand", "Receptor")))
            rownames(annotation) <- colnames(OrderedTopLRpairsCellTypeConnectionMatrix)

            LR <- c("navyblue", "red")
            names(LR) <- c("Ligand", "Receptor")
            anno_colors <- list(LR = LR)

            unique_labels <- sort(unique(FrequencyTopLRs))
            annotationrow <- data.frame(N = factor(c(FrequencyTopLRs),
                                                   labels = unique_labels))
            rownames(annotationrow) <- rownames(OrderedTopLRpairsCellTypeConnectionMatrix)

            if (nrow(OrderedTopLRpairsCellTypeConnectionMatrix) <= 51) {
              if (SaveImage == TRUE) {
                pdf(file="Rplot%03d.pdf")
              }
              if (SaveImage == FALSE){
                #x11()
              }
              pheatmap::pheatmap(OrderedTopLRpairsCellTypeConnectionMatrix,
                                 cluster_rows = F, cluster_cols = F,
                                 color = colorRampPalette(c("#ffeda0","#fed976","#feb24c","#fd8d3c","#fc4e2a","#e31a1c","#bd0026"))(1000),#yellow to red
                                 fontsize=7,fontsize_col=7, fontsize_row = 7,main = "Shared LRpairs",
                                 gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
                                 annotation_col = annotation,
                                 annotation_row = annotationrow,
                                 na_col = c("white"),
                                 annotation_colors = anno_colors)#gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
            } else if (nrow(OrderedTopLRpairsCellTypeConnectionMatrix) > 51 & nrow(OrderedTopLRpairsCellTypeConnectionMatrix) <= 101){

              if (SaveImage == TRUE) {
                pdf(file="Rplot%03d.pdf")
              }

              if (SaveImage == FALSE){
                #x11()
              }
              pheatmap::pheatmap(OrderedTopLRpairsCellTypeConnectionMatrix[1:50,],
                                 cluster_rows = F, cluster_cols = F,
                                 color = colorRampPalette(c("#ffeda0","#fed976","#feb24c","#fd8d3c","#fc4e2a","#e31a1c","#bd0026"))(1000),#yellow to red

                                 fontsize=7,fontsize_col=7, fontsize_row = 7,main = "Shared LRpairs",
                                 gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
                                 annotation_col = annotation,
                                 annotation_row = annotationrow,
                                 na_col = c("white"),
                                 annotation_colors = anno_colors)#gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
              if (SaveImage == FALSE){
                #x11()
              }
              pheatmap::pheatmap(OrderedTopLRpairsCellTypeConnectionMatrix[51:nrow(OrderedTopLRpairsCellTypeConnectionMatrix),],
                                 cluster_rows = F, cluster_cols = F,
                                 color = colorRampPalette(c("#ffeda0","#fed976","#feb24c","#fd8d3c","#fc4e2a","#e31a1c","#bd0026"))(1000),#yellow to red
                                 fontsize=7,fontsize_col=7, fontsize_row = 7,main = "Shared LRpairs",
                                 gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
                                 annotation_col = annotation,
                                 annotation_row = annotationrow,
                                 na_col = c("white"),
                                 annotation_colors = anno_colors)#gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
            }else if (nrow(OrderedTopLRpairsCellTypeConnectionMatrix) > 101 & nrow(OrderedTopLRpairsCellTypeConnectionMatrix) <= 151){
              if (SaveImage == TRUE) {
                pdf(file="Rplot%03d.pdf")
              }
              if (SaveImage == FALSE){
                #x11()
              }
              pheatmap::pheatmap(OrderedTopLRpairsCellTypeConnectionMatrix[1:50,],
                                 cluster_rows = F, cluster_cols = F,
                                 color = colorRampPalette(c("#ffeda0","#fed976","#feb24c","#fd8d3c","#fc4e2a","#e31a1c","#bd0026"))(1000),#yellow to red
                                 fontsize=7,fontsize_col=7, fontsize_row = 7,main = "Shared LRpairs",
                                 gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
                                 annotation_col = annotation,
                                 annotation_row = annotationrow,
                                 na_col = c("white"),
                                 annotation_colors = anno_colors)#gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
              if (SaveImage == FALSE){
                #x11()
              }
              pheatmap::pheatmap(OrderedTopLRpairsCellTypeConnectionMatrix[51:100,],
                                 cluster_rows = F, cluster_cols = F,
                                 color = colorRampPalette(c("#ffeda0","#fed976","#feb24c","#fd8d3c","#fc4e2a","#e31a1c","#bd0026"))(1000),#yellow to red

                                 fontsize=7,fontsize_col=7, fontsize_row = 7,
                                 main = "Shared LRpairs",
                                 gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
                                 annotation_col = annotation,
                                 annotation_row = annotationrow,
                                 na_col = c("white"),
                                 annotation_colors = anno_colors)#gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
              if (SaveImage == FALSE){
                #x11()
              }
              pheatmap::pheatmap(OrderedTopLRpairsCellTypeConnectionMatrix[101:nrow(OrderedTopLRpairsCellTypeConnectionMatrix),],
                                 cluster_rows = F, cluster_cols = F,
                                 color = colorRampPalette(c("#ffeda0","#fed976","#feb24c","#fd8d3c","#fc4e2a","#e31a1c","#bd0026"))(1000),#yellow to red
                                 fontsize=7,fontsize_col=7, fontsize_row = 7,
                                 gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
                                 main = "Shared LRpairs",
                                 annotation_col = annotation,
                                 annotation_row = annotationrow,
                                 na_col = c("white"),
                                 annotation_colors = anno_colors)#gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,

            }else if (nrow(OrderedTopLRpairsCellTypeConnectionMatrix) > 151 & nrow(OrderedTopLRpairsCellTypeConnectionMatrix) <= 201){
              if (SaveImage == TRUE) {
                pdf(file="Rplot%03d.pdf")
              }
              if (SaveImage == FALSE){
                #x11()
              }
              pheatmap::pheatmap(OrderedTopLRpairsCellTypeConnectionMatrix[1:50,],
                                 cluster_rows = F, cluster_cols = F,
                                 color = colorRampPalette(c("#ffeda0","#fed976","#feb24c","#fd8d3c","#fc4e2a","#e31a1c","#bd0026"))(1000),#yellow to red
                                 fontsize=7,fontsize_col=7, fontsize_row = 7,main = "Shared LRpairs",
                                 gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
                                 annotation_col = annotation,
                                 annotation_row = annotationrow,na_col = c("white"),
                                 annotation_colors = anno_colors)#gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
              if (SaveImage == FALSE){
                #x11()
              }
              pheatmap::pheatmap(OrderedTopLRpairsCellTypeConnectionMatrix[51:100,],
                                 cluster_rows = F, cluster_cols = F,
                                 color = colorRampPalette(c("#ffeda0","#fed976","#feb24c","#fd8d3c","#fc4e2a","#e31a1c","#bd0026"))(1000),#yellow to red
                                 fontsize=7,fontsize_col=7, fontsize_row = 7,
                                 main = "Shared LRpairs",
                                 gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
                                 annotation_col = annotation,
                                 annotation_row = annotationrow,na_col = c("white"),
                                 annotation_colors = anno_colors)#gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
              if (SaveImage == FALSE){
                #x11()
              }
              pheatmap::pheatmap(OrderedTopLRpairsCellTypeConnectionMatrix[101:150,],
                                 cluster_rows = F, cluster_cols = F,
                                 color = colorRampPalette(c("#ffeda0","#fed976","#feb24c","#fd8d3c","#fc4e2a","#e31a1c","#bd0026"))(1000),#yellow to red
                                 fontsize=7,fontsize_col=7, fontsize_row = 7,
                                 main = "Shared LRpairs",
                                 gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
                                 annotation_col = annotation,
                                 annotation_row = annotationrow,na_col = c("white"),
                                 annotation_colors = anno_colors)#gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
              if (SaveImage == FALSE){
                #x11()
              }
              pheatmap::pheatmap(OrderedTopLRpairsCellTypeConnectionMatrix[151:nrow(OrderedTopLRpairsCellTypeConnectionMatrix),],
                                 cluster_rows = F, cluster_cols = F,
                                 color = colorRampPalette(c("#ffeda0","#fed976","#feb24c","#fd8d3c","#fc4e2a","#e31a1c","#bd0026"))(1000),#yellow to red
                                 fontsize=7,fontsize_col=7, fontsize_row = 7,
                                 gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
                                 main = "Shared LRpairs",
                                 annotation_col = annotation,na_col = c("white"),
                                 annotation_row = annotationrow,
                                 annotation_colors = anno_colors)#gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
            } else if (nrow(OrderedTopLRpairsCellTypeConnectionMatrix) > 200){
              if (SaveImage == TRUE) {
                pdf(file="Rplot%03d.pdf")
              }
              if (SaveImage == FALSE){
                #x11()
              }
              pheatmap::pheatmap(OrderedTopLRpairsCellTypeConnectionMatrix[1:50,],
                                 cluster_rows = F, cluster_cols = F,
                                 color = colorRampPalette(c("#ffeda0","#fed976","#feb24c","#fd8d3c","#fc4e2a","#e31a1c","#bd0026"))(1000),#yellow to red
                                 fontsize=7,fontsize_col=7, fontsize_row = 7,main = "Shared LRpairs",
                                 gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
                                 annotation_col = annotation,
                                 annotation_row = annotationrow,na_col = c("white"),
                                 annotation_colors = anno_colors)#gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
              if (SaveImage == FALSE){
                #x11()
              }
              pheatmap::pheatmap(OrderedTopLRpairsCellTypeConnectionMatrix[51:100,],
                                 cluster_rows = F, cluster_cols = F,
                                 color = colorRampPalette(c("#ffeda0","#fed976","#feb24c","#fd8d3c","#fc4e2a","#e31a1c","#bd0026"))(1000),#yellow to red
                                 fontsize=7,fontsize_col=7, fontsize_row = 7,
                                 main = "Shared LRpairs",
                                 gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
                                 annotation_col = annotation,
                                 annotation_row = annotationrow,na_col = c("white"),
                                 annotation_colors = anno_colors)#gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
              if (SaveImage == FALSE){
                #x11()
              }
              pheatmap::pheatmap(OrderedTopLRpairsCellTypeConnectionMatrix[101:150,],
                                 cluster_rows = F, cluster_cols = F,
                                 color = colorRampPalette(c("#ffeda0","#fed976","#feb24c","#fd8d3c","#fc4e2a","#e31a1c","#bd0026"))(1000),#yellow to red

                                 fontsize=7,fontsize_col=7, fontsize_row = 7,
                                 main = "Shared LRpairs",
                                 gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
                                 annotation_col = annotation,
                                 annotation_row = annotationrow,na_col = c("white"),
                                 annotation_colors = anno_colors)#gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
              if (SaveImage == FALSE){
                #x11()
              }
              pheatmap::pheatmap(OrderedTopLRpairsCellTypeConnectionMatrix[151:200,],
                                 cluster_rows = F, cluster_cols = F,
                                 color = colorRampPalette(c("#ffeda0","#fed976","#feb24c","#fd8d3c","#fc4e2a","#e31a1c","#bd0026"))(1000),#yellow to red

                                 fontsize=7,fontsize_col=7, fontsize_row = 7,
                                 main = "Shared LRpairs",
                                 gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
                                 annotation_col = annotation,
                                 annotation_row = annotationrow,na_col = c("white"),
                                 annotation_colors = anno_colors)#gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
              if (SaveImage == FALSE){
                #x11()
              }
              pheatmap::pheatmap(OrderedTopLRpairsCellTypeConnectionMatrix[201:nrow(OrderedTopLRpairsCellTypeConnectionMatrix),],
                                 cluster_rows = F, cluster_cols = F,
                                 color = colorRampPalette(c("#ffeda0","#fed976","#feb24c","#fd8d3c","#fc4e2a","#e31a1c","#bd0026"))(1000),#yellow to red

                                 fontsize=7,fontsize_col=7, fontsize_row = 7,
                                 gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
                                 main = "Shared LRpairs",
                                 annotation_col = annotation,
                                 annotation_row = annotationrow,
                                 na_col = c("white"),
                                 annotation_colors = anno_colors)#gaps_col =  ncol(OrderedTopLRpairsCellTypeConnectionMatrix)/2,
              # graphics.off()
            }
          } else if (By_Proportion == TRUE){

            annotation <- data.frame(CellType = factor(rep_vec, labels = c(Unique_CellType)),LR = factor(c(rep(0,ncol(Com_CellTypeConnection_dd)/2),
                                                                                                           rep(1,ncol(Com_CellTypeConnection_dd)/2)),
                                                                                                         labels = c("Ligand", "Receptor")))
            rownames(annotation) <- colnames(Com_CellTypeConnection_dd)

            LR <- c("navyblue", "red")
            names(LR) <- c("Ligand", "Receptor")
            anno_colors <- list(LR = LR)

            unique_labels <- sort(unique(FrequencyTopLRs))
            annotationrow <- data.frame(N = factor(c(FrequencyTopLRs),
                                                   labels = unique_labels))
            rownames(annotationrow) <- rownames(OrderedTopLRpairsCellTypeConnectionMatrix)

            if (nrow(Com_CellTypeConnection_dd)  <= 51 ){
              if (SaveImage == TRUE) {
                pdf(file="CellProportion%03d.pdf")
              }
              if (SaveImage == FALSE){
                #x11()
              }
              pheatmap::pheatmap(Com_CellTypeConnection_dd,
                                 cluster_rows = F, cluster_cols = F,

                                 color = colorRampPalette(c("#ffffff","#000000","#007700","#00bb00","#33cc00","#669900","#996600","#cc3300","#ff0000"))(1000),
                                 main = "Shared LRpairs",
                                 fontsize=7,fontsize_col=7, fontsize_row = 7,
                                 gaps_col =  ncol(Com_CellTypeConnection_dd)/2,
                                 show_colnames = F,
                                 annotation_col = annotation,
                                 annotation_row = annotationrow,
                                 annotation_colors = anno_colors)

            } else if (nrow(Com_CellTypeConnection_dd) > 51 & nrow(Com_CellTypeConnection_dd) <= 101){
              if (SaveImage == TRUE) {
                pdf(file="CellProportion%03d.pdf")
              }
              if (SaveImage == FALSE){
                #x11()
              }
              pheatmap::pheatmap(Com_CellTypeConnection_dd[1:50,],
                                 cluster_rows = F, cluster_cols = F,

                                 color = colorRampPalette(c("#ffffff","#000000","#007700","#00bb00","#33cc00","#669900","#996600","#cc3300","#ff0000"))(1000),
                                 main = "Shared LRpairs",
                                 fontsize=7,fontsize_col=7, fontsize_row = 7,
                                 gaps_col =  ncol(Com_CellTypeConnection_dd)/2,
                                 show_colnames = F,
                                 annotation_col = annotation,
                                 annotation_row = annotationrow,
                                 annotation_colors = anno_colors)
              if (SaveImage == FALSE){
                #x11()
              }
              pheatmap::pheatmap(Com_CellTypeConnection_dd[51:nrow(Com_CellTypeConnection_dd),],
                                 cluster_rows = F, cluster_cols = F,
                                 color = colorRampPalette(c("#ffffff","#000000","#007700","#00bb00","#33cc00","#669900","#996600","#cc3300","#ff0000"))(1000),
                                 main = "Shared LRpairs",
                                 fontsize=7,fontsize_col=7, fontsize_row = 7,
                                 gaps_col =  ncol(Com_CellTypeConnection_dd)/2,
                                 show_colnames = F,
                                 annotation_col = annotation,
                                 annotation_row = annotationrow,
                                 annotation_colors = anno_colors)

            }           else if (nrow(Com_CellTypeConnection_dd) > 151 &nrow(Com_CellTypeConnection_dd) < 201){
              if (SaveImage == TRUE) {
                pdf(file="CellProportion%03d.pdf")
              }
              if (SaveImage == FALSE){
                #x11()
              }
              pheatmap::pheatmap(Com_CellTypeConnection_dd[1:50,],
                                 cluster_rows = F, cluster_cols = F,

                                 color = colorRampPalette(c("#ffffff","#000000","#007700","#00bb00","#33cc00","#669900","#996600","#cc3300","#ff0000"))(1000),
                                 main = "Shared LRpairs",
                                 fontsize=7,fontsize_col=7, fontsize_row = 7,
                                 gaps_col =  ncol(Com_CellTypeConnection_dd)/2,
                                 show_colnames = F,
                                 annotation_col = annotation,
                                 annotation_row = annotationrow,
                                 annotation_colors = anno_colors)
              if (SaveImage == FALSE){
                #x11()
              }
              pheatmap::pheatmap(Com_CellTypeConnection_dd[51:100,],
                                 cluster_rows = F, cluster_cols = F,

                                 color = colorRampPalette(c("#ffffff","#000000","#007700","#00bb00","#33cc00","#669900","#996600","#cc3300","#ff0000"))(1000),
                                 main = "Shared LRpairs",
                                 fontsize=7,fontsize_col=7, fontsize_row = 7,
                                 gaps_col =  ncol(Com_CellTypeConnection_dd)/2,
                                 show_colnames = F,
                                 annotation_col = annotation,
                                 annotation_row = annotationrow,
                                 annotation_colors = anno_colors)
              if (SaveImage == FALSE){
                #x11()
              }
              pheatmap::pheatmap(Com_CellTypeConnection_dd[101:150,],
                                 cluster_rows = F, cluster_cols = F,

                                 color = colorRampPalette(c("#ffffff","#000000","#007700","#00bb00","#33cc00","#669900","#996600","#cc3300","#ff0000"))(1000),
                                 main = "Shared LRpairs",
                                 fontsize=7,fontsize_col=7, fontsize_row = 7,
                                 gaps_col =  ncol(Com_CellTypeConnection_dd)/2,
                                 show_colnames = F,
                                 annotation_col = annotation,
                                 annotation_row = annotationrow,
                                 annotation_colors = anno_colors)
              if (SaveImage == FALSE){
                #x11()
              }
              pheatmap::pheatmap(Com_CellTypeConnection_dd[151:nrow(Com_CellTypeConnection_dd),],
                                 cluster_rows = F, cluster_cols = F,

                                 color = colorRampPalette(c("#ffffff","#000000","#007700","#00bb00","#33cc00","#669900","#996600","#cc3300","#ff0000"))(1000),
                                 main = "Shared LRpairs",
                                 fontsize=7,fontsize_col=7, fontsize_row = 7,
                                 gaps_col =  ncol(Com_CellTypeConnection_dd)/2,
                                 show_colnames = F,
                                 annotation_col = annotation,
                                 annotation_row = annotationrow,
                                 annotation_colors = anno_colors)

            }  else if (nrow(Com_CellTypeConnection_dd) > 201){
              if (SaveImage == TRUE) {
                pdf(file="CellProportion%03d.pdf")
              }
              if (SaveImage == FALSE){
                #x11()
              }
              pheatmap::pheatmap(Com_CellTypeConnection_dd[1:50,],
                                 cluster_rows = F, cluster_cols = F,

                                 color = colorRampPalette(c("#ffffff","#000000","#007700","#00bb00","#33cc00","#669900","#996600","#cc3300","#ff0000"))(1000),
                                 main = "Shared LRpairs",
                                 fontsize=7,fontsize_col=7, fontsize_row = 7,
                                 gaps_col =  ncol(Com_CellTypeConnection_dd)/2,
                                 show_colnames = F,
                                 annotation_col = annotation,
                                 annotation_row = annotationrow,
                                 annotation_colors = anno_colors)
              if (SaveImage == FALSE){
                #x11()
              }
              pheatmap::pheatmap(Com_CellTypeConnection_dd[51:100,],
                                 cluster_rows = F, cluster_cols = F,
                                 color = colorRampPalette(c("#ffffff","#000000","#007700","#00bb00","#33cc00","#669900","#996600","#cc3300","#ff0000"))(1000),
                                 main = "Shared LRpairs",
                                 fontsize=7,fontsize_col=7, fontsize_row = 7,
                                 gaps_col =  ncol(Com_CellTypeConnection_dd)/2,
                                 show_colnames = F,
                                 annotation_col = annotation,
                                 annotation_row = annotationrow,
                                 annotation_colors = anno_colors)
              if (SaveImage == FALSE){
                #x11()
              }
              pheatmap::pheatmap(Com_CellTypeConnection_dd[101:150,],
                                 cluster_rows = F, cluster_cols = F,
                                 color = colorRampPalette(c("#ffffff","#000000","#007700","#00bb00","#33cc00","#669900","#996600","#cc3300","#ff0000"))(1000),
                                 main = "Shared LRpairs",
                                 fontsize=7,fontsize_col=7, fontsize_row = 7,
                                 gaps_col =  ncol(Com_CellTypeConnection_dd)/2,
                                 show_colnames = F,
                                 annotation_col = annotation,
                                 annotation_row = annotationrow,
                                 annotation_colors = anno_colors)
              if (SaveImage == FALSE){
                #x11()
              }
              pheatmap::pheatmap(Com_CellTypeConnection_dd[151:200,],
                                 cluster_rows = F, cluster_cols = F,
                                 color = colorRampPalette(c("#ffffff","#000000","#007700","#00bb00","#33cc00","#669900","#996600","#cc3300","#ff0000"))(1000),
                                 main = "Shared LRpairs",
                                 fontsize=7,fontsize_col=7, fontsize_row = 7,
                                 gaps_col =  ncol(Com_CellTypeConnection_dd)/2,
                                 show_colnames = F,
                                 annotation_col = annotation,
                                 annotation_row = annotationrow,
                                 annotation_colors = anno_colors)
              if (SaveImage == FALSE){
                #x11()
              }
              pheatmap::pheatmap(Com_CellTypeConnection_dd[201:nrow(Com_CellTypeConnection_dd),],
                                 cluster_rows = F, cluster_cols = F,
                                 color = colorRampPalette(c("#ffffff","#000000","#007700","#00bb00","#33cc00","#669900","#996600","#cc3300","#ff0000"))(1000),
                                 main = "Shared LRpairs",
                                 fontsize=7,fontsize_col=7, fontsize_row = 7,
                                 gaps_col =  ncol(Com_CellTypeConnection_dd)/2,
                                 show_colnames = F,
                                 annotation_col = annotation,
                                 annotation_row = annotationrow,
                                 annotation_colors = anno_colors)


            }
          }
        }
      }

    } else if (IndexOfPlot == 3){
      TopUniqueCommunication_RankByMean[which(TopUniqueCommunication_RankByMean== 0)] <- NA

      PlotHeat_Unique <- function(input = TopUniqueCommunication_RankByMean){
        pheatmap::pheatmap(input, cluster_rows = F, cluster_cols = F,
                           color = colorRampPalette(c("#ffeda0","#fed976","#feb24c","#fd8d3c","#fc4e2a","#e31a1c","#bd0026"))(1000),#yellow to red

                           fontsize=7,fontsize_col=7, fontsize_row = 7,
                           gaps_col =  ncol(input)/2,
                           main = "Unique LRpairs",
                           annotation = annotation,

                           annotation_colors = anno_colors,
                           na_col = c("white")
                           ) # gaps_col =  ncol(TopUniqueCommunication)/2,
      }

      if (nrow(UniqueLRpairs) == 0 ){
        print("There is no ligand-receptor pairs used only 1 time among cell types.")
      } else if (nrow(UniqueLRpairs) == 1){
        print("There is only 1 ligand-receptor pair that used only 1 time among cell types.")
      } else if (nrow(UniqueLRpairs) > 1){
        if (By_Proportion == FALSE){
          if (MultiFigure == FALSE){
            annotation <- data.frame(LR = factor(c(rep(0,ncol(TopUniqueCommunication_RankByMean)/2),rep(1,ncol(TopUniqueCommunication_RankByMean)/2)), labels = c("Ligand", "Receptor")))
            rownames(annotation) <- colnames(TopUniqueCommunication_RankByMean)
            LR <- c("navyblue", "red")
            names(LR) <- c("Ligand", "Receptor")
            anno_colors <- list(LR = LR)
            if (SaveImage == TRUE) {
              pdf(file="TopUniqueByMean%03d.pdf")
            }
            if (SaveImage == FALSE){
              #x11()
            }
            PlotHeat_Unique(input = TopUniqueCommunication_RankByMean)
            # gaps_col =  ncol(TopUniqueCommunication)/2,
          }else if (MultiFigure == TRUE){



            # if (nrow(OrderedUniqueLRpairsCellTypeConnectionMatrix_RankByMean) < NumberForTopUniq) {
            #   NumberForTopUniq <- nrow(OrderedUniqueLRpairsCellTypeConnectionMatrix_RankByMean)
            # }
            # TopUniqueCommunication_RankByMean <- OrderedUniqueLRpairsCellTypeConnectionMatrix_RankByMean[1:NumberForTopUniq,]
            # TopUniqueCommunication_RankByMean[which(TopUniqueCommunication_RankByMean== 0)] <- NA

            # TopUniqueCommunication_RankByMean <- OrderedUniqueLRpairsCellTypeConnectionMatrix_RankByMean[(NumberForTopUniq+1):nrow(OrderedUniqueLRpairsCellTypeConnectionMatrix_RankByMean),]
           # TopUniqueCommunication_RankByMean <- OrderedUniqueLRpairsCellTypeConnectionMatrix_RankByMean[(NumberForTopUniq+1):(NumberForTopUniq+50),]

            # TopUniqueCommunication_RankByMean[which(TopUniqueCommunication_RankByMean== 0)] <- NA

            annotation <- data.frame(LR = factor(c(rep(0,ncol(TopUniqueCommunication_RankByMean)/2),rep(1,ncol(TopUniqueCommunication_RankByMean)/2)), labels = c("Ligand", "Receptor")))
            rownames(annotation) <- colnames(TopUniqueCommunication_RankByMean)
            LR <- c("navyblue", "red")
            names(LR) <- c("Ligand", "Receptor")
            anno_colors <- list(LR = LR)

            if (SaveImage == TRUE) {
              pdf(file="TopUniqueByMean%03d.pdf")
            }
            if (SaveImage == FALSE){
              #x11()
            }
            if (nrow(TopUniqueCommunication_RankByMean) <= 50){
              PlotHeat_Unique(input = TopUniqueCommunication_RankByMean)
            } else if (nrow(TopUniqueCommunication_RankByMean) > 50 & nrow(TopUniqueCommunication_RankByMean) <=100){
              PlotHeat_Unique(input = TopUniqueCommunication_RankByMean[1:50,])
              PlotHeat_Unique(input = TopUniqueCommunication_RankByMean[51:nrow(TopUniqueCommunication_RankByMean),])
            } else if (nrow(TopUniqueCommunication_RankByMean)>100 & nrow(TopUniqueCommunication_RankByMean) <150){
              PlotHeat_Unique(input = TopUniqueCommunication_RankByMean[1:50,])
              PlotHeat_Unique(input = TopUniqueCommunication_RankByMean[51:100,])
              PlotHeat_Unique(input = TopUniqueCommunication_RankByMean[101:nrow(TopUniqueCommunication_RankByMean),])
            } else if (nrow(TopUniqueCommunication_RankByMean)>150 & nrow(TopUniqueCommunication_RankByMean) <200){
              PlotHeat_Unique(input = TopUniqueCommunication_RankByMean[1:50,])
              PlotHeat_Unique(input = TopUniqueCommunication_RankByMean[51:100,])
              PlotHeat_Unique(input = TopUniqueCommunication_RankByMean[101:150,])
              PlotHeat_Unique(input = TopUniqueCommunication_RankByMean[151:nrow(TopUniqueCommunication_RankByMean),])

            }else if (nrow(TopUniqueCommunication_RankByMean)>200){
              PlotHeat_Unique(input = TopUniqueCommunication_RankByMean[1:50,])
              PlotHeat_Unique(input = TopUniqueCommunication_RankByMean[51:100,])
              PlotHeat_Unique(input = TopUniqueCommunication_RankByMean[101:150,])
              PlotHeat_Unique(input = TopUniqueCommunication_RankByMean[151:200,])
              PlotHeat_Unique(input = TopUniqueCommunication_RankByMean[200:nrow(TopUniqueCommunication_RankByMean),])
            }

            # pheatmap::pheatmap(TopUniqueCommunication_RankByMean, cluster_rows = F, cluster_cols = F,
            #                    color = colorRampPalette(c("#ffeda0","#fed976","#feb24c","#fd8d3c","#fc4e2a","#e31a1c","#bd0026"))(1000),#yellow to red
            #                    fontsize=7,fontsize_col=7, fontsize_row = 7,
            #                    gaps_col =  ncol(TopUniqueCommunication_RankByMean)/2,
            #                    main = "Unique LRpairs",
            #                    annotation = annotation,na_col = c("white"),
            #                    annotation_colors = anno_colors) # gaps_col =  ncol(TopUniqueCommunication)/2,
          }
        } else if (By_Proportion == TRUE){

          PlotProportion <- function(input = Com_CellTypeConnection_dd){
            pheatmap::pheatmap(input,
                               cluster_rows = F, cluster_cols = F,

                               color = colorRampPalette(c("#ffffff","#000000","#007700","#00bb00","#44cc00","#77aa00","#997700","#cc5500","#ff0000"))(1000),
                               fontsize=7,fontsize_col=7, fontsize_row = 7,
                               main = "Unique LRpairs",
                               gaps_col =  ncol(input)/2,
                               show_colnames = F,
                               annotation_col = annotation,na_col = c("white"),
                               annotation_colors = anno_colors)


          }


          if (MultiFigure == FALSE){

            annotation <- data.frame(CellType = factor(rep_vec, labels = c(Unique_CellType)),LR = factor(c(rep(0,ncol(Com_CellTypeConnection_dd)/2),
                                                                                                           rep(1,ncol(Com_CellTypeConnection_dd)/2)),
                                                                                                         labels = c("Ligand", "Receptor")))
            rownames(annotation) <- colnames(Com_CellTypeConnection_dd)

            LR <- c("navyblue", "red")

            names(LR) <- c("Ligand", "Receptor")
            anno_colors <- list(LR = LR)

            if (SaveImage == TRUE) {
              pdf(file="CellProportion%03d.pdf")
            }
            if (SaveImage == FALSE){
              #x11()
            }
            PlotProportion(input = Com_CellTypeConnection_dd)

          } else if (MultiFigure == TRUE){

            annotation <- data.frame(CellType = factor(rep_vec, labels = c(Unique_CellType)),LR = factor(c(rep(0,ncol(Com_CellTypeConnection_dd)/2),
                                                                                                           rep(1,ncol(Com_CellTypeConnection_dd)/2)),
                                                                                                         labels = c("Ligand", "Receptor")))
            rownames(annotation) <- colnames(Com_CellTypeConnection_dd)

            LR <- c("navyblue", "red")

            names(LR) <- c("Ligand", "Receptor")
            anno_colors <- list(LR = LR)

            if (SaveImage == TRUE) {
              pdf(file="CellProportion%03d.pdf")
            }
            if (SaveImage == FALSE){
              #x11()
            }

            if (nrow(Com_CellTypeConnection_dd) <= 50){
              PlotProportion(input = Com_CellTypeConnection_dd)
            } else if (nrow(Com_CellTypeConnection_dd) > 50 & nrow(Com_CellTypeConnection_dd) <=100){
              PlotProportion(input = Com_CellTypeConnection_dd[1:50,])
              PlotProportion(input = Com_CellTypeConnection_dd[51:nrow(Com_CellTypeConnection_dd),])
            } else if (nrow(Com_CellTypeConnection_dd)>100 & nrow(Com_CellTypeConnection_dd) <150){
              PlotProportion(input = Com_CellTypeConnection_dd[1:50,])
              PlotProportion(input = Com_CellTypeConnection_dd[51:100,])
              PlotProportion(input = Com_CellTypeConnection_dd[101:nrow(Com_CellTypeConnection_dd),])
            } else if (nrow(Com_CellTypeConnection_dd)>150 & nrow(Com_CellTypeConnection_dd) <200){
              PlotProportion(input = Com_CellTypeConnection_dd[1:50,])
              PlotProportion(input = Com_CellTypeConnection_dd[51:100,])
              PlotProportion(input = Com_CellTypeConnection_dd[101:150,])
              PlotProportion(input = Com_CellTypeConnection_dd[151:nrow(Com_CellTypeConnection_dd),])

            }else if (nrow(Com_CellTypeConnection_dd)>200){
              PlotProportion(input = Com_CellTypeConnection_dd[1:50,])
              PlotProportion(input = Com_CellTypeConnection_dd[51:100,])
              PlotProportion(input = Com_CellTypeConnection_dd[101:150,])
              PlotProportion(input = Com_CellTypeConnection_dd[151:200,])
              # PlotHeat_Unique(input = Com_CellTypeConnection_dd[200:nrow(Com_CellTypeConnection_dd),])
            }


            # if (nrow(Com_CellTypeConnection_dd) < NumberForTopUniq){
            #   NumberForTopUniq <- nrow(Com_CellTypeConnection_dd)
            # }
            # pheatmap::pheatmap(Com_CellTypeConnection_dd[1:NumberForTopUniq,],
            #                    cluster_rows = F, cluster_cols = F,
            #
            #                    color = colorRampPalette(c("#ffffff","#000000","#007700","#00bb00","#44cc00","#77aa00","#997700","#cc5500","#ff0000"))(1000),
            #                    fontsize=7,fontsize_col=7, fontsize_row = 7,
            #                    main = "Unique LRpairs",
            #                    gaps_col =  ncol(Com_CellTypeConnection_dd)/2,
            #                    show_colnames = F,
            #                    annotation_col = annotation,na_col = c("white"),
            #                    annotation_colors = anno_colors)
          }
        }

      }
    }
    if (SaveImage == TRUE){
      dev.off ()
    }
  }

  # names(CelltypeConnection_List) <- c("SharedLRpairs","UniqueLRpairs","FrequencyOfSharedLRpairs","TotalLRpairs_RankByMean")

  if(ReturnList == TRUE){
    return(CelltypeConnection_List)
  }

}

#' Visualize the connections for a specific cell type.
#'
#' This function extracts and plots connections of a specific cell type with other cell types.
#' For a certain cell type, ligand to receptor is defined as OUT signal and receptor to ligand is defined as IN signal.
#' By setting parameters, users can chose the type of signals to view.
#' The ligand-receptor pairs are ranked by their mean expression value.
#'
#' @param LRC3_list the input data generated by function LRC3_INF.
#' @param CellType_Index the Index of the queried cell type that users are interested. The maximum number is the number of cell types.
#' @param By_Proportion if \code{By_Proportion = FALSE}, view the connections based on the mean expression of cell types.
#' if \code{By_Proportion = TRUE}, view the connections based on the proportion of cells within a cell type that significantly express the specific ligand or receptor.
#' @param Combine_InOut If \code{Combine_InOut = TRUE}, it gives a figure that combines the IN and OUT signal tegether.
#' @param TargetCellType_IN if \code{TargetCellType_IN = TRUE},it shows the IN signal for the cell type in query.
#' @param TargetCellType_OUT if \code{TargetCellType_OUT = TRUE}, it shows the OUT signal for the cell type in query.
#' @param ViewTop_LRpairs if \code{ViewTop_LRpairs = FALSE}, show all the ligand-receptor pairs in the output figure.
#' if \code{ViewTop_UniqLRpairs = TRUE}, show a number of top highly expressed ligand-receptor pairs in the output figure.
#' The ligand-receptor pairs are ranked by their mean expression value in an decreasing order.
#' @param TOP_Number the number of ligand-receptor pairs to show if \code{ViewTop_LRpairs = TRUE}.
#' @param SaveImage if \code{SaveImage = FALSE},plot without saving the figure. if \code{SaveImage = TRUE},save the figure as pdf.
#' @importFrom pheatmap pheatmap
#' @export


LRC3_CellTypeConnection <- function(LRC3_list,CellType_Index = 1,By_Proportion = FALSE,
                                    Combine_InOut = FALSE,TargetCellType_IN = FALSE,TargetCellType_OUT = FALSE,
                                    ViewTop_LRpairs = FALSE, TOP_Number = 50,SaveImage = FALSE){
  LRpairsTotalUsage <- NULL
  ScaledLRpairsTotalUsage <- NULL
  LigandExpressionLevelINF <- LRC3_list$LigandExpressionLevelINF

  GroupedUniqueLRMatrix <- LRC3_list$GroupedUniqueLRMatrix
  p <- LRC3_list$p
  k <- LRC3_list$k

  CellTypeSpecific_IN <- list()
  CellTypeSpecific_OUT <- list()
  CellTypeSpecific_InOut <- list()
  NoneZero_IN_list <- list()
  NoneZero_OUT_list <- list()
  NoneZero_Total_list <- list()
  List_NUM <- 0
  for (i in 1:nrow(LigandExpressionLevelINF[1,,])){# cell type 1 express ligand and all the other cell type express receptor,for all the LRpairs
    Celltype_SpecificLRpairsTotalUsage_OUT <- apply(LRC3_list$Remaining_contactINF[i,,],2,sum)
    nameVec <- paste0(LRC3_list$rem_LRpairs$Ligand.ApprovedSymbol,"--",LRC3_list$rem_LRpairs$Receptor.ApprovedSymbol)
    names(Celltype_SpecificLRpairsTotalUsage_OUT) <- nameVec
    Celltype_SpecificLRpairsTotalUsage_IN <- apply(LRC3_list$Remaining_contactINF[,i,],2,sum)
    names(Celltype_SpecificLRpairsTotalUsage_IN) <- nameVec

    Celltype_SpecificLRpairsTotalUsage_OUT_Fortotal <- apply(LRC3_list$Remaining_contactINF[i,,],2,function(x){sum(x)-x[i]})
    names(Celltype_SpecificLRpairsTotalUsage_OUT_Fortotal) <- nameVec

    Celltype_SpecificLRpairsTotalUsage_Total <- Celltype_SpecificLRpairsTotalUsage_OUT_Fortotal + Celltype_SpecificLRpairsTotalUsage_IN

    LenNonZeros_IN <- length(which(sort(Celltype_SpecificLRpairsTotalUsage_IN,decreasing = T) > 0))
    TableOfFrequency_IN <- table(Celltype_SpecificLRpairsTotalUsage_IN)
    OrderedIndex_IN <- order(Celltype_SpecificLRpairsTotalUsage_IN,decreasing = T)[1:LenNonZeros_IN]
    NoneZero_IN <- Celltype_SpecificLRpairsTotalUsage_IN[OrderedIndex_IN]
    Ligand_IN_INF <- LRC3_list$ScaledLigandExpressionLevelINF[,,OrderedIndex_IN]
    Receptor_IN_INF <- LRC3_list$ScaledReceptorExpressionLevelINF[,,OrderedIndex_IN]

    if (LenNonZeros_IN == 0){
      IN_Matrix_CellTypeSpecific <- 0
      NoneZero_IN <- 0
    } else if (LenNonZeros_IN ==1 ){
       IN_Matrix_CellTypeSpecific <- 0
       NoneZero_IN <- 0

    } else if (LenNonZeros_IN > 1){
      Ligand_IN <- t(Ligand_IN_INF[,i,])
      colnames_IN <- colnames(LRC3_list$contact_number_Remaining)
      colnames(Ligand_IN) <- paste0(colnames_IN,"-L")
      Target_R <- apply(t(Receptor_IN_INF[,i,]),1,function(x){x[which(x !=0)][1]})
      Target_R[is.na(Target_R)] <- 0

      REM_LRpairs_IN <- LRC3_list$rem_LRpairs[OrderedIndex_IN,]
      REM_rownames_IN <- paste0(REM_LRpairs_IN$Ligand.ApprovedSymbol,"--",REM_LRpairs_IN$Receptor.ApprovedSymbol)

      IN_Matrix_CellTypeSpecific <- cbind(Target_R,Ligand_IN)
      colnames(IN_Matrix_CellTypeSpecific)[1] <- paste0(colnames_IN[CellType_Index],"_R")

      rownames(IN_Matrix_CellTypeSpecific) <- REM_rownames_IN


      Ordered_IN_Matrix_CellTypeSpecific <- NULL
      NameForOne <- NULL
      for (fi in max(NoneZero_IN):min(NoneZero_IN)){
        sub_Matrix <- IN_Matrix_CellTypeSpecific[which(NoneZero_IN == fi),]
        if (length(which(NoneZero_IN ==fi)) != 1){
          vec_NonezeroNum <- apply(sub_Matrix,1,function(x){length(which(x != 0))})
          sub_Matrix <- sub_Matrix[order(rowSums(sub_Matrix)/vec_NonezeroNum,decreasing = T),]
        }
        Ordered_IN_Matrix_CellTypeSpecific <- rbind(Ordered_IN_Matrix_CellTypeSpecific,sub_Matrix)
        if (length(which(NoneZero_IN == fi)) == 1){
          NameForOne <- cbind(NameForOne,names(which(NoneZero_IN == fi)))
        }
      }
      rownames(Ordered_IN_Matrix_CellTypeSpecific)[which(rownames(Ordered_IN_Matrix_CellTypeSpecific)== "sub_Matrix")] <- c(NameForOne)
      IN_Matrix_CellTypeSpecific <- Ordered_IN_Matrix_CellTypeSpecific
    }

    LenNonZeros_OUT <- length(which(sort(Celltype_SpecificLRpairsTotalUsage_OUT,decreasing = T) > 0))
    if (LenNonZeros_OUT == 0){
      OUT_Matrix_CellTypeSpecific <- 0
      NoneZero_OUT <- 0
    } else if (LenNonZeros_OUT == 1){
     OUT_Matrix_CellTypeSpecific <- 0
     NoneZero_OUT <- 0
     }else if (LenNonZeros_OUT > 1){

        TableOfFrequency_OUT <- table(Celltype_SpecificLRpairsTotalUsage_OUT)
        OrderedIndex_OUT <- order(Celltype_SpecificLRpairsTotalUsage_OUT,decreasing = T)[1:LenNonZeros_OUT]
        NoneZero_OUT <- Celltype_SpecificLRpairsTotalUsage_OUT[OrderedIndex_OUT]
        Ligand_OUT_INF <- LRC3_list$ScaledLigandExpressionLevelINF[,,OrderedIndex_OUT]
        Receptor_OUT_INF <- LRC3_list$ScaledReceptorExpressionLevelINF[,,OrderedIndex_OUT]

        REM_LRpairs_OUT <- LRC3_list$rem_LRpairs[OrderedIndex_OUT,]
        REM_rownames_OUT <- paste0(REM_LRpairs_OUT$Ligand.ApprovedSymbol,"--",REM_LRpairs_OUT$Receptor.ApprovedSymbol)

        Receptor_OUT <- t(Receptor_OUT_INF[i,,])
        colnames_OUT <- colnames(LRC3_list$contact_number_Remaining)
        colnames(Receptor_OUT) <- paste0(colnames_OUT,"-R")
        Target_L <- apply(t(Ligand_OUT_INF[i,,]),1,function(x){x[which(x !=0)][1]})
        Target_L[is.na(Target_L)] <- 0
        OUT_Matrix_CellTypeSpecific <- cbind(Target_L,Receptor_OUT)
        colnames(OUT_Matrix_CellTypeSpecific)[1] <- paste0(colnames_OUT[CellType_Index],"_L")
        rownames(OUT_Matrix_CellTypeSpecific) <- REM_rownames_OUT
        Ordered_OUT_Matrix_CellTypeSpecific <- NULL
        NameForOne <- NULL
        for (fi in max(NoneZero_OUT):min(NoneZero_OUT)){
          sub_Matrix <- OUT_Matrix_CellTypeSpecific[which(NoneZero_OUT == fi),]
          if (length(which(NoneZero_OUT ==fi)) != 1){
            vec_NonezeroNum <- apply(sub_Matrix,1,function(x){length(which(x != 0))})
            sub_Matrix <- sub_Matrix[order(rowSums(sub_Matrix)/vec_NonezeroNum,decreasing = T),]
          }
          Ordered_OUT_Matrix_CellTypeSpecific <- rbind(Ordered_OUT_Matrix_CellTypeSpecific,sub_Matrix)
          if (length(which(NoneZero_OUT ==fi)) == 1){
            NameForOne <- cbind(NameForOne,names(which(NoneZero_OUT == fi)))
          }
        }
        rownames(Ordered_OUT_Matrix_CellTypeSpecific)[which(rownames(Ordered_OUT_Matrix_CellTypeSpecific)== "sub_Matrix")] <- c(NameForOne)
        OUT_Matrix_CellTypeSpecific <- Ordered_OUT_Matrix_CellTypeSpecific
      }

    LenNonZeros_total <- length(which(sort(Celltype_SpecificLRpairsTotalUsage_Total,decreasing = T) > 0))
    TableOfFrequency_total <- table(Celltype_SpecificLRpairsTotalUsage_Total)
    OrderedIndex_total <- order(Celltype_SpecificLRpairsTotalUsage_Total,decreasing = T)[1:LenNonZeros_total]
    NoneZero_Total <- Celltype_SpecificLRpairsTotalUsage_Total[OrderedIndex_total]
    Ligand_Total_INF <- LRC3_list$ScaledLigandExpressionLevelINF[,,OrderedIndex_total]
    Receptor_Total_INF <- LRC3_list$ScaledReceptorExpressionLevelINF[,,OrderedIndex_total]

    REM_LRpairs <- LRC3_list$rem_LRpairs[OrderedIndex_total,]
    REM_rownames <- paste0(REM_LRpairs$Ligand.ApprovedSymbol,"--",REM_LRpairs$Receptor.ApprovedSymbol)

    if (LenNonZeros_total == 0){ #in case there is no signal out for this cell type.
      TotalMatrix_CellTypeSpecific <- 0
      NoneZero_Total <- 0
    } else if (LenNonZeros_total !=0){
      if (LenNonZeros_total == 1){
        TotalMatrix_CellTypeSpecific <- 0
        NoneZero_Total <- 0
      } else if (LenNonZeros_total > 1){
        Ligand_IN <- t(Ligand_Total_INF[,i,])
        colnames(Ligand_IN) <- paste0(colnames(LRC3_list$contact_number_Remaining),"-L")

        Receptor_OUT <- t(Receptor_Total_INF[i,,])
        colnames(Receptor_OUT) <- paste0(colnames(LRC3_list$contact_number_Remaining),"-R")
        Target_L <- apply(t(Ligand_Total_INF[i,,]),1,function(x){x[which(x !=0)][1]})
        Target_L[is.na(Target_L)] <- 0

        Target_R <- apply(t(Receptor_Total_INF[,i,]),1,function(x){x[which(x !=0)][1]})
        Target_R[is.na(Target_R)] <- 0

        TotalMatrix_CellTypeSpecific <- cbind(Target_L,Receptor_OUT,Ligand_IN,Target_R)

        colnames_targets <- c(paste0(sapply(strsplit(colnames(Ligand_IN)[CellType_Index],"-"),"[[",1),"_L"),paste0(sapply(strsplit(colnames(Ligand_IN)[CellType_Index],"-"),"[[",1),"_R"))
        colnames(TotalMatrix_CellTypeSpecific)[c(1,ncol(TotalMatrix_CellTypeSpecific))] <- colnames_targets
        rownames(TotalMatrix_CellTypeSpecific) <- REM_rownames

        Ordered_TotalCellTypeSpecific <- NULL
        Ordered_TopCommonLRpairsConnection_Frequency <- NULL
        NameForOne <- NULL
        for (fi in max(NoneZero_Total):min(NoneZero_Total)){
          sub_Matrix <- TotalMatrix_CellTypeSpecific[which(NoneZero_Total == fi),]
          if (length(which(NoneZero_Total ==fi)) != 1){
            vec_NonezeroNum <- apply(sub_Matrix,1,function(x){length(which(x != 0))})
            sub_Matrix <- sub_Matrix[order(rowSums(sub_Matrix)/vec_NonezeroNum,decreasing = T),]
          }
          Ordered_TotalCellTypeSpecific <- rbind(Ordered_TotalCellTypeSpecific,sub_Matrix)
          if (length(which(NoneZero_Total == fi)) == 1){
            NameForOne <- cbind(NameForOne,names(which(NoneZero_Total == fi)))
          }
        }
        rownames(Ordered_TotalCellTypeSpecific)[which(rownames(Ordered_TotalCellTypeSpecific)== "sub_Matrix")] <- c(NameForOne)
        TotalMatrix_CellTypeSpecific <- Ordered_TotalCellTypeSpecific
      }
    }


    CellTypeSpecific_IN[[i]] <- IN_Matrix_CellTypeSpecific
    NoneZero_IN_list[[i]] <- NoneZero_IN
    CellTypeSpecific_OUT[[i]] <- OUT_Matrix_CellTypeSpecific
    NoneZero_OUT_list[[i]] <- NoneZero_OUT
    CellTypeSpecific_InOut[[i]] <- TotalMatrix_CellTypeSpecific
    NoneZero_Total_list[[i]] <- NoneZero_Total
  }

  M_IN <- CellTypeSpecific_IN[[CellType_Index]]
  M_OUT <- CellTypeSpecific_OUT[[CellType_Index]]
  M_InOut <- CellTypeSpecific_InOut[[CellType_Index]]
  Index_IN <- NoneZero_IN_list[[CellType_Index]]
  Index_OUT<- NoneZero_OUT_list[[CellType_Index]]
  Index_InOut<-NoneZero_Total_list[[CellType_Index]]

###################################################################

  if (By_Proportion == TRUE){
    if (TargetCellType_IN == TRUE){
      TopMatrix_IN <- M_IN
      if (sum(TopMatrix_IN) == 0) {
        CellType_IN <- 0
      }else if (sum(TopMatrix_IN) != 0){

        rem_LRpairs <- rownames(TopMatrix_IN)
        rem_ComL <- sapply(strsplit(rem_LRpairs,"--"),"[[",1)
        rem_ComR <- sapply(strsplit(rem_LRpairs,"--"),"[[",2)
        ComReceptor_Matrix <- TopMatrix_IN[,1]
        ComLigand_Matrix <- TopMatrix_IN[,2:ncol(TopMatrix_IN)]

        ComLigand_Binary <- ComLigand_Matrix
        ComLigand_Binary[which(ComLigand_Binary != 0)] <- 1
        ComReceptor_Binary <- ComReceptor_Matrix
        ComReceptor_Binary[which(ComReceptor_Binary !=0)] <- 1

        index_ComLigand <- sapply(as.list(rem_ComL),function(x){which(rownames(GroupedUniqueLRMatrix) == x)})
        index_ComReceptor <- sapply(as.list(rem_ComR),function(x){which(rownames(GroupedUniqueLRMatrix) == x)})

        ComReceptor_dd <- GroupedUniqueLRMatrix[index_ComReceptor, which(colnames(GroupedUniqueLRMatrix) ==CellType_Index)]
        ComLigand_dd <- GroupedUniqueLRMatrix[index_ComLigand,]

        if (is.null(ncol(ComReceptor_dd))){
          gap_In <- length(ComReceptor_dd)
          ComReceptor_colnames <- CellType_Index
          ComLigand_colnames <- sapply(strsplit(names(ComLigand_Binary),"-"),"[[",1)
          Index_P <- sapply(as.list(rem_ComL),function(x){which(names(p) == x)}[1])
          Index_K <- sapply(as.list(rem_ComR),function(x){which(names(k) == x)}[1])

          rem_P <- p[Index_P]
          rem_K <- k[Index_K]

          ComReceptor_dd[which(ComReceptor_dd < rem_K)] <- 0
          ComReceptor_dd <- ComReceptor_dd/rem_K
          ComReceptor_dd[which(ComReceptor_dd > 7)] <- 7

          TargetReceptor_dd_IN <- ComReceptor_dd
          ComLigand_dd_Connection_IN <- NULL
          for (i in 1:length(ComLigand_colnames)){

            sub_ComLigand_dd <- ComLigand_dd[which(names(ComLigand_dd)== i)]
            sub_ComLigand_Binary <- ComLigand_Binary[which(ComLigand_colnames == ComLigand_colnames[i])]

            tem_ComLigand_dd <- sub_ComLigand_dd*sub_ComLigand_Binary
            tem_ComLigand_dd[which(tem_ComLigand_dd < rem_P)] <- 0
            tem_ComLigand_dd <- tem_ComLigand_dd/rem_P
            tem_ComLigand_dd[which(tem_ComLigand_dd > 7)] <- 7

            names(tem_ComLigand_dd) <- paste0(names(tem_ComLigand_dd),".",c(1:length(tem_ComLigand_dd)))
            ComLigand_dd_Connection_IN <- c(ComLigand_dd_Connection_IN,tem_ComLigand_dd)

          }
          names(TargetReceptor_dd_IN) <- paste0(names(TargetReceptor_dd_IN),".",c(1:length(TargetReceptor_dd_IN)))
          CellType_IN <- c(TargetReceptor_dd_IN,ComLigand_dd_Connection_IN)

          names(CellType_IN)[1:length(TargetReceptor_dd_IN)] <- paste0(names(TargetReceptor_dd_IN),"-R")
          names(CellType_IN)[(length(TargetReceptor_dd_IN)+1):length(CellType_IN)] <- paste0(names(CellType_IN)[(length(TargetReceptor_dd_IN)+1):length(CellType_IN)],"-L")

          ####################################################
          Unique_CellType <- sort(unique(names(ComLigand_dd)))

          rep_Num_IN <- NULL
          rep_vec_IN <- NULL
          for (i in 1:length(Unique_CellType)){
            rep_Num_IN[i] <- length(which(names(ComLigand_dd) == i))
            rep_subvec <- rep(i-1,rep_Num_IN[i])
            rep_vec_IN <- c(rep_vec_IN,rep_subvec)

          }

          Target_rep <- rep(CellType_Index-1,length(ComReceptor_dd))
          rep_vec_IN <- c(Target_rep,rep_vec_IN)#for the annotation of pheatmap

        } else if (!is.null(ncol(ComReceptor_dd))) {
          gap_In <- ncol(ComReceptor_dd)

          ComReceptor_colnames <- CellType_Index
          ComLigand_colnames <- sapply(strsplit(colnames(ComLigand_Binary),"-"),"[[",1)
          Index_P <- sapply(as.list(rem_ComL),function(x){which(names(p) == x)}[1])#get the positions of remaining p
          Index_K <- sapply(as.list(rem_ComR),function(x){which(names(k) == x)}[1])

          rem_P <- p[Index_P]
          rem_K <- k[Index_K]

          ComReceptor_dd[which(ComReceptor_dd < rem_K)] <- 0
          ComReceptor_dd <- ComReceptor_dd/rem_K
          ComReceptor_dd[which(ComReceptor_dd > 7)] <- 7

          TargetReceptor_dd_IN <- ComReceptor_dd
          ComLigand_dd_Connection_IN <- NULL


          for (i in 1:length(ComLigand_colnames)){

            sub_ComLigand_dd <- ComLigand_dd[,which(colnames(ComLigand_dd)== i)]
            sub_ComLigand_Binary <- ComLigand_Binary[,which(ComLigand_colnames == ComLigand_colnames[i])]

            tem_ComLigand_dd <- sub_ComLigand_dd*sub_ComLigand_Binary
            tem_ComLigand_dd[which(tem_ComLigand_dd < rem_P)] <- 0
            tem_ComLigand_dd <- tem_ComLigand_dd/rem_P
            tem_ComLigand_dd[which(tem_ComLigand_dd > 7)] <- 7
            colnames(tem_ComLigand_dd) <- paste0(colnames(tem_ComLigand_dd),".",c(1:ncol(tem_ComLigand_dd)))
            ComLigand_dd_Connection_IN <- cbind(ComLigand_dd_Connection_IN,tem_ComLigand_dd)
            colnames(TargetReceptor_dd_IN) <- paste0(colnames(TargetReceptor_dd_IN),".",c(1:ncol(TargetReceptor_dd_IN)))
            CellType_IN <- cbind(TargetReceptor_dd_IN,ComLigand_dd_Connection_IN)

          }
          rownames(CellType_IN) <- paste0(rownames(ComLigand_dd_Connection_IN),"--",rownames(TargetReceptor_dd_IN))
          colnames(CellType_IN)[1:ncol(TargetReceptor_dd_IN)] <- paste0(colnames(TargetReceptor_dd_IN),"-R")
          colnames(CellType_IN)[(ncol(TargetReceptor_dd_IN)+1):ncol(CellType_IN)] <- paste0(colnames(CellType_IN)[(ncol(TargetReceptor_dd_IN)+1):ncol(CellType_IN)],"-L")
          Unique_CellType <- sort(unique(colnames(ComLigand_dd)))

          rep_Num_IN <- NULL
          rep_vec_IN <- NULL
          for (i in 1:length(Unique_CellType)){
            rep_Num_IN[i] <- length(which(colnames(ComLigand_dd) == i))
            rep_subvec <- rep(i-1,rep_Num_IN[i])
            rep_vec_IN <- c(rep_vec_IN,rep_subvec)

          }

          Target_rep <- rep(CellType_Index-1,ncol(ComReceptor_dd))
          rep_vec_IN <- c(Target_rep,rep_vec_IN)#for the annotation of pheatmap

        }
      }
    }

    if (TargetCellType_OUT == TRUE){
      TopMatrix_OUT <- M_OUT

      if (sum(TopMatrix_OUT) == 0) {
        CellType_OUT <- 0
      }else if (sum(TopMatrix_OUT) != 0){
        if (is.null(dim(TopMatrix_OUT))){

          rem_LRpairs <- rownames(TopMatrix_OUT)
          rem_ComL <- sapply(strsplit(rem_LRpairs,"--"),"[[",1)
          rem_ComR <- sapply(strsplit(rem_LRpairs,"--"),"[[",2)
          ComLigand_Matrix <- TopMatrix_OUT[,1]
          ComReceptor_Matrix <- TopMatrix_OUT[,2:ncol(TopMatrix_OUT)]

          ComLigand_Binary <- ComLigand_Matrix
          ComLigand_Binary[which(ComLigand_Binary != 0)] <- 1
          ComReceptor_Binary <- ComReceptor_Matrix
          ComReceptor_Binary[which(ComReceptor_Binary !=0)] <- 1

          index_ComLigand <- sapply(as.list(rem_ComL),function(x){which(rownames(GroupedUniqueLRMatrix) == x)})
          index_ComReceptor <- sapply(as.list(rem_ComR),function(x){which(rownames(GroupedUniqueLRMatrix) == x)})

          ComLigand_dd <- GroupedUniqueLRMatrix[index_ComLigand, which(colnames(GroupedUniqueLRMatrix) ==CellType_Index)]
          ComReceptor_dd <- GroupedUniqueLRMatrix[index_ComReceptor,]

          ComLigand_colnames <- CellType_Index
          ComReceptor_colnames <- sapply(strsplit(names(ComReceptor_Binary),"-"),"[[",1)


          Index_P <- sapply(as.list(rem_ComL),function(x){which(names(p) == x)}[1])
          Index_K <- sapply(as.list(rem_ComR),function(x){which(names(k) == x)}[1])

          rem_P <- p[Index_P]
          rem_K <- k[Index_K]

          ComLigand_dd[which(ComLigand_dd < rem_P)] <- 0
          ComLigand_dd <- ComLigand_dd/rem_P
          ComLigand_dd[which(ComLigand_dd > 7)] <- 7

          TargetLigand_dd_OUT <- ComLigand_dd
          ComReceptor_dd_Connection_OUT <- NULL


          for (i in 1:length(ComReceptor_colnames)){

            sub_ComReceptor_dd <- ComReceptor_dd[which(names(ComReceptor_dd)== i)]
            sub_ComReceptor_Binary <- ComReceptor_Binary[which(ComReceptor_colnames == ComReceptor_colnames[i])]
            tem_ComReceptor_dd <- sub_ComReceptor_dd*sub_ComReceptor_Binary
            tem_ComReceptor_dd[which(tem_ComReceptor_dd < rem_K)] <- 0
            tem_ComReceptor_dd <- tem_ComReceptor_dd/rem_P #use scaled expression level
            tem_ComReceptor_dd[which(tem_ComReceptor_dd > 7)] <- 7#if the scaled level is larger than 7, set as 7
            names(tem_ComReceptor_dd) <- paste0(names(tem_ComReceptor_dd),".",c(1:length(tem_ComReceptor_dd)))
            ComReceptor_dd_Connection_OUT <- c(ComReceptor_dd_Connection_OUT,tem_ComReceptor_dd)
          }


          names(TargetLigand_dd_OUT) <- paste0(names(TargetLigand_dd_OUT),".",c(1:length(TargetLigand_dd_OUT)))

          CellType_OUT <- c(TargetLigand_dd_OUT,ComReceptor_dd_Connection_OUT)
          names(CellType_OUT)[1:length(TargetLigand_dd_OUT)] <- paste0(names(TargetLigand_dd_OUT),"-L")
          names(CellType_OUT)[(length(TargetLigand_dd_OUT)+1):length(CellType_OUT)] <- paste0(names(CellType_OUT)[(length(TargetLigand_dd_OUT)+1):length(CellType_OUT)],"-R")


          rep_Num_OUT <- NULL
          rep_vec_OUT <- NULL
          for (i in 1:length(Unique_CellType)){
            rep_Num_OUT[i] <- length(which(names(ComReceptor_dd) == i))
            rep_subvec <- rep(i-1,rep_Num_OUT[i])
            rep_vec_OUT <- c(rep_vec_OUT,rep_subvec)

          }

          Target_rep <- rep(CellType_Index-1,length(ComLigand_dd))
          rep_vec_OUT <- c(Target_rep,rep_vec_OUT)


        } else if (!is.null(dim(TopMatrix_OUT))){

          rem_LRpairs <- rownames(TopMatrix_OUT)
          rem_ComL <- sapply(strsplit(rem_LRpairs,"--"),"[[",1)
          rem_ComR <- sapply(strsplit(rem_LRpairs,"--"),"[[",2)
          ComLigand_Matrix <- TopMatrix_OUT[,1]
          ComReceptor_Matrix <- TopMatrix_OUT[,2:ncol(TopMatrix_OUT)]

          ComLigand_Binary <- ComLigand_Matrix
          ComLigand_Binary[which(ComLigand_Binary != 0)] <- 1
          ComReceptor_Binary <- ComReceptor_Matrix
          ComReceptor_Binary[which(ComReceptor_Binary !=0)] <- 1

          index_ComLigand <- sapply(as.list(rem_ComL),function(x){which(rownames(GroupedUniqueLRMatrix) == x)})
          index_ComReceptor <- sapply(as.list(rem_ComR),function(x){which(rownames(GroupedUniqueLRMatrix) == x)})

          ComLigand_dd <- GroupedUniqueLRMatrix[index_ComLigand, which(colnames(GroupedUniqueLRMatrix) ==CellType_Index)]
          ComReceptor_dd <- GroupedUniqueLRMatrix[index_ComReceptor,]

          ComLigand_colnames <- CellType_Index
          if (is.null(ncol(ComReceptor_Binary))){
            ComReceptor_colnames <- sapply(strsplit(names(ComReceptor_Binary),"-"),"[[",1)
          }else {
            ComReceptor_colnames <- sapply(strsplit(colnames(ComReceptor_Binary),"-"),"[[",1)
          }

          Index_P <- sapply(as.list(rem_ComL),function(x){which(names(p) == x)}[1])
          Index_K <- sapply(as.list(rem_ComR),function(x){which(names(k) == x)}[1])

          rem_P <- p[Index_P]
          rem_K <- k[Index_K]

          ComLigand_dd[which(ComLigand_dd < rem_P)] <- 0
          ComLigand_dd <- ComLigand_dd/rem_P
          ComLigand_dd[which(ComLigand_dd > 7)] <- 7

          TargetLigand_dd_OUT <- ComLigand_dd
          ComReceptor_dd_Connection_OUT <- NULL


          for (i in 1:length(ComReceptor_colnames)){

            sub_ComReceptor_dd <- ComReceptor_dd[,which(colnames(ComReceptor_dd)== i)]
            sub_ComReceptor_Binary <- ComReceptor_Binary[,which(ComReceptor_colnames == ComReceptor_colnames[i])]
            tem_ComReceptor_dd <- sub_ComReceptor_dd*sub_ComReceptor_Binary
            tem_ComReceptor_dd[which(tem_ComReceptor_dd < rem_K)] <- 0
            tem_ComReceptor_dd <- tem_ComReceptor_dd/rem_P
            tem_ComReceptor_dd[which(tem_ComReceptor_dd > 7)] <- 7
            colnames(tem_ComReceptor_dd) <- paste0(colnames(tem_ComReceptor_dd),".",c(1:ncol(tem_ComReceptor_dd)))
            ComReceptor_dd_Connection_OUT <- cbind(ComReceptor_dd_Connection_OUT,tem_ComReceptor_dd)
          }


          colnames(TargetLigand_dd_OUT) <- paste0(colnames(TargetLigand_dd_OUT),".",c(1:ncol(TargetLigand_dd_OUT)))

          CellType_OUT <- cbind(TargetLigand_dd_OUT,ComReceptor_dd_Connection_OUT)
          rownames(CellType_OUT) <- paste0(rownames(TargetLigand_dd_OUT),"--",rownames(ComReceptor_dd_Connection_OUT))
          colnames(CellType_OUT)[1:ncol(TargetLigand_dd_OUT)] <- paste0(colnames(TargetLigand_dd_OUT),"-L")
          colnames(CellType_OUT)[(ncol(TargetLigand_dd_OUT)+1):ncol(CellType_OUT)] <- paste0(colnames(CellType_OUT)[(ncol(TargetLigand_dd_OUT)+1):ncol(CellType_OUT)],"-R")

          Unique_CellType <- sort(unique(colnames(ComReceptor_dd)))
          rep_Num_OUT <- NULL
          rep_vec_OUT <- NULL
          for (i in 1:length(Unique_CellType)){
            rep_Num_OUT[i] <- length(which(colnames(ComReceptor_dd) == i))
            rep_subvec <- rep(i-1,rep_Num_OUT[i])
            rep_vec_OUT <- c(rep_vec_OUT,rep_subvec)

          }
          Target_rep <- rep(CellType_Index-1,ncol(ComLigand_dd))
          rep_vec_OUT <- c(Target_rep,rep_vec_OUT)

        }
      }

    }
    if (Combine_InOut == TRUE){

      TopMatrix_InOut <- M_InOut
      if (sum(TopMatrix_InOut) == 0){
        CellType_InOut_Proportion <- 0
      }else if (sum(TopMatrix_InOut) != 0){

        if (is.null(dim(TopMatrix_InOut))){

          rem_LRpairs <- rownames(TopMatrix_InOut)
          rem_ComL <- sapply(strsplit(rem_LRpairs,"--"),"[[",1)
          rem_ComR <- sapply(strsplit(rem_LRpairs,"--"),"[[",2)
          Target_Ligand <- TopMatrix_InOut[,1]
          Receptor_OUT <- TopMatrix_InOut[,2:(ncol(TopMatrix_InOut)/2)]
          Ligand_IN <- TopMatrix_InOut[,((ncol(TopMatrix_InOut)/2)+1):(ncol(TopMatrix_InOut)-1)]
          Target_Receptor <- TopMatrix_InOut[,ncol(TopMatrix_InOut)]

          Ligand_Binary <- Target_Ligand
          Ligand_Binary[which(Ligand_Binary != 0)] <- 1
          Receptor_Binary <- Target_Receptor
          Receptor_Binary[which(Receptor_Binary != 0)] <- 1
          Receptor_OUT_Binary <- Receptor_OUT
          Receptor_OUT_Binary[which(Receptor_OUT_Binary != 0)] <- 1
          Ligand_IN_Binary <- Ligand_IN
          Ligand_IN_Binary[which(Ligand_IN_Binary != 0)] <- 1


          index_ComLigand <- sapply(as.list(rem_ComL),function(x){which(rownames(GroupedUniqueLRMatrix) == x)})
          index_ComReceptor <- sapply(as.list(rem_ComR),function(x){which(rownames(GroupedUniqueLRMatrix) == x)})

          TargetLigand_dd <- GroupedUniqueLRMatrix[index_ComLigand, which(colnames(GroupedUniqueLRMatrix) ==CellType_Index)]
          TargetReceptor_dd <- GroupedUniqueLRMatrix[index_ComReceptor,which(colnames(GroupedUniqueLRMatrix) ==CellType_Index)]
          Ligand_IN_dd <- GroupedUniqueLRMatrix[index_ComLigand,]
          Receptor_OUT_dd <- GroupedUniqueLRMatrix[index_ComReceptor,]

          Index_P <- sapply(as.list(rem_ComL),function(x){which(names(p) == x)}[1])
          Index_K <- sapply(as.list(rem_ComR),function(x){which(names(k) == x)}[1])

          rem_P <- p[Index_P]
          rem_K <- k[Index_K]

          TargetLigand_dd <- TargetLigand_dd*Ligand_Binary
          TargetReceptor_dd <- TargetReceptor_dd*Receptor_Binary
          TargetLigand_dd[which(TargetLigand_dd < rem_P)] <- 0
          TargetLigand_dd <- TargetLigand_dd/rem_P
          TargetLigand_dd[which(TargetLigand_dd > 7)] <- 7#if bigger than 7 set as 7
          TargetReceptor_dd[which(TargetReceptor_dd < rem_K)] <- 0
          TargetReceptor_dd <- TargetReceptor_dd/rem_K
          TargetReceptor_dd[which(TargetReceptor_dd > 7)] <- 7

          Ligand_colnames <- sapply(strsplit(names(Ligand_IN_Binary),"-"),"[[",1)
          Receptor_colnames <- sapply(strsplit(names(Receptor_OUT_Binary),"-"),"[[",1)
          Unique_CellType <- unique(names(Ligand_IN_dd))


          Receptor_OUT_Connection <- NULL
          Ligand_IN_Connection <- NULL


          for (i in 1:length(Unique_CellType)){

            sub_Ligand_dd <- Ligand_IN_dd[which(names(Ligand_IN_dd)== i)]
            sub_Ligand_Binary <- Ligand_IN_Binary[which(Ligand_colnames == Ligand_colnames[i])]
            tem_Ligand_dd <- sub_Ligand_dd*sub_Ligand_Binary
            tem_Ligand_dd[which(tem_Ligand_dd < rem_P)] <- 0
            tem_Ligand_dd <- tem_Ligand_dd/rem_P
            tem_Ligand_dd[which(tem_Ligand_dd > 7)] <- 7
            names(tem_Ligand_dd) <- paste0(names(tem_Ligand_dd),".",c(1:length(tem_Ligand_dd)))
            Ligand_IN_Connection <- c(Ligand_IN_Connection,tem_Ligand_dd)

            sub_ComReceptor_dd <- Receptor_OUT_dd[which(names(Receptor_OUT_dd)== i)]
            sub_ComReceptor_Binary <- Receptor_OUT_Binary[which(Receptor_colnames == Receptor_colnames[i])]
            tem_ComReceptor_dd <- sub_ComReceptor_dd*sub_ComReceptor_Binary
            tem_ComReceptor_dd[which(tem_ComReceptor_dd < rem_K)] <- 0
            tem_ComReceptor_dd <- tem_ComReceptor_dd/rem_K
            tem_ComReceptor_dd[which(tem_ComReceptor_dd > 7)] <- 7
            names(tem_ComReceptor_dd) <- paste0(names(tem_ComReceptor_dd),".",c(1:length(tem_ComReceptor_dd)))
            Receptor_OUT_Connection <- c(Receptor_OUT_Connection,tem_ComReceptor_dd)

          }


          names(TargetLigand_dd) <- paste0(names(TargetLigand_dd),".",c(1:length(TargetLigand_dd)))
          names(TargetReceptor_dd) <- paste0(names(TargetReceptor_dd),".",c(1:length(TargetReceptor_dd)))


          CellType_InOut_Proportion <- c(TargetLigand_dd,Receptor_OUT_Connection,Ligand_IN_Connection,TargetReceptor_dd)
          names(CellType_InOut_Proportion)[1:length(TargetLigand_dd)] <- paste0(names(TargetLigand_dd),"_L")
          names(CellType_InOut_Proportion)[(length(TargetLigand_dd)+1):(length(CellType_InOut_Proportion)/2)] <- paste0(names(CellType_InOut_Proportion)[(length(TargetLigand_dd)+1):(length(CellType_InOut_Proportion)/2)],"-R")
          names(CellType_InOut_Proportion)[((length(CellType_InOut_Proportion)/2)+1):(length(CellType_InOut_Proportion)-length(TargetReceptor_dd))] <- paste0(names(CellType_InOut_Proportion)[((length(CellType_InOut_Proportion)/2)+1):(length(CellType_InOut_Proportion)-length(TargetReceptor_dd))],"-L")
          names(CellType_InOut_Proportion)[((length(CellType_InOut_Proportion)-length(TargetReceptor_dd))+1):(length(CellType_InOut_Proportion))] <- paste0(names(CellType_InOut_Proportion)[((length(CellType_InOut_Proportion)-length(TargetReceptor_dd))+1):(length(CellType_InOut_Proportion))],"_R")

          rep_Num_INOUT <- NULL
          rep_vec_INOUT <- NULL
          for (i in 1:length(Unique_CellType)){
            rep_Num_INOUT[i] <- length(which(names(Ligand_IN_dd) == i))
            rep_subvec <- rep(i-1,rep_Num_INOUT[i])
            rep_vec_INOUT <- c(rep_vec_INOUT,rep_subvec)

          }

          Target_rep <- rep(CellType_Index-1,length(TargetLigand_dd))
          rep_vec_INOUT <- c(Target_rep,rep_vec_INOUT,rep_vec_INOUT,Target_rep)

        }else if (!is.null(dim(TopMatrix_InOut))){
          rem_LRpairs <- rownames(TopMatrix_InOut)
          rem_ComL <- sapply(strsplit(rem_LRpairs,"--"),"[[",1)
          rem_ComR <- sapply(strsplit(rem_LRpairs,"--"),"[[",2)

          Target_Ligand <- TopMatrix_InOut[,1]
          Receptor_OUT <- TopMatrix_InOut[,2:(ncol(TopMatrix_InOut)/2)]
          Ligand_IN <- TopMatrix_InOut[,((ncol(TopMatrix_InOut)/2)+1):(ncol(TopMatrix_InOut)-1)]
          Target_Receptor <- TopMatrix_InOut[,ncol(TopMatrix_InOut)]

          Ligand_Binary <- Target_Ligand
          Ligand_Binary[which(Ligand_Binary != 0)] <- 1
          Receptor_Binary <- Target_Receptor
          Receptor_Binary[which(Receptor_Binary != 0)] <- 1
          Receptor_OUT_Binary <- Receptor_OUT
          Receptor_OUT_Binary[which(Receptor_OUT_Binary != 0)] <- 1
          Ligand_IN_Binary <- Ligand_IN
          Ligand_IN_Binary[which(Ligand_IN_Binary != 0)] <- 1


          index_ComLigand <- sapply(as.list(rem_ComL),function(x){which(rownames(GroupedUniqueLRMatrix) == x)})
          index_ComReceptor <- sapply(as.list(rem_ComR),function(x){which(rownames(GroupedUniqueLRMatrix) == x)})

          TargetLigand_dd <- GroupedUniqueLRMatrix[index_ComLigand, which(colnames(GroupedUniqueLRMatrix) ==CellType_Index)]
          TargetReceptor_dd <- GroupedUniqueLRMatrix[index_ComReceptor,which(colnames(GroupedUniqueLRMatrix) ==CellType_Index)]
          Ligand_IN_dd <- GroupedUniqueLRMatrix[index_ComLigand,]
          Receptor_OUT_dd <- GroupedUniqueLRMatrix[index_ComReceptor,]

          Index_P <- sapply(as.list(rem_ComL),function(x){which(names(p) == x)}[1])
          Index_K <- sapply(as.list(rem_ComR),function(x){which(names(k) == x)}[1])

          rem_P <- p[Index_P]
          rem_K <- k[Index_K]

          TargetLigand_dd <- TargetLigand_dd*Ligand_Binary
          TargetReceptor_dd <- TargetReceptor_dd*Receptor_Binary
          TargetLigand_dd[which(TargetLigand_dd < rem_P)] <- 0
          TargetLigand_dd <- TargetLigand_dd/rem_P
          TargetLigand_dd[which(TargetLigand_dd > 7)] <- 7
          TargetReceptor_dd[which(TargetReceptor_dd < rem_K)] <- 0
          TargetReceptor_dd <- TargetReceptor_dd/rem_K
          TargetReceptor_dd[which(TargetReceptor_dd > 7)] <- 7

          Ligand_colnames <- sapply(strsplit(colnames(Ligand_IN_Binary),"-"),"[[",1)
          Receptor_colnames <- sapply(strsplit(colnames(Receptor_OUT_Binary),"-"),"[[",1)
          Unique_CellType <- unique(colnames(Ligand_IN_dd))


          Receptor_OUT_Connection <- NULL
          Ligand_IN_Connection <- NULL


          for (i in 1:length(Unique_CellType)){

            sub_Ligand_dd <- Ligand_IN_dd[,which(colnames(Ligand_IN_dd) == i)]
            sub_Ligand_Binary <- Ligand_IN_Binary[,which(Ligand_colnames == Ligand_colnames[i])]
            tem_Ligand_dd <- sub_Ligand_dd*sub_Ligand_Binary
            tem_Ligand_dd[which(tem_Ligand_dd < rem_P)] <- 0
            tem_Ligand_dd <- tem_Ligand_dd/rem_P
            tem_Ligand_dd[which(tem_Ligand_dd > 7)] <- 7
            colnames(tem_Ligand_dd) <- paste0(colnames(tem_Ligand_dd),".",c(1:ncol(tem_Ligand_dd)))
            Ligand_IN_Connection <- cbind(Ligand_IN_Connection,tem_Ligand_dd)

            sub_ComReceptor_dd <- Receptor_OUT_dd[,which(colnames(Receptor_OUT_dd)== i)]
            sub_ComReceptor_Binary <- Receptor_OUT_Binary[,which(Receptor_colnames == Receptor_colnames[i])]
            tem_ComReceptor_dd <- sub_ComReceptor_dd*sub_ComReceptor_Binary
            tem_ComReceptor_dd[which(tem_ComReceptor_dd < rem_K)] <- 0
            tem_ComReceptor_dd <- tem_ComReceptor_dd/rem_K
            tem_ComReceptor_dd[which(tem_ComReceptor_dd > 7)] <- 7
            colnames(tem_ComReceptor_dd) <- paste0(colnames(tem_ComReceptor_dd),".",c(1:ncol(tem_ComReceptor_dd)))
            Receptor_OUT_Connection <- cbind(Receptor_OUT_Connection,tem_ComReceptor_dd)

          }


          colnames(TargetLigand_dd) <- paste0(colnames(TargetLigand_dd),".",c(1:ncol(TargetLigand_dd)))
          colnames(TargetReceptor_dd) <- paste0(colnames(TargetReceptor_dd),".",c(1:ncol(TargetReceptor_dd)))


          CellType_InOut_Proportion <- cbind(TargetLigand_dd,Receptor_OUT_Connection,Ligand_IN_Connection,TargetReceptor_dd)
          rownames(CellType_InOut_Proportion) <- paste0(rownames(Ligand_IN_Connection),"--",rownames(Receptor_OUT_Connection))

          colnames(CellType_InOut_Proportion)[1:ncol(TargetLigand_dd)] <- paste0(colnames(TargetLigand_dd),"_L")
          colnames(CellType_InOut_Proportion)[(ncol(TargetLigand_dd)+1):(ncol(CellType_InOut_Proportion)/2)] <- paste0(colnames(CellType_InOut_Proportion)[(ncol(TargetLigand_dd)+1):(ncol(CellType_InOut_Proportion)/2)],"-R")
          colnames(CellType_InOut_Proportion)[((ncol(CellType_InOut_Proportion)/2)+1):(ncol(CellType_InOut_Proportion)-ncol(TargetReceptor_dd))] <- paste0(colnames(CellType_InOut_Proportion)[((ncol(CellType_InOut_Proportion)/2)+1):(ncol(CellType_InOut_Proportion)-ncol(TargetReceptor_dd))],"-L")

          colnames(CellType_InOut_Proportion)[((ncol(CellType_InOut_Proportion)-ncol(TargetReceptor_dd))+1):(ncol(CellType_InOut_Proportion))] <- paste0(colnames(CellType_InOut_Proportion)[((ncol(CellType_InOut_Proportion)-ncol(TargetReceptor_dd))+1):(ncol(CellType_InOut_Proportion))],"_R")

          rep_Num_INOUT <- NULL
          rep_vec_INOUT <- NULL
          for (i in 1:length(Unique_CellType)){
            rep_Num_INOUT[i] <- length(which(colnames(Ligand_IN_dd) == i))
            rep_subvec <- rep(i-1,rep_Num_INOUT[i])
            rep_vec_INOUT <- c(rep_vec_INOUT,rep_subvec)

          }

          Target_rep <- rep(CellType_Index-1,ncol(TargetLigand_dd))
          rep_vec_INOUT <- c(Target_rep,rep_vec_INOUT,rep_vec_INOUT,Target_rep)#for the annotation of pheatmap

        }
      }
    }
    Unique_CellType <- colnames(LRC3_list$contact_number_Remaining)
  }

  ###########################################

  if (TargetCellType_IN == TRUE) {

    if (sum(M_IN) == 0){
      print("There is no IN signal to this cell type!")
    } else if (sum(M_IN)!= 0){
        if (nrow(M_IN) == 1){
        print("There is only 1 Ligand-Receptor pair used in this cell type and it is neglected! Please check from the file of total List Of LRpairs!")
      }
      else if(nrow(M_IN) != 1){
        if (SaveImage == TRUE) {
          pdf(file="InOut_Seperate%03d.pdf")
          # png(file="Rplot%03d.png")
        }

        if (By_Proportion == FALSE){

          if (SaveImage == FALSE){
            #x11()
          }

          annotation_IN <- data.frame(LR = factor(c(rep(0,1),
                                                    rep(1,ncol(M_IN )-1)),
                                                  labels = c("Receptor", "Ligand")))
          rownames(annotation_IN) <- colnames(M_IN)

          LR <- c("navyblue", "red")
          names(LR) <- c("Receptor", "Ligand")
          anno_colors <- list(LR = LR)
          unique_labels <- sort(unique(Index_IN))

          annotationrow_IN <- data.frame(N = factor(c(Index_IN),
                                                    labels = unique_labels))
          rownames(annotationrow_IN) <- rownames(M_IN )


          if (ViewTop_LRpairs == TRUE){
            if (nrow(M_IN) < TOP_Number){
              M_IN <- M_IN[1:nrow(M_IN),]
            } else if (nrow(M_IN) >= TOP_Number)
            M_IN <- M_IN[1:TOP_Number,]
          }

          pheatmap::pheatmap(M_IN,
                             cluster_rows = F, cluster_cols = F,
                             # colorRampPalette(c("black","green","red"))(100),
                             fontsize=7,fontsize_col=7, fontsize_row = 7,main = "IN",
                             gaps_col = 1,
                             annotation_col = annotation_IN,
                             annotation_row = annotationrow_IN,
                             annotation_colors = anno_colors)



        } else if (By_Proportion == TRUE){

          annotation <- data.frame(CellType = factor(rep_vec_IN,labels = c(Unique_CellType)),LR = factor(c(rep(0,ncol(TargetReceptor_dd_IN)),
                                                                                                           rep(1,ncol(ComLigand_dd_Connection_IN))),
                                                                                                         labels = c("Receptor", "Ligand")))
          rownames(annotation) <- colnames(CellType_IN)

          LR <- c("red", "navyblue")

          names(LR) <- c("Receptor","Ligand")
          anno_colors <- list(LR = LR)

          unique_labels <- sort(unique(Index_IN))
          annotationrow_IN <- data.frame(N = factor(c(Index_IN),
                                                    labels = unique_labels))
          rownames(annotationrow_IN) <- rownames(M_IN)


          if (SaveImage == FALSE){
            #x11()
          }

          if (ViewTop_LRpairs == TRUE){
            if (nrow(CellType_IN) < TOP_Number){
              # TOP_Number <- nrow(CellType_IN)
              CellType_IN <- CellType_IN[1:nrow(CellType_IN),]
            } else if (nrow(CellType_IN) >= TOP_Number){
              CellType_IN <- CellType_IN[1:TOP_Number,]
              # CellType_OUT <- CellType_OUT[1:TOP_Number,]
            }
          }

          pheatmap::pheatmap(CellType_IN,
                             cluster_rows = F, cluster_cols = F,
                             color = colorRampPalette(c("#ffffff","#000000","#007700","#00bb00","#44cc00","#77aa00","#997700","#cc5500","#ff0000"))(1000),
                             fontsize=7,fontsize_col=7, fontsize_row = 7,main = "IN",
                             gaps_col =  gap_In,
                             show_colnames = F,
                             annotation_col = annotation,
                             annotation_row = annotationrow_IN,
                             annotation_colors = anno_colors)

        }
      }
      }
      }

    if (TargetCellType_OUT == TRUE){
      if (sum(M_OUT) == 0){
        print("There is no OUT signal to this cell type!")
      } else if (sum(M_OUT)!= 0){
        # if (is.null(dim(M_OUT))){
        if (nrow(M_OUT) == 1){
          print("There is only 1 Ligand-Receptor pair used in this cell type and it is neglected! Please check from the file of total List Of LRpairs!")
        }
        else if(nrow(M_OUT) != 1){

          if (SaveImage == TRUE) {
            pdf(file="InOut_Seperate%03d.pdf")
            # png(file="Rplot%03d.png")
          }

          if (By_Proportion == FALSE){

            # if (SaveImage == FALSE){
            #   #x11()
            # }

            ################

            annotation_OUT <- data.frame(LR = factor(c(rep(0,1),
                                                       rep(1,ncol(M_OUT)-1)),
                                                     labels = c("Ligand", "Receptor")))
            rownames(annotation_OUT) <- colnames(M_OUT)

            LR <- c("red", "navyblue")
            names(LR) <- c("Ligand", "Receptor")
            anno_colors <- list(LR = LR)
            unique_labels <- sort(unique(Index_OUT))

            annotationrow_OUT <- data.frame(N = factor(c(Index_OUT),
                                                       labels = unique_labels))
            rownames(annotationrow_OUT) <- rownames(M_OUT)

            if (ViewTop_LRpairs == TRUE){
              if (nrow(M_OUT) < TOP_Number){
                # TOP_Number <- nrow(M_OUT)
                M_OUT <- M_OUT[1:nrow(M_OUT),]
              } else if (nrow(M_OUT) >= TOP_Number){
                # M_IN <- M_IN[1:TOP_Number,]
                M_OUT <- M_OUT[1:TOP_Number,]
              }
            }


            ######## plot the OUT

            if (SaveImage == FALSE){
              #x11()
            }
            pheatmap::pheatmap(M_OUT,
                               cluster_rows = F, cluster_cols = F,
                               fontsize=7,fontsize_col=7, fontsize_row = 7,
                               main = "OUT",
                               gaps_col =  1,
                               annotation_col = annotation_OUT,
                               annotation_row = annotationrow_OUT,
                               annotation_colors = anno_colors)

          } else if (By_Proportion == TRUE){

            annotation <- data.frame(CellType = factor(rep_vec_OUT,labels = c(Unique_CellType)),LR = factor(c(rep(0,ncol(TargetLigand_dd_OUT)),
                                                                                                              rep(1,ncol(ComReceptor_dd_Connection_OUT))),
                                                                                                            labels = c("Ligand", "Receptor")))
            rownames(annotation) <- colnames(CellType_OUT)

            LR <- c("navyblue","red")

            names(LR) <- c("Ligand", "Receptor")
            anno_colors <- list(LR = LR)
            unique_labels <- sort(unique(Index_OUT))

            annotationrow_OUT <- data.frame(N = factor(c(Index_OUT),
                                                       labels = unique_labels))
            rownames(annotationrow_OUT) <- rownames(CellType_OUT)

            if (ViewTop_LRpairs == TRUE){
              if (nrow(CellType_OUT) < TOP_Number){
                # TOP_Number <- nrow(CellType_OUT)
                CellType_OUT <- CellType_OUT[1:nrow(CellType_OUT),]

              } else if (nrow(CellType_OUT) >= TOP_Number){
                # CellType_IN <- CellType_IN[1:TOP_Number,]
                CellType_OUT <- CellType_OUT[1:TOP_Number,]
              }

            }

            if (SaveImage == FALSE){
              #x11()
            }
            pheatmap::pheatmap(CellType_OUT,
                               cluster_rows = F, cluster_cols = F,
                               main = "OUT",
                               color = colorRampPalette(c("#ffffff","#000000","#007700","#00bb00","#44cc00","#77aa00","#997700","#cc5500","#ff0000"))(1000),
                               fontsize=7,fontsize_col=7, fontsize_row = 7,
                               gaps_col =  ncol(ComLigand_dd),
                               show_colnames = F,
                               annotation_col = annotation,
                               annotation_row = annotationrow_OUT,
                               annotation_colors = anno_colors)
          }
        }

        }

    }


   if (Combine_InOut == TRUE){

    if (sum(M_InOut) == 0){
      print("There is no IN or OUT signal in this cell type!")
    } else if (sum(M_InOut)!= 0){
      # if (is.null(dim(M_InOut))){
      if (nrow(M_InOut) ==1){
        print("There is only 1 Ligand-Receptor pair used in this cell type and it is neglected! Please check from the file of total List Of LRpairs!")
      } else if(nrow(M_InOut) !=1){
        if (SaveImage == TRUE) {
          pdf(file="InOut_Combine%03d.pdf")
          # png(file="Rplot%03d.png")
        }
        if (SaveImage == FALSE){
          #x11()
        }
        if (By_Proportion == FALSE){
          annotation <- data.frame(LR = factor(c(rep(0,1),
                                                 rep(1,ncol(M_InOut )/2-1),rep(0,ncol(M_InOut )/2-1),1),
                                               labels = c("Ligand", "Receptor")))
          rownames(annotation) <- colnames(M_InOut)
          LR <- c("red", "navyblue")
          names(LR) <- c("Receptor", "Ligand")
          anno_colors <- list(LR = LR)
          unique_labels <- sort(unique(Index_InOut))

          annotationrow <- data.frame(N = factor(c(Index_InOut),
                                                 labels = unique_labels))
          rownames(annotationrow) <- rownames(M_InOut)

          if (ViewTop_LRpairs == TRUE){
            if (nrow(M_InOut) < TOP_Number){
              # TOP_Number <- nrow(M_InOut)
              M_InOut <- M_InOut[1:nrow(M_InOut),]
            } else if (nrow(M_InOut) >= TOP_Number) {
              M_InOut <- M_InOut[1:TOP_Number,]
              # CellType_InOut_Proportion <- CellType_InOut_Proportion[1:TOP_Number,]
            }
          }
          pheatmap::pheatmap(M_InOut,
                             cluster_rows = F, cluster_cols = F,
                             # colorRampPalette(c("black","green","red"))(100),
                             fontsize=7,fontsize_col=7, fontsize_row = 7,
                             gaps_col =  c(1,ncol(M_InOut)/2,ncol(M_InOut)-1),
                             annotation_col = annotation,
                             main = "InOut",
                             annotation_row = annotationrow,
                             annotation_colors = anno_colors)
        }
        else if (By_Proportion == TRUE){
          annotation <- data.frame(CellType = factor(rep_vec_INOUT,labels = c(Unique_CellType)),LR = factor(c(rep(0,ncol(TargetLigand_dd)),
                                                                                                              rep(1,ncol(Receptor_OUT_Connection)),rep(0,ncol(Ligand_IN_Connection)),rep(1,ncol(TargetReceptor_dd))),
                                                                                                            labels = c("Ligand", "Receptor")))
          rownames(annotation) <- colnames(CellType_InOut_Proportion)
          LR <- c("navyblue","red")
          names(LR) <- c("Ligand", "Receptor")
          anno_colors <- list(LR = LR)
          unique_labels <- sort(unique(Index_InOut))
          annotationrow <- data.frame(N = factor(c(Index_InOut),
                                                 labels = unique_labels))
          rownames(annotationrow) <- rownames(M_InOut)
          if (ViewTop_LRpairs == TRUE){
            if (nrow(CellType_InOut_Proportion) < TOP_Number){
              # TOP_Number <- nrow(CellType_InOut_Proportion)
              CellType_InOut_Proportion <- CellType_InOut_Proportion[1:nrow(CellType_InOut_Proportion),]

            } else if (nrow(CellType_InOut_Proportion) >= TOP_Number){
              # M_InOut <- M_InOut[1:TOP_Number,]
              CellType_InOut_Proportion <- CellType_InOut_Proportion[1:TOP_Number,]

            }
          }
          pheatmap::pheatmap(CellType_InOut_Proportion,
                             cluster_rows = F, cluster_cols = F,
                             color = colorRampPalette(c("#ffffff","#000000","#007700","#00bb00","#44cc00","#77aa00","#997700","#cc5500","#ff0000"))(1000),
                             fontsize=7,fontsize_col=7, fontsize_row = 7,
                             gaps_col = c(ncol(TargetLigand_dd),ncol(CellType_InOut_Proportion)/2,(ncol(CellType_InOut_Proportion)-ncol(TargetReceptor_dd))),
                             show_colnames = F,
                             main = "InOut",
                             annotation_col = annotation,
                             # annotation_colors = anno_colors,
                             annotation_row = annotationrow,
                             annotation_colors = anno_colors)
        }

      }
    }
  }
  if (SaveImage == TRUE){
    dev.off ()
  }
}

#' Search the ligand-receptor pathways between two cell types.
#'
#' This function searches the ligand and receptor pairs that are used between two cell types.
#' It will give all pairs that being used between two cell types in one direction: CellType1 express ligand, CellType2 express receptor.
#' The ligand-receptor pairs are ranked by the signal of the pathway,
#' which is the mean of the two scaled expression value of both ligand and receptor in the corresponding pairwise cell types.
#' The scaled expression value is gained through the original mean expression value of a gene in a cell type devied by the threshold of the gene.
#' @param LRC3_list the input data generated by function LRC3_INF.
#' @param CellType1 the Index of one cell type in query. CellType1 is the cell type that significantly express ligand.
#' @param CellType2 the Index of one cell type in query. CellType2 is the cell type that significantly express receptor.
#' @export

LRC3_LRPairsCellTypeSearch <- function(LRC3_list,CellType1,CellType2){
  Remaining_scaled_contactINF <- LRC3_list$Remaining_scaled_contactINF
  rem_LRpairs <- LRC3_list$rem_LRpairs

  Scaled_LRPairs_Vec <- Remaining_scaled_contactINF[CellType1,CellType2,]
  InUse_Positions <- which(Scaled_LRPairs_Vec > 0)
  Scaled_Signal <- Scaled_LRPairs_Vec[InUse_Positions[order(Scaled_LRPairs_Vec[InUse_Positions],decreasing = T)]]

  Scaled_LRPairs_InUse <-rem_LRpairs[InUse_Positions[order(Scaled_LRPairs_Vec[InUse_Positions],decreasing = T)],]
  Scaled_LRPairs_InUse[,3] <- round(Scaled_Signal,1)
  if (length(InUse_Positions) != 0){
    rownames(Scaled_LRPairs_InUse) <- c(1:nrow(Scaled_LRPairs_InUse))
    colnames(Scaled_LRPairs_InUse) <- c("Ligand","Receptor","Signal")
  }
  Scaled_LRPairs_InUse_Transpose <- t(Scaled_LRPairs_InUse)
  # Scaled_LRPairs_InUse_Transpose <- t(Scaled_LRPairs_InUse[1:300,])
  Table_Scaled_LRpairs_InUse <- as.table(Scaled_LRPairs_InUse_Transpose)
  return(Table_Scaled_LRpairs_InUse)
}


#' Total list of ligand-receptor pairs.
#'
#' This function gives a table that contains the information of ligand-receptor pairs used between pairwised cell types.
#' It includes the name of the ligands and receptors, the signal of the pathways,
#' the original mean expression value of a gene, the threshold for a gene, and cell type names that uses the pathways.
#' The pairs are ranked by their signal of the pathways, which is the mean of the two expression value of ligand and receptor in the pairwised cell types.
#'
#' @param LRC3_list the input data generated by function LRC3_INF.
#' @return ListOfLRpairs a table contains information of ligand-receptor pairs used between pairwised cell types.
#' In the table, the first two colomn are the names of ligand and receptor pairs.
#' Signal is the mean of the two expression value of ligand and receptor in the pairwised cell types.
#' L_Signal is the signal for the ligand, which is the factor of mean expression value of the ligand in a cell type devided by the threshold of the ligand.
#' R_Signal is the signal for the receptor, which is the factor of the mean expression value of the receptor in a cell type devided by the threshold of the receptor.
#' L_Origin is the mean expression value of the ligand in a cell type.
#' R_Origin is the mean expression value of the receptor in a cell type.
#' L_Thresh is the threshold of the ligand.
#' R_Thresh is the threshold of the receptor.
#' Celltypes are the two interacting cell types that are using the pathway.
#'
#' @export
LRC3_TotalListOfLRpairs <- function(LRC3_list){
  Remaining_scaled_contactINF <- LRC3_list$Remaining_scaled_contactINF
  rem_LRpairs <- LRC3_list$rem_LRpairs
  p <- LRC3_list$p
  k <- LRC3_list$k
  ScaledLigandExpressionLevelINF <- LRC3_list$ScaledLigandExpressionLevelINF
  ScaledReceptorExpressionLevelINF <- LRC3_list$ScaledReceptorExpressionLevelINF
  LigandExpressionLevelINF <- LRC3_list$LigandExpressionLevelINF
  ReceptorExpressionLevelINF <- LRC3_list$ReceptorExpressionLevelINF

  ListOfLRpairs <- NULL
  for (i in 1:ncol(Remaining_scaled_contactINF[,,1])){
    for (j in 1:ncol(Remaining_scaled_contactINF[,,1])){
      Scaled_LRPairs_Vec <- Remaining_scaled_contactINF[i,j,]
      InUse_Positions <- which(Scaled_LRPairs_Vec > 0)
      orderInusePosition <- InUse_Positions[order(Scaled_LRPairs_Vec[InUse_Positions],decreasing = T)]
      Scaled_Signal <- Scaled_LRPairs_Vec[orderInusePosition]

      Scaled_LigandSignal <- ScaledLigandExpressionLevelINF[i,j,][orderInusePosition]
      Scaled_ReceptorSignal <- ScaledReceptorExpressionLevelINF[i,j,][orderInusePosition]

      Pval_Ligand <- p[orderInusePosition]
      Pval_Receptor <- k[orderInusePosition]

      Ligand_ExpressionLevel <- LRC3_list$LigandExpressionLevelINF[i,j,][orderInusePosition]
      Receptor_ExpressionLevel <- LRC3_list$ReceptorExpressionLevelINF[i,j,][orderInusePosition]

      Scaled_LRPairs_InUse <-rem_LRpairs[InUse_Positions[order(Scaled_LRPairs_Vec[InUse_Positions],decreasing = T)],]
      Scaled_LRPairs_InUse[,3] <- round(Scaled_Signal,2)
      Scaled_LRPairs_InUse[,4] <- round(Scaled_LigandSignal,2)
      Scaled_LRPairs_InUse[,5] <- round(Scaled_ReceptorSignal,2)
      Scaled_LRPairs_InUse[,6] <- round(Ligand_ExpressionLevel,4)
      Scaled_LRPairs_InUse[,7] <- round(Receptor_ExpressionLevel,4)
      Scaled_LRPairs_InUse[,8] <- round(Pval_Ligand,4)
      Scaled_LRPairs_InUse[,9] <- round(Pval_Receptor,4)

      celltype_names <- colnames(LRC3_list$contact_number_Remaining)
      if (length(Scaled_Signal) != 0){
        Scaled_LRPairs_InUse[,10] <- paste0(celltype_names[i],"->",celltype_names[j])
      }
      if (length(InUse_Positions) != 0){
        rownames(Scaled_LRPairs_InUse) <- c(1:nrow(Scaled_LRPairs_InUse))
        colnames(Scaled_LRPairs_InUse) <- c("Ligand","Receptor","Signal","L_Signal","R_Signal","L_Origin","R_Origin","L_Thresh","R_Thresh","Celltypes")
      }

      ListOfLRpairs <- rbind(ListOfLRpairs,Scaled_LRPairs_InUse)
    }
  }
  # write.csv(ListOfLRpairs,file = paste0("/Users/xs2/R/Liver_Gurdon/",paste(NameOfData,".csv")))
  return(ListOfLRpairs)
}