Last active
August 24, 2025 12:30
-
-
Save alexg9010/312daf779614686aba56f14400250a83 to your computer and use it in GitHub Desktop.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| # Combine methylKit annotation outputs into a single comprehensive dataframe | |
| library(methylKit) | |
| library(genomation) | |
| ## Generating setup with test data ------- | |
| # this follows https://bioconductor.org/packages/devel/bioc/vignettes/methylKit/inst/doc/methylKit.html#4_Annotating_differentially_methylated_bases_or_regions | |
| data("methylKit") | |
| myDiff <- methylDiff.obj | |
| myDiff25p <- methylKit::getMethylDiff(myDiff, difference = 25) | |
| # read the gene BED file | |
| gene.obj=genomation::readTranscriptFeatures(system.file("extdata", "refseq.hg18.bed.txt", | |
| package = "methylKit")) | |
| diffAnn=genomation::annotateWithGeneParts(as(myDiff25p,"GRanges"),gene.obj) | |
| cpg.obj=genomation::readFeatureFlank(system.file("extdata", "cpgi.hg18.bed.txt", | |
| package = "methylKit"), | |
| feature.flank.name=c("CpGi","shores")) | |
| diffCpGann=genomation::annotateWithFeatureFlank(as(myDiff25p,"GRanges"), | |
| cpg.obj$CpGi,cpg.obj$shores, | |
| feature.name="CpGi",flank.name="shores") | |
| ## Extract data from objects ------- | |
| # Assuming you have these objects from your workflow: | |
| # myDiff - ouptut from calculateDiffMeth() | |
| # myDiff25p - output from getMethylDiff() | |
| # gene.obj - output from annotateWithGeneParts() | |
| # cpg.obj - output from annotateWithFeatureFlank() for CpG islands | |
| ## convert methylDiff to data.frame | |
| df_diff <- methylKit::getData(myDiff25p) | |
| dim(df_diff) | |
| ## get gene part annotation via getMembers function | |
| ## matrix showing overlap of target features with annotation features. | |
| ## 1 for overlap, 0 for non-overlap | |
| df_genepart_assoc <- genomation::getMembers(diffAnn) | |
| dim(df_genepart_assoc) | |
| ## get tss distance for each association via getAssociationWithTSS | |
| df_tss_assoc <- genomation::getAssociationWithTSS(diffAnn) | |
| dim(df_tss_assoc) | |
| ## get CpGi or Shore annotation via getMembers function | |
| df_cpgi_assoc <- genomation::getMembers(diffCpGann) | |
| dim(df_cpgi_assoc) | |
| ## join data.frames by column | |
| df_annot <- cbind(df_diff, | |
| df_genepart_assoc, | |
| df_tss_assoc, | |
| df_cpgi_assoc | |
| ) | |
| ## remove target.row | |
| df_annot$target.row <- NULL | |
| # Display first few rows to verify the result | |
| print("First 5 rows of comprehensive dataframe:") | |
| print(head(df_annot,n=5)) | |
| # Save summary | |
| write.csv(df_annot, "diffmeth_annotation.csv", row.names = FALSE) | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment