diff --git a/getSignifActivities.R b/getSignifActivities.R index fefe949..d2b87eb 100644 --- a/getSignifActivities.R +++ b/getSignifActivities.R @@ -1,57 +1,59 @@ getSignifActivities <- function(E_centered, N, treatment_group, control_group) { # Estimates and returns the activities of TE subfamilies by multiple linear regression # Each activity coefficient is tested agains the null hypothesis that it is equal to zero # A Benjamini-Hochberg procedure to control the FDR is employed. # E centered is the matrix of gene expression values with genes as rows and samples as columns # N is the predictor matrix with genes as rows and TE subfamilies as columns # treatment_group is vector with the name of treated samples (correspond to columns in E_centered) # control_group: same for control samples + print('Starting new sample') + # defensive programming stopifnot(dim(E_centered)[1]==dim(N)[1]) stopifnot(all(sapply(treatment_group, function(x) x %in% colnames(E_centered)))) stopifnot(all(sapply(control_group, function(x) x %in% colnames(E_centered)))) stopifnot(all(sapply(control_group, function(x) ! x %in% treatment_group))) # generating the difference in expression vector control = c(sapply(control_group, (function(x) E_centered[, x]))) treatment = c(sapply(treatment_group, (function(x) E_centered[, x]))) delta = treatment-control # generating the augmented N matrix to match the dimensions of delta n_replicates = length(control_group) N_augmented = N i = 1 while (i < n_replicates) { N_augmented = rbind(N_augmented, N) i = i+1 } # assembling the design matrix, regression delta as a function of all columns in N_augmented design_matrix = as.data.frame(cbind(delta, N_augmented)) fit = lm(delta ~ ., data = design_matrix) coefs = as.data.frame(summary(fit)$coefficients) # FDR adjustment # ignoring the intercept (-1) p_adj = p.adjust(coefs[-1, 4], method='BH') # adding an NA to match the dimensions of `coefs` and appending p_adj = append(p_adj, NA, 0) coefs$p_adj = p_adj res = list() res$summary = summary(fit) res$coefs = coefs res$control_group = control_group res$treatment_group = treatment_group return(res) }