R/ligrec.R

In Morriseylab/ligrec:

ligrec <- function(scrna,grp.var,org,perc){
  #get grouping variable
  var=as.character(grp.var)
  tt=rownames(GetAssayData(object = scrna, slot = "counts"))

  #genes=fread("data/ligrecgenes.txt",header = TRUE)
  if(org=="mouse"){
    data('Mm_PairsLigRec')
    rl = mm
  }else if(org=="human"){
    data('Hs_PairsLigRec')
    rl=hs
    }
  genes=unique(c(as.character(rl$ligand),as.character(rl$receptor)))
  genes2=tt[tt %in% genes]
  #For all unique genes in the ligrec list, get their expression value for all cells and the groups the cells belong to
  my.data=FetchData(scrna,c(var,genes2))
  colnames(my.data)[1]= "clust"
  perc=perc/100
  result=data.frame()
  res=data.frame()
  #loop over each cluster to find pairs
  for(i in 1:(length(levels(my.data$clust)))){
    for(j in 1:(length(levels(my.data$clust)))){
      #from the large martix, subselect receptor and lig subgoups (if i=1 and j=2, keep cells in grps 1 and 2)
      test=my.data[my.data$clust==levels(my.data$clust)[i] | my.data$clust==levels(my.data$clust)[j],]
      #Subselect genes in receptor list in cells in rec subgroup (say 1)
      R_c1=test[test$clust==levels(my.data$clust)[i] ,(colnames(test) %in% rl$receptor)]
      #Subselect genes in ligand list in cells in lig subgroup (say 2)
      L_c2=test[test$clust==levels(my.data$clust)[j] , (colnames(test) %in% rl$ligand)]
      if(nrow(R_c1)!=0 &nrow(L_c2)!=0){
        #keep genes that are expressed in more than user-input percent of the cells
        keep1 = colSums(R_c1>0)>=perc*dim(R_c1)[1]
        keep2 = colSums(L_c2>0)>=perc*dim(L_c2)[1]
        R_c1=R_c1[,keep1]
        L_c2=L_c2[,keep2]
        #get list of lig-rec pairs
        res=rl[(rl$ligand %in% colnames(L_c2)) & (rl$receptor %in% colnames(R_c1)),]
      }else{}
      if(nrow(res)!=0){
        res$Receptor_cluster=levels(my.data$clust)[i]
        res$Lig_cluster=levels(my.data$clust)[j]
        result=rbind(result,res)
      }else{result=result}
    }
  }
  # get final list of all lig-rec pairs
  #result=result[result$Receptor_cluster!=result$Lig_cluster,]
  return(result)
}