.Rhistory

saveWorkbook(wb, file = "../results/TableS1.xlsx")
options(java.parameters = "-Xmx4g")
library(xlsx)
res <- readRDS("../tmp/SteleMicroarray.rda")
#write.table(res$necorrHub_nodes, "tmp.txt", sep = "\t")
wb <- createWorkbook()
sheet1 <- createSheet(wb, sheetName = "necorrHub_nodes")
sheet2 <- createSheet(wb, sheetName = "necorrReg")
sheet3 <- createSheet(wb, sheetName = "necorrEdges")
sheet4 <- createSheet(wb, sheetName = "netstat")
sheet5 <- createSheet(wb, sheetName = "coexpres")
sheet6 <- createSheet(wb, sheetName = "hub.m.param")
addDataFrame(res$necorrHub_nodes, sheet = sheet1)
rm sheet1
wb <- createWorkbook()
sheet1 <- createSheet(wb, sheetName = "necorrHub_nodes")
sheet2 <- createSheet(wb, sheetName = "necorrReg")
sheet3 <- createSheet(wb, sheetName = "necorrEdges")
sheet4 <- createSheet(wb, sheetName = "netstat")
sheet5 <- createSheet(wb, sheetName = "coexpres")
sheet6 <- createSheet(wb, sheetName = "hub.m.param")
addDataFrame(res$necorrHub_nodes, sheet = sheet1)
addDataFrame(res$necorrReg, sheet = sheet2)
addDataFrame(res$necorrEdges, sheet = sheet3)
options(java.parameters = "-Xmx4g")
library(xlsx)
wb <- createWorkbook()
sheet1 <- createSheet(wb, sheetName = "necorrHub_nodes")
sheet2 <- createSheet(wb, sheetName = "necorrReg")
sheet3 <- createSheet(wb, sheetName = "necorrEdges")
sheet4 <- createSheet(wb, sheetName = "netstat")
sheet5 <- createSheet(wb, sheetName = "coexpres")
sheet6 <- createSheet(wb, sheetName = "hub.m.param")
addDataFrame(res$necorrHub_nodes, sheet = sheet1)
SteleHub_array <- res$necorrHub_nodes
SteleHub_array <- SteleHub_array[grep("AT[M|C|[1-5]G\\d{5}", rownames(SteleHub_array)),]
# the reg genes from the stele microarray
SteleReg_array <- res$necorrReg
SteleReg_array <- SteleReg_array[grep("AT[M|C|[1-5]G\\d{5}", rownames(SteleReg_array)),]
write.csv(SteleHub_array, "../results/Microarray_Stele_Hub.csv")
write.csv(SteleReg_array, "../results/Microarray_Stele_Reg.csv")
POanatomical <- read.delim("~/Documents/DataAnalysis/NECorr_analysis/data/po_anatomy_gene_arabidopsis_tair.assoc",
skip = 1, header = F)
POanatomical <- read.delim("~/Documents/1_projects/NECorr_analysis/data/po_anatomy_gene_arabidopsis_tair.assoc",
skip = 1, header = F)
POtime <- read.delim("~/Documents/1_projects/NECorr_analysis/data/po_temporal_gene_arabidopsis_tair.assoc",
skip = 1, header = F)
# extract Phloem annotation from the data
PO1 <- POanatomical[grep("phloem|PHLOEM|Phloem|xylem|XYLEM|vascular|VASCULAR", POanatomical$V11),]
PO2 <- POtime[grep("phloem|PHLOEM|Phloem|xylem|XYLEM|vascular|VASCULAR", POtime$V11),]
# get all genes with annotation that are linked to phloem
vascular.annot.genes <- noquote(unique(c(PO1$V10, PO2$V10)))
##
vascular.GO <- read.delim("../data/Ath_vascular_TAIR.txt", header = F)
vascular.GO.genes <- noquote(unique(vascular.GO$V1))
vascular.genes <- unique(c(vascular.GO.genes, vascular.annot.genes))
vascular.genes <- vascular.genes[grep("^AT", vascular.genes)]
library(ROCR)
MAHub <-SteleHub_array
#MAHub <- resPhloem_meristem$necorrHub_nodes
MAHub <- MAHub[order(MAHub$Gene_rank_h, decreasing = T),]
View(MAHub)
library(ROCR)
MAHub <-SteleHub_array
#MAHub <- resPhloem_meristem$necorrHub_nodes
MAHub <- MAHub[order(MAHub$Score, decreasing = T),]
MAHub <- MAHub[which(MAHub$Score > 0),]
label <- rep(0, nrow(MAHub))
label[which(rownames(MAHub) %in% vascular.genes)] <- 1
# calculate the AUC
pred <- prediction(MAHub$Score, label)
roc.perf = performance(pred, measure = "tpr", x.measure = "fpr")
auc_ROCR <- performance(pred, measure = "auc")
auc_ROCR <- auc_ROCR@y.values[[1]]
#plot
#png("../results/ROC_AUC_validation.png", res=120)
plot(roc.perf, measure = "tpr", x.measure = "fpr",
avg= "threshold",
colorize=FALSE,
lwd= 1,
main= "NECorr Stele data ")
abline(a=0, b= 1,col="red")
text( 7, 4, x = .75, y = .25, label = paste("AUC =", round(auc_ROCR, 2)))
#dev.off()
MAHublabel <- MAHub
MAHublabel$label <- label
head(MAHublabel)
MAReg <- res$necorrReg
labelReg <- rep(0, nrow(MAReg))
labelReg[which(rownames(MAReg) %in% vascular.genes)] <- 1
MAReg$labelReg <- labelReg
head(MAReg)
Edgesrank <- res$necorrEdges
head(Edgesrank)
# study the start of the ROC curve
library(ROCR)
# Assuming MAHub and label are already defined as in your initial code
# Calculate prediction and performance
pred <- prediction(MAHub$Score, label)
roc.perf <- performance(pred, measure = "tpr", x.measure = "fpr")
# Extract FPR and TPR values
fpr_values <- unlist(roc.perf@x.values)
tpr_values <- unlist(roc.perf@y.values)
# Find the TPR corresponding to the closest FPR to 0.1
index_closest <- which.min(abs(fpr_values - 0.1))
average_tpr_at_0.1 <- tpr_values[index_closest]
# Plot the ROC curve with highlighted TPR at FPR = 0.1
plot(roc.perf, main = "NECorr Stele data")
abline(a = 0, b = 1, col = "red")
points(fpr_values[index_closest], tpr_values[index_closest], col = "blue", pch = 19)
text(fpr_values[index_closest], tpr_values[index_closest], labels = paste("TPR at FPR=0.1: ", round(average_tpr_at_0.1, 2)), pos = 3)
# Add AUC text
auc_ROCR <- performance(pred, measure = "auc")
auc_ROCR <- auc_ROCR@y.values[[1]]
text(x = 0.75, y = 0.25, labels = paste("AUC =", round(auc_ROCR, 2)))
library(ROCR)
# Assuming MAHub and label are already defined as in your initial code
# Calculate prediction and performance
pred <- prediction(MAHub$Score, label)
roc.perf <- performance(pred, measure = "tpr", x.measure = "fpr")
# Extract FPR and TPR values
fpr_values <- unlist(roc.perf@x.values)
tpr_values <- unlist(roc.perf@y.values)
# Find the TPR corresponding to the closest FPR to 0.1
index_closest <- which.min(abs(fpr_values - 0.1))
average_tpr_at_0.1 <- tpr_values[index_closest]
# Plot the ROC curve with a different color and zoom into the steep increase
plot(roc.perf, col = "blue", main = "NECorr Stele data", xlim = c(0, 0.1), ylim = c(0, max(tpr_values[1:index_closest])))
abline(a = 0, b = 1, col = "red", lty = 2) # Adding a dashed line for reference
points(fpr_values[index_closest], tpr_values[index_closest], col = "darkgreen", pch = 19)
text(fpr_values[index_closest], tpr_values[index_closest], labels = paste("TPR at FPR=0.1: ", round(average_tpr_at_0.1, 2)), pos = 4)
# Add AUC text
auc_ROCR <- performance(pred, measure = "auc")
auc_ROCR <- auc_ROCR@y.values[[1]]
text(x = 0.05, y = 0.9, labels = paste("AUC =", round(auc_ROCR, 2)))
# Calculate the number of true positives with almost no true negatives
threshold_fpr <- 0.01 # Define a threshold for "almost no true negatives"
true_positives_initial <- sum(tpr_values[fpr_values <= threshold_fpr])
# Print the number of true positive genes with almost no true negatives
cat("Number of true positive genes with FPR <=", threshold_fpr, ":", true_positives_initial, "\n")
library(ROCR)
# Assuming MAHub and label are already defined as in your initial code
# Calculate prediction and performance
pred <- prediction(MAHub$Score, label)
roc.perf <- performance(pred, measure = "tpr", x.measure = "fpr")
# Extract FPR and TPR values
fpr_values <- unlist(roc.perf@x.values)
tpr_values <- unlist(roc.perf@y.values)
# Find the TPR corresponding to the closest FPR to 0.1
index_closest <- which.min(abs(fpr_values - 0.1))
average_tpr_at_0.1 <- tpr_values[index_closest]
# Set up a multi-panel plot
par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
# Full ROC curve plot
plot(roc.perf, col = "blue", main = "Full ROC Curve", xlab = "False Positive Rate", ylab = "True Positive Rate")
abline(a = 0, b = 1, col = "red", lty = 2)
text(x = 0.75, y = 0.25, labels = paste("AUC =", round(auc_ROCR, 2)))
points(fpr_values[index_closest], tpr_values[index_closest], col = "darkgreen", pch = 19)
# Zoomed-in ROC curve plot
plot(roc.perf, col = "blue", main = "Zoomed-in ROC Curve", xlab = "False Positive Rate", ylab = "True Positive Rate", xlim = c(0, 0.1), ylim = c(0, max(tpr_values[1:index_closest])))
abline(a = 0, b = 1, col = "red", lty = 2)
points(fpr_values[index_closest], tpr_values[index_closest], col = "darkgreen", pch = 19)
text(fpr_values[index_closest], tpr_values[index_closest], labels = paste("TPR at FPR=0.1: ", round(average_tpr_at_0.1, 2)), pos = 4)
# Calculate the number of true positives with almost no true negatives
threshold_fpr <- 0.01 # Define a threshold for "almost no true negatives"
true_positives_initial <- sum(tpr_values[fpr_values <= threshold_fpr])
# Print the number of true positive genes with almost no true negatives
cat("Number of true positive genes with FPR <=", threshold_fpr, ":", true_positives_initial, "\n")
library(ROCR)
# Assuming MAHub and label are already defined as in your initial code
# Calculate prediction and performance
pred <- prediction(MAHub$Score, label)
roc.perf <- performance(pred, measure = "tpr", x.measure = "fpr")
# Extract FPR and TPR values
fpr_values <- unlist(roc.perf@x.values)
tpr_values <- unlist(roc.perf@y.values)
# Find the TPR corresponding to the closest FPR to 0.1
index_closest <- which.min(abs(fpr_values - 0.1))
average_tpr_at_0.1 <- tpr_values[index_closest]
# Set up a multi-panel plot
par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
# Full ROC curve plot
plot(roc.perf, col = "blue", main = "Full ROC Curve", xlab = "False Positive Rate", ylab = "True Positive Rate")
abline(a = 0, b = 1, col = "red", lty = 2)
text(x = 0.75, y = 0.25, labels = paste("AUC =", round(auc_ROCR, 2)))
points(fpr_values[index_closest], tpr_values[index_closest], col = "darkgreen", pch = 19)
# Zoomed-in ROC curve plot
plot(roc.perf, col = "blue", main = "Zoomed-in ROC Curve", xlab = "False Positive Rate", ylab = "True Positive Rate", xlim = c(0, 0.1), ylim = c(0, max(tpr_values[1:index_closest])))
abline(a = 0, b = 1, col = "red", lty = 2)
points(fpr_values[index_closest], tpr_values[index_closest], col = "darkgreen", pch = 19)
text(fpr_values[index_closest], tpr_values[index_closest], labels = paste("TPR at FPR=0.1: ", round(average_tpr_at_0.1, 2)), pos = 4)
# Calculate the number of true positives with almost no true negatives
threshold_fpr <- 0.01 # Define a threshold for "almost no true negatives"
true_positives_initial <- sum(tpr_values[fpr_values <= threshold_fpr])
# Print the number of true positive genes with almost no true negatives
cat("Number of true positive genes with FPR <=", threshold_fpr, ":", true_positives_initial, "\n")
# Calculate the number of true positives with almost no true negatives
threshold_fpr <- 0.01 # Define a threshold for "almost no true negatives"
true_positives_initial <- sum(tpr_values[fpr_values <= threshold_fpr])
false_positives_initial <- sum(fpr_values[index_threshold] * length(label[label == 0]))
# Define the threshold for "almost no true negatives" (e.g., FPR <= 0.01)
threshold_fpr <- 0.01
# Calculate the number of true positives and false positives at the threshold FPR
index_threshold <- which(fpr_values <= threshold_fpr)
true_positives_initial <- sum(tpr_values[index_threshold] * length(label[label == 1]))
false_positives_initial <- sum(fpr_values[index_threshold] * length(label[label == 0]))
# Print the number of true positive genes with almost no true negatives
cat("Number of true positive genes with FPR <=", threshold_fpr, ":", true_positives_initial, "\n")
cat("Number of false positive genes with FPR <=", threshold_fpr, ":", false_positives_initial, "\n")
# Define the threshold for "almost no true negatives" (e.g., FPR <= 0.01)
threshold_fpr <- 0.001
# Calculate the number of true positives and false positives at the threshold FPR
index_threshold <- which(fpr_values <= threshold_fpr)
true_positives_initial <- sum(tpr_values[index_threshold] * length(label[label == 1]))
false_positives_initial <- sum(fpr_values[index_threshold] * length(label[label == 0]))
# Print the number of true positive genes with almost no true negatives
cat("Number of true positive genes with FPR <=", threshold_fpr, ":", true_positives_initial, "\n")
cat("Number of false positive genes with FPR <=", threshold_fpr, ":", false_positives_initial, "\n")
# Define the threshold for "almost no true negatives" (e.g., FPR <= 0.01)
threshold_fpr <- 0.0001
# Calculate the number of true positives and false positives at the threshold FPR
index_threshold <- which(fpr_values <= threshold_fpr)
true_positives_initial <- sum(tpr_values[index_threshold] * length(label[label == 1]))
false_positives_initial <- sum(fpr_values[index_threshold] * length(label[label == 0]))
# Print the number of true positive genes with almost no true negatives
cat("Number of true positive genes with FPR <=", threshold_fpr, ":", true_positives_initial, "\n")
cat("Number of false positive genes with FPR <=", threshold_fpr, ":", false_positives_initial, "\n")
# Define the threshold for "almost no true negatives" (e.g., FPR <= 0.01)
threshold_fpr <- 0.001
# Calculate the number of true positives and false positives at the threshold FPR
index_threshold <- which(fpr_values <= threshold_fpr)
true_positives_initial <- sum(tpr_values[index_threshold] * length(label[label == 1]))
false_positives_initial <- sum(fpr_values[index_threshold] * length(label[label == 0]))
# Print the number of true positive genes with almost no true negatives
cat("Number of true positive genes with FPR <=", threshold_fpr, ":", true_positives_initial, "\n")
cat("Number of false positive genes with FPR <=", threshold_fpr, ":", false_positives_initial, "\n")
categories <- c("TF", "Kinase", "Enzyme", "Receptor", "Other")
description.file <- data.table(
gene_id = genes,
Category = sample(categories, 50, replace = TRUE),
Description = paste("Description for", genes)
)
# ---- 1. Expression data ----
genes <- paste0("Gene", 1:50)
categories <- c("TF", "Kinase", "Enzyme", "Receptor", "Other")
description.file <- data.table(
gene_id = genes,
Category = sample(categories, 50, replace = TRUE),
Description = paste("Description for", genes)
)
View(description.file)
View(description.file)
getwd()
description_df <- fread("data/Ath.GeneDesc_20221123.csv", header = TRUE, sep = ",",
nThread = n_threads, showProgress = FALSE)
description_df <- fread("data/Ath.GeneDesc_20221123.csv", header = TRUE, sep = ",",
nThread = 2, showProgress = FALSE)
View(description_df)
desc_ids <- description_df[[1]]
desc_ids
description_df[, 1 := NULL]
description_df[, Gene := desc_ids]
setkey(description_df, Gene)
rm(list=ls())
# ============================
# NECorr Package Test Script
# ============================
remove.packages("NECorr")
# rm -rf /Users/liseronc/Library/R/arm64/4.3/library/NECorr
setwd("~/Documents/1_projects/NECorr/")
devtools::clean_dll()
#devtools::compileAttributes()
devtools::check()
devtools::install()
set.seed(21)
library(NECorr)
setwd("/Users/liseronc/Documents/1_projects/NECorr_analysis")
resPhloem_meristem <- NECorr(
networkFile     = "data/metanetwork_20250909.txt",
expression      = "data/pseudobulk_At_root_norm.txt",
description.file= "data/Ath.GeneDesc_20221123.csv",
condition       = "Meristem_Phloem",
meta_data       = "data/metadata_20221122.txt",
method = c("GCC"),
visualize = TRUE,
permutation = 1000,
save_results = TRUE,
interactive_net = FALSE,
output_dir = "results/results_Phloem_meristem_2025_09_09",
useBestGCC = TRUE,
asymmetricGCC = FALSE
)
saveRDS(resPhloem_meristem, "/Users/liseronc/Documents/1_projects/NECorr_analysis/tmp/resPhloem.rda")
results <- NECorr(
networkFile     = necorr_example$network,
expression      = necorr_example$expression,
description.file= necorr_example$description,
condition       = "f1",
meta_data       = necorr_example$meta_data,
permutation     = 100,
save_results = FALSE,
interactive_net = FALSE
)
print(results)
# ============================
# NECorr Package Test Script
# ============================
remove.packages("NECorr")
rm(list=ls())
# rm -rf /Users/liseronc/Library/R/arm64/4.3/library/NECorr
setwd("~/Documents/1_projects/NECorr/")
# rm -rf /Users/liseronc/Library/R/arm64/4.3/library/NECorr
setwd("~/Documents/1_projects/NECorr/")
devtools::clean_dll()
#devtools::compileAttributes()
devtools::check()
devtools::install()
# 1. Load your package
library(NECorr)
#run the tutorial
#vignette("NECorr_tutorial")
# 2. Load example data
source("data-raw/make_necorr_example.R")
data(necorr_example)
# 3. Run NECorr
results <- NECorr(
networkFile     = necorr_example$network,
expression      = necorr_example$expression,
description.file= necorr_example$description,
condition       = "f1",
meta_data       = necorr_example$meta_data,
permutation     = 100,
save_results = FALSE,
interactive_net = FALSE
)
print(results)
remove.packages("NECorr")
rm(list=ls())
devtools::clean_dll()
devtools::install()
# 1. Load your package
library(NECorr)
source("data-raw/make_necorr_example.R")
data(necorr_example)
# 3. Run NECorr
results <- NECorr(
networkFile     = necorr_example$network,
expression      = necorr_example$expression,
description.file= necorr_example$description,
condition       = "f1",
meta_data       = necorr_example$meta_data,
permutation     = 100,
save_results = FALSE,
interactive_net = FALSE
)
print(results)
set.seed(21)
library(NECorr)
setwd("/Users/liseronc/Documents/1_projects/NECorr_analysis")
resPhloem_meristem <- NECorr(
networkFile     = "data/metanetwork_20250909.txt",
expression      = "data/pseudobulk_At_root_norm.txt",
description.file= "data/Ath.GeneDesc_20221123.csv",
condition       = "Meristem_Phloem",
meta_data       = "data/metadata_20221122.txt",
method = c("GCC"),
visualize = TRUE,
permutation = 1000,
save_results = TRUE,
interactive_net = FALSE,
output_dir = "results/results_Phloem_meristem_2025_09_09",
useBestGCC = TRUE,
asymmetricGCC = FALSE
)
remove.packages("NECorr")
rm(list=ls())
setwd("~/Documents/1_projects/NECorr/")
devtools::clean_dll()
# !! restart session within the RStudio IDE (Session > Restart R)
devtools::install()
library(NECorr)
set.seed(21)
library(NECorr)
setwd("/Users/liseronc/Documents/1_projects/NECorr_analysis")
resPhloem_meristem <- NECorr(
networkFile     = "data/metanetwork_20250909.txt",
expression      = "data/pseudobulk_At_root_norm.txt",
description.file= "data/Ath.GeneDesc_20221123.csv",
condition       = "Meristem_Phloem",
meta_data       = "data/metadata_20221122.txt",
method = c("GCC"),
visualize = TRUE,
permutation = 1000,
save_results = TRUE,
interactive_net = FALSE,
output_dir = "results/results_Phloem_meristem_2025_09_09",
useBestGCC = TRUE,
asymmetricGCC = FALSE
)
saveRDS(resPhloem_meristem, "/Users/liseronc/Documents/1_projects/NECorr_analysis/tmp/resPhloem.rda")
library(readxl)
flower1 <- read_excel("~/Documents/1_projects/NECorr_analysis/data/eraa057_suppl_supplementary_tables_s1_s6.xlsx",
sheet = "Suppl Table 1", skip = 10)
flower2 <- read_excel("~/Documents/1_projects/NECorr_analysis/data/eraa057_suppl_supplementary_tables_s1_s6.xlsx",
sheet = "Suppl Table 2", skip = 1)
flower3 <- read_excel("~/Documents/1_projects/NECorr_analysis/data/eraa057_suppl_supplementary_tables_s1_s6.xlsx",
sheet = "Suppl Table 3", skip = 1)
flower4 <- read_excel("~/Documents/1_projects/NECorr_analysis/data/eraa057_suppl_supplementary_tables_s1_s6.xlsx",
sheet = "Suppl Table 4", skip = 1)
flower5 <- read_excel("~/Documents/1_projects/NECorr_analysis/data/eraa057_suppl_supplementary_tables_s1_s6.xlsx",
sheet = "Suppl Table 5", skip = 1)
flower6 <- read_excel("~/Documents/1_projects/NECorr_analysis/data/eraa057_suppl_supplementary_tables_s1_s6.xlsx",
skip = 10)
View(flower3)
View(flower1)
View(flower2)
View(flower4)
View(flower5)
View(flower6)
View(flower1)
data_list <- list(flower1, flower2, flower3, flower4, flower5, flower6)
merged_data <- do.call(rbind, data_list)
View(flower3)
View(flower3)
View(flower3)
data_list <- list(flower1, flower2, flower3[,c(1..7)], flower4, flower5, flower6)
data_list <- list(flower1, flower2, flower3[,c(1:7)], flower4, flower5, flower6)
merged_data <- do.call(rbind, data_list)
flower1 <- read_excel("~/Documents/1_projects/NECorr_analysis/data/eraa057_suppl_supplementary_tables_s1_s6.xlsx",
sheet = "Suppl Table 1", skip = 10)
flower2 <- read_excel("~/Documents/1_projects/NECorr_analysis/data/eraa057_suppl_supplementary_tables_s1_s6.xlsx",
sheet = "Suppl Table 2", skip = 1)
colnames(flower2) <- colnames(flower1)
flower3 <- read_excel("~/Documents/1_projects/NECorr_analysis/data/eraa057_suppl_supplementary_tables_s1_s6.xlsx",
sheet = "Suppl Table 3", skip = 1)
flower3 <- flower3[,c(1:7)]
colnames(flower3) <- colnames(flower1)
flower4 <- read_excel("~/Documents/1_projects/NECorr_analysis/data/eraa057_suppl_supplementary_tables_s1_s6.xlsx",
sheet = "Suppl Table 4", skip = 1)
colnames(flower4) <- colnames(flower1)
flower5 <- read_excel("~/Documents/1_projects/NECorr_analysis/data/eraa057_suppl_supplementary_tables_s1_s6.xlsx",
sheet = "Suppl Table 5", skip = 1)
colnames(flower5) <- colnames(flower1)
flower6 <- read_excel("~/Documents/1_projects/NECorr_analysis/data/eraa057_suppl_supplementary_tables_s1_s6.xlsx",
skip = 10)
colnames(flower6) <- colnames(flower1)
# merge all flower files
data_list <- list(flower1, flower2, flower3[,c(1:7)], flower4, flower5, flower6)
merged_data <- do.call(rbind, data_list)
merged_unique <- merged_data[!duplicated(merged_data$`Gene number`), ]
View(merged_unique)
library(data.table)
library(matrixStats)
library(igraph)
library(parallel)
library(stats)
meta_data <- "PRJNA314076_sample_metadata.txt"
expression <- "PRJNA314076_norm_counts.txt"
eset <- fread(expression, header = TRUE, sep = "\t",
nThread = n_threads, showProgress = FALSE)
n_threads = 3
eset <- fread(expression, header = TRUE, sep = "\t",
nThread = n_threads, showProgress = FALSE)
gene_ids <- eset[[1]]
eset[, 1 := NULL]
eset[, gene_id := gene_ids]
setkey(eset, gene_id)
View(eset)
eset <- fread(expression, header = TRUE, sep = "\t",
nThread = n_threads, showProgress = FALSE)
View(eset)
fread
?fread
eset <- fread(expression, header = TRUE, sep = "\t",
nThread = n_threads, showProgress = FALSE)
View(eset)
eset <- fread(expression, sep = "\t",
nThread = n_threads, showProgress = FALSE)
eset <- fread(expression, sep = "\t",
nThread = n_threads, showProgress = FALSE)
eset <- fread(expression, header = TRUE, sep = "\t",
nThread = n_threads, showProgress = FALSE)
gene_ids <- eset[[1]]
eset[, 1 := NULL]
eset[, gene_id := gene_ids]
setkey(eset, gene_id)
factortab <- fread(meta_data, header = TRUE, sep = "\t",
nThread = n_threads, showProgress = FALSE)
# First column is sample IDs
sample_ids <- factortab[[1]]
# Remove the first column from the data.table
factortab[, 1 := NULL]
# Rename the current first column to "condition"
setnames(factortab, 1, "condition")
# Add sample_id as a separate column
factortab[, sample_id := sample_ids]
# Set rownames for internal use
rownames(factortab) <- sample_ids
View(factortab)
expr_sample_cols <- setdiff(colnames(eset), "gene_id")
if (!setequal(rownames(factortab), expr_sample_cols)) {
stop("Mismatch between metadata sample IDs and expression data column names.")
}
getwd("../NECorr_analysis/)
""
)
getwd("../NECorr_analysis/"
)
getwd("~/Documents/1_projects/NECorr_analysis/")
setwd("~/Documents/1_projects/NECorr_analysis/")
library(NECorr)
NECorr