Merge pull request #30 from uclahs-cds/jarbet-illumina-850K-support

Support for CpGs in Illumina EPIC 850k human methylation array

Author: jarbet

DESCRIPTION

@@ -1,9 +1,9 @@
 Package: PrCaMethy
 Title: Prostate Cancer Methylation
-Version: 0.1.0
-Date: 2025-01-16
+Version: 0.2.0
+Date: 2025-05-15
 Maintainer: Jaron Arbet <[email protected]>
-Description: Resources for predicting clinical and molecular features using prostate cancer DNA methylation data generated from the Illumina 450K array.
+Description: Resources for predicting clinical and molecular features using prostate cancer DNA methylation data.
 Authors@R: c(
     person("Jaron Arbet", role = c("aut", "cre"), email = "[email protected]"),
     person("Paul C. Boutros", role = "aut"))

NEWS.md

@@ -1,3 +1,13 @@
+## PrCaMethy 0.2.0 (2025-05-15)
+
+### New Features
+
+* Support CpGs from Illumina 850K array when using `gene.methylation` to calculate gene-level methylation.  Previously only 450K array was supported.
+
+### Enhancements
+
+* Documentation updated to explain both Illumina 450K and 850K arrays are supported (previously only 450K was supported)
+
 ## PrCaMethy 0.1.0 (2025-02-12)
 
 * First release of the package.

R/gene.methylation.R

@@ -1,7 +1,8 @@
 #' Gene-level methylation
 #'
-#' Calculate gene-level methylation for a dataset containing CpGs from the Illumina 450k methylation array.
+#' Calculate gene-level methylation for a dataset containing CpGs from the Illumina 450k and 850k human methylation arrays.
 #' Gene-level methylation is calculated as the median beta-value among CpG islands in the gene promoter region.
+#' Note that gene names that originally contained a hyphen (-) are replaced with a period (.)
 #'
 #' @param methy methylation dataset where rownames give patient ids and columns use CpG ids
 #' @param print.progress TRUE/FALSE to show progress bar

R/gene.promoter.cpgi.R

@@ -1,4 +1,4 @@
 #' CpG islands in the gene promoter region
 #'
-#' List of genes where each element gives a vector of CpG islands in the promoter region of that gene.
+#' List of genes where each element gives a vector of CpG islands in the promoter region of that gene. The Illumina cg IDs come from the 450k and 850k arrays.
 'gene.promoter.cpgi'

README.md

@@ -1,6 +1,8 @@
 # PrCaMethy
 
-The `PrCaMethy` R package offers tools for predicting clinical and molecular features using prostate cancer DNA methylation data generated from the Illumina 450K array.
+The `PrCaMethy` R package offers tools for predicting clinical and molecular features using prostate cancer DNA methylation data.
+
+The user inputs CpG beta-values from CpGs contained in the Illumina 450K and/or 850K human methylation arrays.  See the tutorial below for more details.
 
 ## Tutorial
 

data/example.data.gene.methy.rda

@@ -0,0 +1,86 @@
+library(rtracklayer);
+library(GenomicRanges);
+library(R.utils);
+
+devtools::load_all();
+source('config.R') # see project PRAD-000101-MethySubtypes/PrCaMethy/config.R
+
+cpg.450k <- readRDS(arg$path.cpg.annotation.450k);
+cpg.850k <- readRDS(arg$path.cpg.annotation.850k);
+cpg.850k <- cpg.850k[!cpg.850k$cpg %in% cpg.450k$cpg, ];
+cpg.850k <- cpg.850k[!duplicated(cpg.850k$cpg), ];
+
+cpg.450k$array <- '450k';
+cpg.850k$array <- '850k';
+
+colnames(cpg.450k)[which(colnames(cpg.450k) == 'chr')] <- 'chr.hg19';
+colnames(cpg.450k)[which(colnames(cpg.450k) == 'pos')] <- 'pos.hg19';
+colnames(cpg.850k)[which(colnames(cpg.850k) == 'chr')] <- 'chr.hg38';
+colnames(cpg.850k)[which(colnames(cpg.850k) == 'pos')] <- 'pos.hg38';
+
+
+cpg.450k.sub <- cpg.450k[, c('cpg', 'chr.hg19', 'pos.hg19', 'strand', 'cpg.type', 'promoter', 'array','UCSC_RefGene_Name')];
+cpg.850k.sub <- cpg.850k[, c('cpg', 'chr.hg38', 'pos.hg38', 'strand', 'cpg.type', 'promoter', 'array','UCSC_RefGene_Name')];
+
+levels(cpg.850k.sub$cpg.type)[levels(cpg.850k.sub$cpg.type) == 'OpenSea'] <- 'Open sea';
+stopifnot(all(levels(cpg.450k.sub$cpg.type) == levels(cpg.850k.sub$cpg.type)));
+for (i in 1:ncol(cpg.450k.sub)) {
+    if (is.factor(cpg.450k.sub[, i])) {
+        #print(i);
+        stopifnot(all(levels(cpg.450k.sub[, i]) == levels(cpg.850k.sub[, i])));
+        }
+    }
+
+
+########## convert 450k coordinates to hg38
+# Create GRanges object from hg19 coordinates
+gr.hg19 <- GRanges(
+    seqnames = paste0('chr', as.character(cpg.450k.sub$chr.hg19)),
+    ranges = IRanges(start = cpg.450k.sub$pos.hg19, end = cpg.450k.sub$pos.hg19),
+    strand = cpg.450k.sub$strand,
+    cpg = cpg.450k.sub$cpg
+    );
+
+# Download the liftover chain file
+# path.file <- file.path(arg$path.save.annot, paste0(chain.file, '.gz'));
+# download.file(
+#     'http://hgdownload.cse.ucsc.edu/goldenPath/hg19/liftOver/hg19ToHg38.over.chain.gz',
+#     destfile = path.file
+#     );
+# gunzip(path.file, remove = FALSE);
+
+
+# Load the chain file and perform liftOver
+chain <- import.chain(arg$path.hg19.to.hg38.chain);
+gr.hg38 <- liftOver(gr.hg19, chain)
+gr.hg38 <- unlist(gr.hg38)
+
+# Convert back to data frame and merge with original
+lifted.df <- data.frame(
+    cpg = gr.hg38$cpg,
+    chr.hg38 = gsub('chr', '', as.character(seqnames(gr.hg38))),
+    pos.hg38 = start(gr.hg38)
+    );
+
+# Merge hg38 coordinates back into original data
+cpg.450k.hg38 <- merge(cpg.450k.sub, lifted.df, by = 'cpg', all.x = TRUE)
+mean(is.na(cpg.450k.hg38$pos.hg38));
+
+
+
+# pool 450k and 850k annotations
+cpg.annotation <- data.frame(data.table::rbindlist(list(cpg.450k.hg38, cpg.850k.sub), fill = TRUE), check.names = FALSE)
+cpg.annotation <- cpg.annotation[, c('cpg', 'chr.hg38', 'pos.hg38', 'chr.hg19', 'pos.hg19', 'cpg.type', 'promoter', 'array','UCSC_RefGene_Name')];
+mean(cpg.annotation$chr.hg38 == cpg.annotation$chr.hg19, na.rm = TRUE);
+
+saveRDS(
+    cpg.annotation,
+    file = file.path(arg$path.save.annot, paste0(Sys.Date(), '_cpg_annotation_450k-850k.rds')),
+    );
+
+data.table::fwrite(
+    cpg.annotation,
+    file = file.path(arg$path.save.annot, paste0(Sys.Date(), '_cpg_annotation_450k-850k.csv.gz')),
+    row.names = FALSE,
+    col.names = TRUE
+    );

data/gene.promoter.cpgi.rda

@@ -1,20 +1,15 @@
-# https://www.bioconductor.org/packages/release/data/annotation/html/IlluminaHumanMethylation450kanno.ilmn12.hg19.html
-# v 0.6.1
 devtools::load_all();
 library(IlluminaHumanMethylation450kanno.ilmn12.hg19);
+source('config.R') # see project PRAD-000101-MethySubtypes/PrCaMethy/config.R
 
 data(Other);
 cpg.annotation <- as.data.frame(Other);
-cpg.annotation$genomic.location <- factor(ifelse(
-    test = cpg.annotation$UCSC_RefGene_Name == '',
-    yes = 'Intergenic',
-    no = ifelse(
+cpg.annotation$promoter <- ifelse(
         test = grepl('TSS200|TSS1500|5\'UTR', cpg.annotation$UCSC_RefGene_Group) | cpg.annotation$Regulatory_Feature_Group == 'Promoter_Associated',
-        yes = 'Promoter',
-        no = 'Body'
-        )
-    ));
-table(cpg.annotation$genomic.location);
+        yes = 'yes',
+        no = 'no'
+        );
+table(cpg.annotation$promoter);
 
 cpg.annotation$cpg <- rownames(cpg.annotation);
 
@@ -61,4 +56,9 @@ cpg.annotation <- merge(
     all.x = TRUE
     );
 
-usethis::use_data(cpg.annotation, overwrite = TRUE, compress = 'xz');
+cpg.annotation.450k <- cpg.annotation;
+#usethis::use_data(cpg.annotation.450k, overwrite = TRUE, compress = 'xz');
+saveRDS(
+    cpg.annotation.450k,
+    file = file.path(arg$path.save.annot, paste0(Sys.Date(), '_cpg_annotation_450k.rds')),
+    );

inst/cpg.annotation.450k.850k.R

@@ -0,0 +1,54 @@
+devtools::load_all();
+library(IlluminaHumanMethylationEPICv2anno.20a1.hg38);
+source('config.R') # see project PRAD-000101-MethySubtypes/PrCaMethy/config.R
+
+data(Other);
+cpg.annotation <- as.data.frame(Other);
+cpg.annotation$promoter <- ifelse(
+        test = grepl('TSS200|TSS1500|5\'UTR', cpg.annotation$UCSC_RefGene_Group) | cpg.annotation$Regulatory_Feature_Group == 'Promoter_Associated',
+        yes = 'yes',
+        no = 'no'
+        );
+
+cpg.annotation$cpg.id.full <- rownames(cpg.annotation);
+cpg.annotation$cpg <- sapply(strsplit(cpg.annotation$cpg.id.full, '_'), function(x) x[1]);
+
+data(Locations);
+locations <- as.data.frame(Locations);
+locations$chr <- gsub('chr', '',locations$chr);
+stopifnot(all(as.character(c(1:22, 'X', 'Y')) %in% locations$chr));
+locations$chr <- factor(
+    x = locations$chr,
+    levels = as.character(c(1:22, 'X', 'Y'))
+    );
+locations$cpg.id.full <- rownames(locations);
+locations$cpg <- sapply(strsplit(locations$cpg.id.full, '_'), function(x) x[1]);
+stopifnot(all(locations$cpg %in% cpg.annotation$cpg));
+cpg.annotation <- merge(
+    x = cpg.annotation,
+    y = locations,
+    by = 'cpg.id.full'
+    );
+stopifnot(nrow(cpg.annotation) == nrow(locations));
+
+# island, open sea, shelf, shore
+data(Islands.UCSC);
+islands <- as.data.frame(Islands.UCSC);
+islands$cpg.type <- islands$Relation_to_Island
+islands$cpg.type <- factor(islands$cpg.type);
+islands$cpg.id.full <- rownames(islands);
+islands$cpg <- sapply(strsplit(islands$cpg.id.full, '_'), function(x) x[1]);
+cpg.annotation <- merge(
+    x = cpg.annotation,
+    y = islands,
+    by = 'cpg.id.full',
+    all.x = TRUE
+    );
+stopifnot(nrow(cpg.annotation) == nrow(islands));
+
+cpg.annotation.850k <- cpg.annotation;
+#usethis::use_data(cpg.annotation.450k, overwrite = TRUE, compress = 'xz');
+saveRDS(
+    cpg.annotation.850k,
+    file = file.path(arg$path.save.annot, paste0(Sys.Date(), '_cpg_annotation_850k.rds')),
+    );