Load packages

library(viper)
library(aracne.networks)
library(dplyr)
library(plyr)
library(stringr)
library(Biobase)
library(EnsDb.Hsapiens.v86)
dir.create("output/DME/",showWarnings = FALSE)
dir.create("output/methy/",showWarnings = FALSE)

Define functions

read_exp <- function(file_name) {
    expr <- read.table(file_name, header = TRUE, sep = "\t", row.names = 1,
        as.is = TRUE)
    expr[is.na(expr)] <- 0
    n_0 <- count_zeros_in_rows(as.matrix(expr))
    # Remove features with more 20% zero/missing values
    expr <- expr[n_0 >= ncol(expr)/5, ]
    # Split rownames by |
    meta <- data.frame(rownames(expr))
    colnames(meta) <- c("ENSG")
    meta$ENSG <- sub("\\..*$", "", meta$ENSG)

    geneID <- ensembldb::select(EnsDb.Hsapiens.v86, keys = meta$ENSG, keytype = "GENEID",
        columns = c("SYMBOL", "UNIPROTID", "GENEID"))
    meta$UNIPROTID <- plyr::mapvalues(meta$ENSG, from = geneID$GENEID,
        to = geneID$UNIPROTID, warn_missing = FALSE)
    meta$SYMBOL <- plyr::mapvalues(meta$ENSG, from = geneID$GENEID, to = geneID$SYMBOL,
        warn_missing = FALSE)

    # Remove duplicated indexs
    binary_unique_index <- (!duplicated(meta$SYMBOL)) & (!is.na(meta$SYMBOL) &
        (!sapply(meta$SYMBOL, function(x) startsWith(x, "ENSG"))))

    print(table(binary_unique_index))

    meta <- meta[binary_unique_index, ]

    expr <- expr[binary_unique_index, ]
    rownames(expr) <- meta$SYMBOL

    rownames(meta) <- rownames(expr)
    return(list(expr, meta))
}

calculate_log_fold_change_and_pvalue <- function(data_matrix, group, adjust_method = "BH") {
    # Check if the length of the group variable matches the number of
    # columns in the data matrix
    if (length(group) != ncol(data_matrix)) {
        stop("The length of the group variable must match the number of columns in the data matrix.")
    }
    # Ensure the group variable contains only 'normal' and 'tumor'
    if (!all(group %in% c("normal", "tumor"))) {
        stop("The group variable must only contain 'normal' and 'tumor' values.")
    }

    # Calculate the mean for each row in the tumor and normal groups
    mean_tumor <- rowMeans(data_matrix[, group == "tumor"], na.rm = TRUE)
    mean_normal <- rowMeans(data_matrix[, group == "normal"], na.rm = TRUE)

    # Calculate the fold change
    fold_change <- mean_tumor - mean_normal

    # Initialize a vector to store p-values
    p_values <- numeric(nrow(data_matrix))

    # Perform t-test for each row
    for (i in 1:nrow(data_matrix)) {
        normal_values <- data_matrix[i, group == "normal"]
        tumor_values <- data_matrix[i, group == "tumor"]
        wilcox_test_result <- wilcox.test(normal_values, tumor_values,
            paired = FALSE)
        p_values[i] <- wilcox_test_result$p.value
    }

    # Adjust the p-values
    adjusted_p_values <- p.adjust(p_values, method = adjust_method)

    # Create a data frame with fold change, p-values, and adjusted
    # p-values
    results <- data.frame(Fold_Change = fold_change, P_Value = p_values,
        Adjusted_P_Value = adjusted_p_values)
    # Return the results data frame
    return(results)
}

count_zeros_in_rows <- function(mat) {
    # Ensure the input is a matrix
    if (!is.matrix(mat)) {
        stop("Input must be a matrix.")
    }

    # Use rowSums to count zeros in each row
    zero_counts <- rowSums(mat != 0)

    return(zero_counts)
}

Extract and process Methylation data

df <- read.csv("data/omics_regulon_pairs.csv")
labels <- c("kirc", "kirc", "hnsc", "hnsc", "lusc", "lusc", "luad", "luad",
    "paad", "paad")
for (i in c(1, 3, 5, 7, 9)) {
    normal <- read_exp(df$methy[i])
    expr_n <- normal[[1]]
    meta_n <- normal[[2]]

    tumor <- read_exp(df$methy[i+1])
    expr_t <- tumor[[1]]
    meta_t <- tumor[[2]]

    common_terms <- intersect(rownames(expr_t), rownames(expr_n))

    expr_t[!is.finite(as.matrix(expr_t))] <- 0
    expr_n[!is.finite(as.matrix(expr_n))] <- 0

    expr <- cbind(expr_t[common_terms, ], expr_n[common_terms, ])

    saveRDS(expr_n[common_terms, ], paste("output/methy/count_matrix_", labels[i], "_normal.RDS",
        sep = ""))
    saveRDS(expr_t[common_terms, ], paste("output/methy/count_matrix_", labels[i], "_tumor.RDS",
        sep = ""))

    meta <- meta_n[common_terms, ]

    fc <- calculate_log_fold_change_and_pvalue(data_matrix = as.matrix(expr),
        group = c(rep("tumor", ncol(expr_t)), rep("normal", ncol(expr_n))))
    fc <- cbind(meta, fc)
    write.csv(fc, paste("output/DME/", labels[i],
        "_fc.csv", sep = ""), row.names = F)
    fc <- fc[(fc$Adjusted_P_Value < 0.05)&(abs(fc$Fold_Change)>0.1), ]
    write.csv(fc, paste("output/DME/", labels[i],
        "_fc_0.05.csv", sep = ""), row.names = F)
}

binary_unique_index
FALSE  TRUE 
  115 12841 
binary_unique_index
FALSE  TRUE 
  115 12843 
binary_unique_index
FALSE  TRUE 
  115 12842 
binary_unique_index
FALSE  TRUE 
  115 12843 
binary_unique_index
FALSE  TRUE 
  115 12842 
binary_unique_index
FALSE  TRUE 
  115 12842 
binary_unique_index
FALSE  TRUE 
  115 12840 
binary_unique_index
FALSE  TRUE 
  115 12843 
binary_unique_index
FALSE  TRUE 
  115 12842 
binary_unique_index
FALSE  TRUE 
  115 12843 

Session information

sessionInfo()

R version 4.4.0 (2024-04-24)
Platform: aarch64-apple-darwin20
Running under: macOS Sonoma 14.5

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: America/New_York
tzcode source: internal

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] EnsDb.Hsapiens.v86_2.99.0 ensembldb_2.28.0         
 [3] AnnotationFilter_1.28.0   GenomicFeatures_1.56.0   
 [5] AnnotationDbi_1.66.0      GenomicRanges_1.56.1     
 [7] GenomeInfoDb_1.40.1       IRanges_2.38.0           
 [9] S4Vectors_0.42.0          stringr_1.5.1            
[11] plyr_1.8.9                dplyr_1.1.4              
[13] aracne.networks_1.30.0    viper_1.38.0             
[15] Biobase_2.64.0            BiocGenerics_0.50.0      

loaded via a namespace (and not attached):
  [1] bitops_1.0-7                DBI_1.2.3                  
  [3] rlang_1.1.4                 magrittr_2.0.3             
  [5] git2r_0.33.0                matrixStats_1.3.0          
  [7] e1071_1.7-14                compiler_4.4.0             
  [9] RSQLite_2.3.7               png_0.1-8                  
 [11] vctrs_0.6.5                 ProtGenerics_1.36.0        
 [13] pkgconfig_2.0.3             crayon_1.5.2               
 [15] fastmap_1.2.0               XVector_0.44.0             
 [17] utf8_1.2.4                  Rsamtools_2.20.0           
 [19] promises_1.3.0              rmarkdown_2.27             
 [21] UCSC.utils_1.0.0            purrr_1.0.2                
 [23] bit_4.0.5                   xfun_0.45                  
 [25] zlibbioc_1.50.0             cachem_1.1.0               
 [27] jsonlite_1.8.8              blob_1.2.4                 
 [29] later_1.3.2                 DelayedArray_0.30.1        
 [31] BiocParallel_1.38.0         parallel_4.4.0             
 [33] R6_2.5.1                    bslib_0.7.0                
 [35] stringi_1.8.4               rtracklayer_1.64.0         
 [37] jquerylib_0.1.4             SummarizedExperiment_1.34.0
 [39] Rcpp_1.0.12                 knitr_1.47                 
 [41] mixtools_2.0.0              httpuv_1.6.15              
 [43] Matrix_1.7-0                splines_4.4.0              
 [45] tidyselect_1.2.1            abind_1.4-5                
 [47] rstudioapi_0.16.0           yaml_2.3.8                 
 [49] codetools_0.2-20            curl_5.2.1                 
 [51] lattice_0.22-6              tibble_3.2.1               
 [53] KEGGREST_1.44.0             evaluate_0.24.0            
 [55] survival_3.7-0              proxy_0.4-27               
 [57] kernlab_0.9-32              Biostrings_2.72.1          
 [59] pillar_1.9.0                MatrixGenerics_1.16.0      
 [61] KernSmooth_2.23-24          plotly_4.10.4              
 [63] generics_0.1.3              rprojroot_2.0.4            
 [65] RCurl_1.98-1.14             ggplot2_3.5.1              
 [67] munsell_0.5.1               scales_1.3.0               
 [69] class_7.3-22                glue_1.7.0                 
 [71] lazyeval_0.2.2              tools_4.4.0                
 [73] BiocIO_1.14.0               data.table_1.15.4          
 [75] GenomicAlignments_1.40.0    fs_1.6.4                   
 [77] XML_3.99-0.16.1             grid_4.4.0                 
 [79] tidyr_1.3.1                 colorspace_2.1-0           
 [81] nlme_3.1-165                GenomeInfoDbData_1.2.12    
 [83] restfulr_0.0.15             cli_3.6.2                  
 [85] workflowr_1.7.1             fansi_1.0.6                
 [87] S4Arrays_1.4.1              segmented_2.1-0            
 [89] viridisLite_0.4.2           gtable_0.3.5               
 [91] sass_0.4.9                  digest_0.6.35              
 [93] SparseArray_1.4.8           rjson_0.2.21               
 [95] htmlwidgets_1.6.4           memoise_2.0.1              
 [97] htmltools_0.5.8.1           lifecycle_1.0.4            
 [99] httr_1.4.7                  bit64_4.0.5                
[101] MASS_7.3-61