Welcome to sesame Discussions!

[zwdzwd]

👋 Welcome!

We’re using Discussions as a place to connect with other members of our community. We hope that you:

• Ask questions you’re wondering about.
• Share ideas.
• Engage with other community members.
• Welcome others and are open-minded. Remember that this is a community we
  build together 💪.

To get started, comment below with an introduction of yourself and tell us about what you do with this community.

[SofyaMa]

I would like to read in Data form muliple file and extract the beta values from them. So far, I have run this command : sdfs = lapply(
searchIDATprefixes(system.file("extdata", package = "sesameData")),
readIDATpair) . However : getBetas(sdfs) does not seem to work. Could you maybe help me out here? thank you very much !

[Pancreas-Pratik]

@SofyaMa

I would like to read in Data form muliple file and extract the beta values from them. So far, I have run this command : sdfs = lapply(
searchIDATprefixes(system.file("extdata", package = "sesameData")),
readIDATpair) . However : getBetas(sdfs) does not seem to work. Could you maybe help me out here? thank you very much !

If you see the output of:

sdfs = lapply(
searchIDATprefixes(system.file("extdata", package = "sesameData")),
readIDATpair)

A list of SigSet objects are made. getBetas() takes an input of one SigSet object, so you would have to run each of the SIgSet objects that are within the sdfs list separately and then cbind the the outputs to get a beta value matrix.

So you should do something like this:

library(sesame)
sdfs = lapply(
  searchIDATprefixes(system.file("extdata", package = "sesameData")),
  readIDATpair)

sample1 = getBetas(sdfs$`4207113116_A`)
sample2 = getBetas(sdfs$`4207113116_B`)
betas <- cbind(sample1, sample2)

Alternatively, or I would say more mainstream, @zwdzwd has made the openSesame() pipeline function that will do all the of the recommended pre-processing all together in one shot. So you could do something like this:

idat_dir = system.file("extdata/", package = "sesameData")
betas = openSesame(idat_dir, BPPARAM=BiocParallel::MulticoreParam(2))

[SofyaMa]

👍 thank you!!!

[SofyaMa]

@Pancreas-Pratik As you suggested, I applied the openSesame() pipeline. However now I end up with way to many NAs. I used platform =EPIC in the function, as my Data is from an EPIC array. I think something is wrong with the mask ... could you maybe think of what have gone wrong in this step or suggest what to do about it? thank you!

[Pancreas-Pratik]

You're welcome @SofyaMa

@Pancreas-Pratik As you suggested, I applied the openSesame() pipeline. However now I end up with way to many NAs. I used platform =EPIC in the function, as my Data is from an EPIC array. I think something is wrong with the mask ... could you maybe think of what have gone wrong in this step or suggest what to do about it? thank you!

From the Sesame Basic Usage vignette:

We recommend two types of probe masking:

Experiment-dependent Probe Masking based on signal detection p-value (Zhou et al. 2018). Probes with p-value higher than a threshold (default: 0.05) are masked (see following for detection p-value calculation using the pOOBAH method).

Experiment-independent Probe Masking due to design issues. This is typically designated in the mask column of the manifest (see Zhou et al. 2017): This masking supports EPIC, MM285, HM450 and HM27 and is turned on by default and can also be explicitly added by the function qualityMask.

These are the two types of masking that are done to sort of generally filter out low quality probes. You can also manually adjust the p-value cut-off (See here on the right before the Background subtraction section) so that you get less NAs.

However, you should also see here: #42 (comment)

@zwdzwd in that issue is basically just saying that as long as there is one value (non-NA value) in the groups for the comparisons you are making for differentially methylated loci (DMLs) or differentially methylated regions (DMRs), the comparison will still be made. the checkLevels() function will check for this in your SummarizedExperiment object.

See here for checkLevels() and a note by @zwdzwd:

CRITICAL: If your data is NA-masked (see the Basic Usage tutorial), it is critical to preclude CpGs with missing levels. For example, if a CpG does not contain any non-NA values for all male mice, one cannot assess sex-specific DNA methylation for this CpG. One needs to preclude such CpGs from the modeling of sex-specific differential methylation. This can be done using the checkLevels function.

In the following, we will test for both sex and tissue as co-variates using the checkLevels function:

se_ok = (checkLevels(assay(se), colData(se)$sex) &
    checkLevels(assay(se), colData(se)$tissue))
sum(se_ok)
## [1] 995
se = se[se_ok,]

Hope this helps. It took me some time to make the SummarizedExperiment object because this was new for me, but this might help.

Good luck!

[Pancreas-Pratik]

Also this one might help too: #27 (comment)

[SofyaMa]

Î would like to produce such a cell Copostion graphic that shows the estimated cellcopostion per sample. The estimateCellComposition() is not in included in the Namespace. Does someone know a work around?

[Pancreas-Pratik]

@SofyaMa you might have better luck reaching out by making an issue? Perhaps someone who is actively involved in the coding of sesame could help you?

Alternatively, you could see what the authors of that cell composition image have used to create that image, and you could see if you could either transfer your data analysis to their package or if you could emulate their coding of the function into sesame?

[Jmase0921]

Hello Sesame Community,

I hope everyone is doing well!

I am very new to R and I am trying to determine the age with some of my idats by using the predictAge. But I am receiving error messages, : object 'exit' not found or object 'predictAge' not found. Does anyone have any suggestions of how to get the age without any errors? Any information will greatly be appreciated.

library(sesame)

sdftest = openSesame("primary_dir", func=NULL, package = "sesameData")
betas <-openSesame(sdftest, mask =FALSE)
predictAge(betas)
Error in predictAge(betas) : object 'exit' not found

sessionInfo()
R version 4.3.0 (2023-04-21 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19045)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.utf8 LC_CTYPE=English_United States.utf8
[3] LC_MONETARY=English_United States.utf8 LC_NUMERIC=C
[5] LC_TIME=English_United States.utf8

time zone: America/New_York
tzcode source: internal

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

other attached packages:
[1] sesame_1.18.3 sesameData_1.18.0 ExperimentHub_2.8.0 AnnotationHub_3.8.0
[5] BiocFileCache_2.8.0 dbplyr_2.3.2 BiocGenerics_0.46.0 BiocManager_1.30.20

loaded via a namespace (and not attached):
[1] tidyselect_1.2.0 dplyr_1.1.2
[3] blob_1.2.4 filelock_1.0.2
[5] Biostrings_2.68.0 bitops_1.0-7
[7] fastmap_1.1.1 RCurl_1.98-1.12
[9] promises_1.2.0.1 digest_0.6.31
[11] mime_0.12 lifecycle_1.0.3
[13] ellipsis_0.3.2 base64_2.0.1
[15] KEGGREST_1.40.0 interactiveDisplayBase_1.38.0
[17] RSQLite_2.3.1 magrittr_2.0.3
[19] compiler_4.3.0 rlang_1.1.1
[21] tools_4.3.0 utf8_1.2.3
[23] yaml_2.3.7 askpass_1.1
[25] S4Arrays_1.0.1 bit_4.0.5
[27] curl_5.0.0 DelayedArray_0.26.1
[29] plyr_1.8.8 RColorBrewer_1.1-3
[31] BiocParallel_1.34.0 withr_2.5.0
[33] purrr_1.0.1 grid_4.3.0
[35] stats4_4.3.0 preprocessCore_1.62.0
[37] fansi_1.0.4 wheatmap_0.2.0
[39] xtable_1.8-4 colorspace_2.1-0
[41] ggplot2_3.4.2 MASS_7.3-59
[43] scales_1.2.1 SummarizedExperiment_1.30.1
[45] cli_3.6.1 crayon_1.5.2
[47] generics_0.1.3 rstudioapi_0.14
[49] reshape2_1.4.4 httr_1.4.5
[51] tzdb_0.3.0 DBI_1.1.3
[53] cachem_1.0.7 stringr_1.5.0
[55] zlibbioc_1.46.0 parallel_4.3.0
[57] AnnotationDbi_1.62.0 XVector_0.40.0
[59] matrixStats_0.63.0 vctrs_0.6.2
[61] Matrix_1.5-4 IRanges_2.34.0
[63] hms_1.1.3 S4Vectors_0.38.0
[65] bit64_4.0.5 glue_1.6.2
[67] codetools_0.2-19 stringi_1.7.12
[69] gtable_0.3.3 BiocVersion_3.17.1
[71] later_1.3.0 GenomeInfoDb_1.36.0
[73] GenomicRanges_1.52.0 munsell_0.5.0
[75] tibble_3.2.1 pillar_1.9.0
[77] illuminaio_0.42.0 rappdirs_0.3.3
[79] htmltools_0.5.5 openssl_2.0.6
[81] GenomeInfoDbData_1.2.10 R6_2.5.1
[83] shiny_1.7.4 Biobase_2.60.0
[85] lattice_0.21-8 readr_2.1.4
[87] png_0.1-8 memoise_2.0.1
[89] httpuv_1.6.9 Rcpp_1.0.10
[91] MatrixGenerics_1.12.0 pkgconfig_2.0.3

[zwdzwd]

Glad it worked! We have a video on the previous version. https://www.youtube.com/playlist?list=PLKRu7cmBQlaiauS9JadpsacqFrI9xSXpL
But the best documentation is the software vignette.
https://bioconductor.org/packages/release/bioc/vignettes/sesame/inst/doc/sesame.html

[Jmase0921]

Thank you again! I have one more question. Does the output age shown above give the correct/est age? The output number above (20.91063) seems to be a bit off. The known age of the person idats that I used is 31 yrs old and the output from R is not that. Thank you again for your response.

[Jmase0921]

Hi @zwdzwd ,

I hope you are doing well! I wanted to follow up to my question above?

Thank you again! I have one more question. Does the output age shown above give the correct/est age? The output number above (20.91063) seems to be a bit off. The known age of the person idats that I used is 31 yrs old and the output from R is not that. Thank you again for your response.

Thank you for your help again!

[zwdzwd]

It is hard for me to tell whether this is biological or technical. Have you tested it on more public datasets?

[Jmase0921]

Hi @zwdzwd ,

Tank you for responding. No I have not tried on public data. Do you have a recommendation on where to get this public data that would have the right age? I can also look around for this as well. Thank you again!

[neverseverette]

Hi zwdzwd! I was having issues creating vcfs (Forum #93) and tried what you suggested, however. I'm still getting the same error message with verbose=TRUE and vcf = NULL. When running that for loop, it only executes a certain number of vcfs for each sample. The samples that are showing the error message are samples that passed in Sesame and in GS. I'm not sure why only a handful of samples aren't creating vcfs. Is there another way to produce vcfs for all files?

[Jmase0921]

Hey, Have you received a response?

[Jmase0921]

Hey @neverseverette ,

Did you receive a message from @zwdzwd ?

[zwdzwd]

I replied in #93 . We have a new formatVCF function that generalize the genotype output. I still haven't run into this error myself. Maybe some info on how to reproduce the error on a public dataset would be very helpful.

[slenhard]

Hi zwdzwd! Is there an updated version of SeSAMe to handle EPICv2? I'm able to generate vcfs for EPICv1, but the same process is not working for EPICv2 idats. This is the error:

Error: BiocParallel errors
1 remote errors, element index: 1
46 unevaluated and other errors
first remote error:
Error in .sesameDataGet(title, use_alternative = use_alternative): Data EPICplus.address not found.

Thank you!

[zwdzwd]

If possible, could you update your sesame to the following versions:

R: 4.3.0
Bioconductor: 3.17
sesame: 1.19.7
sesameData: 1.19.0
ExperimentHub: 2.8.1

EPICv2 should be supported in the 1.18+s

[slenhard]

Updating the versions definitely helped! Thank you!

Now I am able to create the sdf, but there is an error when trying to create the vcfs. Here is my input:

for(name in names(sdfTEST)){tryCatch({formatVCF(sdfTEST[[name]], vcf = paste("D:\Methylation\TEST_vcfs\", name, ".vcf", sep=''),genome = "hg19", annoS = NULL, annoI = NULL, verbose = TRUE)}, error=function(e){cat("ERROR :",conditionMessage(e), "\n",closeAllConnections())})}

And here is the error:

ERROR : not all 'rnames' found or unique.
adding rname 'https://github.com/zhou-lab/InfiniumAnnotationV1/raw/main/Anno/EPICv2/EPICv2.hg19.snp_overlap_b151.rds'

[zwdzwd]

I updated the Sesame vcf production to support EPICv2. Here is the new documentation (you need the latest sesame).
Briefly, the following code should work.

library(tidyverse)
anno = read_tsv(sesameAnno_download("EPICv2.hg38.snp.tsv.gz"))
vcf = formatVCF(sdf, anno)

A major update here is we now report genotype on all Infinium-I probes. So the output will be bigger. The rsID overlap is put to the INFO column. REF and ALT factor in the base ambiguity (Hs and Ds)
Let me know if this works or if you need any info.

> sesame_checkVersion()
SeSAMe requires matched versions of R, sesame, sesameData and ExperimentHub.
Here is the current versions installed:
R: 4.3.0
Bioconductor: 3.18
sesame: 1.19.10
sesameData: 1.19.0
ExperimentHub: 2.9.1

[slenhard]

This was very helpful! Thank you so much!

[slenhard]

I did need to use R version 4.3.1, but all the other versions match what you recommended on my end.

[tamefelis]

Hi sesame users and Dr. Wanding,

may i ask, if anyone had tried incorporate the sesame data into epigenetic clocks??
I'm trying to convert the sesame:SigDF annotations back to the RGsets format from minfi RGChannelSets.

I tried to look at the documentation/codes on the https://github.com/zwdzwd/sesamize/blob/main/sesamize.R , from line no. 185,

appreciate if anyone active here could shed some light on me. 👯‍♂️