LRinteractions/1_Multiple_Imputation.R at master · TransBioInfoLab/LRinteractions · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
###################################################################
### <June 20 2022> Min Lu and Yifan Sha
###
###
###
###################################################################
###
### An R program for identifying ligand-receptor gene interactions for scRNA-seq
###  using Random Forest based methods
###
### R code for "LR hunting: a Random Forests based cell-cell interaction
###             discovery method for single-cell gene expression data"
### ------------------------------------------------------------
###  Min Lu, PhD                     m.lu6@umiami.edu
###                                  phone: 305-243-5473
###  Research Assistant Professor, Division of Biostatistics
###  Clinical Research Building
###  1120 NW 14th Street, Room 1059
###  University of Miami, Miami FL 33136
###
###  -------------------------------------------------------------
###  THIS PROGRAM SHOULD NOT BE COPIED, USED, MODIFIED, OR
###  DISSEMINATED IN ANY WAY WITHOUT SPECIFIC WRITTEN PERMISSION
###  FROM THE AUTHOR.
###
###################################################################
library(randomForestSRC)

load("LogNormalized_9types.data")
nam <- rownames(LogNormalized.data)

Lig.cell <- "Myeloid"
Recp.cell <-"Epithelial_Basal"
pair <- 1
patient <- 1
setwd(paste("/Volumes/Extreme SSD/all/folders/pair", pair, "patient", patient, sep = ""))

Lig.id <- which((info$patientID == "P1") & (info$celltype_final == Lig.cell) )
Recp.id <- which((info$patientID == "P1") & (info$celltype_final == Recp.cell) )
setwd("/Volumes/HDD/pair7patient1")

Lig.n <- length(Lig.id)
Recp.n <- length(Recp.id)

lig <- unique(LRdb$ligand)
recp <- unique(LRdb$receptor)

varsel.lig <- which((rownames(LogNormalized.data) %in% lig)==TRUE)
varsel.recp <- which((rownames(LogNormalized.data) %in% recp)==TRUE)

lig.p <- length(varsel.lig)
recp.p <- length(varsel.recp)

dat <- cbind(t(LogNormalized.data[varsel.lig, Lig.id]), matrix(NA, Lig.n, recp.p))
colnames(dat) <- c(paste(nam[varsel.lig], ".Lig", sep = ""),
                   paste(nam[varsel.recp], ".Recp", sep = "") )
dat <- rbind(dat, cbind(matrix(NA,Recp.n, lig.p), t(LogNormalized.data[varsel.recp, Recp.id])) )
dat <- dat[,-which(colSums(dat,na.rm = TRUE) <= 0)]
dat <- as.data.frame(dat)

rm(LogNormalized.data)

for (kk in 1:20){
  obj <- rfsrc(Unsupervised() ~., data = dat, na.action = "na.impute")
  options(rf.cores = 1)
  options(mc.cores = 1)
  pmd <- find.interaction(obj, method="maxsubtree", sorted=FALSE, verbose=FALSE)

  lig.P <- length(which(grepl(".Lig",colnames(pmd))==TRUE))
  recp.P <- length(which(grepl(".Recp",colnames(pmd))==TRUE))


  save(pmd, obj, Lig.n, Recp.n, lig.P, recp.P, Lig.cell, Recp.cell, file = paste(c(c(paste("Ligand (",Lig.cell,")", sep =""),
                                                                                     paste("Receptor (", Recp.cell,")", sep =""),"PMDvalue"),kk, ".RData"), collapse = "_"))
}