update limma analysis and run script

analyses/methylation_analysis/01-limma_analysis.R

@@ -6,6 +6,7 @@ suppressPackageStartupMessages({
   library(limma)
   library(optparse)
   library(data.table)
+  library(dplyr)
 })
 
 # parse command line options
@@ -22,19 +23,16 @@ output_dir <- opt$output_dir
 dir.create(output_dir, showWarnings = F, recursive = T)
 
 # read m-values
-methyl_m_values_full <- fread(opt$methyl_mat, data.table = FALSE, nrows = 500000) %>% 
+methyl_m_values_full <- readRDS(opt$methyl_mat) %>% 
   dplyr::slice_head(n = 500000) %>%
   na.omit() %>%
-  dplyr::distinct(Probe_ID)
-
-rownames(methyl_m_values_full) <- NULL
-
+  dplyr::filter(!duplicated(Probe_ID))
 methyl_m_values_full <- methyl_m_values_full %>%
-  tibble::column_to_rownames(var = 'Probe_ID')
+  tibble::column_to_rownames('Probe_ID')
 
 # read annotation
 methyl_annot_full <- data.table::fread(opt$methyl_annot) %>%
-  dplyr::distinct('Probe_ID', .keep_all = TRUE)
+  dplyr::filter(!duplicated(Probe_ID))
 
 # create generalized function for gene feature analysis
 run_analysis <- function(methyl_m_values_full,

analyses/methylation_analysis/run_analysis.sh

@@ -15,7 +15,7 @@ echo "R_MAX_VSIZE=100Gb" > .Renviron
 # Define directory and input files
 data_dir="/sbgenomics/project-files/opc-v15"
 count_file="${data_dir}/v15/gene-counts-rsem-expected_count-collapsed.rds"
-methyl_m_file="${data_dir}/v15/methyl-m-values.rds" 
+methyl_m_file="${data_dir}/v15/methyl-m-values-hgat.rds" 
 cluster_file="../intNMF/results/intnmf_clusters.tsv"
 methyl_annot_file="${data_dir}/v15/infinium.gencode.v39.probe.annotations.tsv.gz"