Example_script.R

In FedericoMelograna/ISN: ISN modeling

source("R/Script_ISNs_library.R")
source("R/simulation_data")
# data path: where I import the data
data_path = "..."
# Result path: where I save the file
result_path = "..."


# First, generate simulation data -----------------------------------------

dataset_sim = data_simulation(nind = 20,ngeni = 10, seed = 16)
write.table(dataset_sim, file = "example_dataset.txt",quote = F, sep = "\t", row.names = T, col.names = T)


# Import the data and start with the analysis -----------------------------

dataset_sim = read.table(file = "example_dataset.txt")
head(dataset_sim)
t_data_sim = t(dataset_sim)


# Liones wants data in form features x individuals



# Prepare summarized experiment -------------------------------------------


rowData <- S4Vectors::DataFrame(row.names = rownames(t_data_sim), gene = rownames(t_data_sim))
colData <- S4Vectors::DataFrame(col.names = colnames(t_data_sim), sample =colnames(t_data_sim))

se <- SummarizedExperiment::SummarizedExperiment(assays = list(counts = t_data_sim), 
                                                 colData = colData, rowData = rowData)




# Apply functions ---------------------------------------------------------

# Depending on you data you can go for a correlation function or a bicor function with a prespecified esp

# WGCNA
esp = 2
lionessResults  <- Lioness_computer(se, f = wgcna_bicor_function(x,esp))

# Correlation
lionessResults  <- Lioness_computer(se, f = cor_function)




# Analyze results ---------------------------------------------------------


summary(lionessResults)
all_assays = assays(lionessResults)

# LOO network - intermediate step 
LOO_net = all_assays$perturbed

# Individual-specific network as in LIONESS 
ISN_net = all_assays$lioness

# Statistics
dim(ISN_net)
head(str(ISN_net))
head(rownames(ISN_net) )
ISN_net[1:5,1:5]

# Structure: in each column there are the individuals and in each row are depicted
# the individual-specific interaction for given individual. 



# Eliminate duplicate rows ------------------------------------------------

names_r_net = rownames(ISN_net)
to_eliminate_r = eliminate_repeated_row(names_r_net)

ISN_net = ISN_net[-to_eliminate_r,]
LOO_net = LOO_net[-to_eliminate_r,]



# GLOBAL network calculation --------------------------------------


net = wgcna_bicor_function(x,esp)(t_data_sim)
net = cor_function(t_data_sim)




# Trasforming global net as LOO net and ISN net
lionessOutput <- matrix(NA, nrow(net) * ncol(net),3)
lionessOutput[, 1] <- rep(row.names(net), ncol(net))
lionessOutput[, 2] <- rep(colnames(net), each = nrow(net))
lionessOutput <- as.data.frame(lionessOutput, stringsAsFactors = FALSE)
lionessOutput[,3] <- c(net)
global_net = data.frame("global" = lionessOutput[,3] )
rownames(global_net) = paste(lionessOutput[,1], lionessOutput[,2], sep = "_")



# Writing results ---------------------------------------------------------


setwd(result_path)

# ISNs
write.table(ISN_net,file="ISN_net_only_correct_rows_1000.txt") # keeps the rownames

# LOO networks
write.table(LOO_net,file="LOO_net_only_correct_rows_1000.txt") # keeps the rownamesù

# GLOBAL NETWORK in tabular format
write.table(net,file="Complete_net_only_correct_rows_1000.txt") # keeps the rownames

# GLOBAL NETWORK in line format
write.table(global_net,file="Global_net_only_correct_rows_1000.txt") # keeps the rownames
ISN_net[1:5,1:5]