methylKit script to run on klone

[yaaminiv]

@sr320 Here's my methylKit script to run on klone: https://github.com/RobertsLab/project-oyster-oa/blob/master/code/Haws/04-methylKit-parameter-testing.Rmd

Raw data here: https://gannet.fish.washington.edu/spartina/project-oyster-oa/Haws/bismark-2/

Let me know if there are any errors/formatting that needs to be modified

[sr320]

It might be me... but that Rmd seemed fraught. - feels like might have worked with other files in environment and believe pH is not coded correctly..
I did of series of edits at https://github.com/RobertsLab/project-oyster-oa/blob/master/code/Haws/04.2-methylKit-parameter-testing.Rmd

Ran but to the extent I trust is limited.

I did pull off from scratch - https://gannet.fish.washington.edu/v1_web/temp/haws-bs/analysis_report.html
Figures are not rendering but available at https://gannet.fish.washington.edu/v1_web/temp/haws-bs/

[yaaminiv]

feels like might have worked with other files in environment and believe pH is not coded correctly..
I did of series of edits at https://github.com/RobertsLab/project-oyster-oa/blob/master/code/Haws/04.2-methylKit-parameter-testing.Rmd

Can you explain what you mean by this?

[sr320]

Main aspect was chunks would often not run because object being calling was not created prior in Rmd.

…

On Tue, Jul 14, 2026 at 4:24 PM Dr. Yaamini R. Venkataraman < ***@***.***> wrote: *yaaminiv* left a comment (RobertsLab/resources#2472) <https://urldefense.com/v3/__https://github.com/RobertsLab/resources/issues/2472*issuecomment-4974984030__;Iw!!K-Hz7m0Vt54!mN-Nda8ck6-9BLxC6RaN4iNU1L7i2do64IuP-T_Hp0EhGnNJMAFc8-qrgIwS0e_DG8ZYJfUAk9Lo8u5cVEgpyWo$> feels like might have worked with other files in environment and believe pH is not coded correctly.. I did of series of edits at https://github.com/RobertsLab/project-oyster-oa/blob/master/code/Haws/04.2-methylKit-parameter-testing.Rmd <https://urldefense.com/v3/__https://github.com/RobertsLab/project-oyster-oa/blob/master/code/Haws/04.2-methylKit-parameter-testing.Rmd__;!!K-Hz7m0Vt54!mN-Nda8ck6-9BLxC6RaN4iNU1L7i2do64IuP-T_Hp0EhGnNJMAFc8-qrgIwS0e_DG8ZYJfUAk9Lo8u5cM5AtyaE$> Can you explain what you mean by this? — Reply to this email directly, view it on GitHub <https://urldefense.com/v3/__https://github.com/RobertsLab/resources/issues/2472?email_source=notifications&email_token=ABB4PN2DXVXH3I2FYJBXJPT5E26LLA5CNFSNUABFM5UWIORPF5TWS5BNNB2WEL2JONZXKZKDN5WW2ZLOOQXTIOJXGQ4TQNBQGMYKM4TFMFZW63VHNVSW45DJN5XKKZLWMVXHJLDGN5XXIZLSL5RWY2LDNM*issuecomment-4974984030__;Iw!!K-Hz7m0Vt54!mN-Nda8ck6-9BLxC6RaN4iNU1L7i2do64IuP-T_Hp0EhGnNJMAFc8-qrgIwS0e_DG8ZYJfUAk9Lo8u5chTBpJgo$>, or unsubscribe <https://urldefense.com/v3/__https://github.com/notifications/unsubscribe-auth/ABB4PN7VE6GURGE4TSFAZ5T5E26LLAVCNFSNUABEKJSXA33TNF2G64TZHM4TSNRZHA3DQOB3JFZXG5LFHM2DQNJZGEZTSOJSHCQXMAQ__;!!K-Hz7m0Vt54!mN-Nda8ck6-9BLxC6RaN4iNU1L7i2do64IuP-T_Hp0EhGnNJMAFc8-qrgIwS0e_DG8ZYJfUAk9Lo8u5cwZNGZdw$> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[kubu4]

believe pH is not coded correctly..

I think this is the correct assessment.

When the metadata file is loaded and manipulated (https://github.com/RobertsLab/project-oyster-oa/blob/c9cf2718a6c5ae3e05ba7304201ea2dd08101825/code/Haws/04-methylKit-parameter-testing.Rmd#L93-L98), it's looking for a column named pH, but there is no column pH. The column name is actually ph.

  dplyr::rename(pH = ph) %>%
  mutate(ploidyTreatment = case_when(ploidy == "3N" ~ 1,
                                     ploidy == "2N" ~ 0)) %>%
  mutate(pHTreatment = case_when(pH == "high" ~ 0,
                                 pH == "low" ~ 0)) #Import metadata and modify for use with methylKit
head(sampleMetadata) #Confirm dataframe creation

Steven's edits also modified the binary assignment for pH in this step, since this original code sets both "high" and "low" conditions to the same binary value of 0.

[kubu4]

To clarify, HOWEVER, Steven's edits did not fix the pH label issue, Here's his code (https://github.com/RobertsLab/project-oyster-oa/blob/c9cf2718a6c5ae3e05ba7304201ea2dd08101825/code/Haws/04.2-methylKit-parameter-testing.Rmd#L93-L98):

  dplyr::rename(pH = ph) %>%
  mutate(ploidyTreatment = case_when(ploidy == "3N" ~ 1,
                                     ploidy == "2N" ~ 0)) %>%
  mutate(pHTreatment = case_when(pH == "high" ~ 0,
                                 pH == "low" ~ 1)) #Import metadata and modify for use with methylKit
head(sampleMetadata) #Confirm dataframe creation

[yaaminiv]

When the metadata file is loaded and manipulated (https://github.com/RobertsLab/project-oyster-oa/blob/c9cf2718a6c5ae3e05ba7304201ea2dd08101825/code/Haws/04-methylKit-parameter-testing.Rmd#L93-L98), it's looking for a column named pH, but there is no column pH. The column name is actually ph.

Right, which is why I use dplyr::rename(pH = ph) to rename the column in the sampleMetadata dataframe.

Steven's edits also modified the binary assignment for pH in this step, since this original code sets both "high" and "low" conditions to the same binary value of 0.

YIKES I fixed that in another script but I guess not this one.

Main aspect was chunks would often not run because object being calling was
not created prior in Rmd.

Could you specify which ones? I checked this prior to sending but I clearly missed some!

[kubu4]

Right, which is why I use dplyr::rename(pH = ph) to rename the column in the sampleMetadata dataframe.

Except, the code is looking for pH, not ph like the column is named.

mutate(pHTreatment = case_when(pH == "high" ~ 0,

[yaaminiv]

I did pull off from scratch - https://gannet.fish.washington.edu/v1_web/temp/haws-bs/analysis_report.html

How the interaction was done is very interesting...may consider doing that myself but also not quite sold on the validity of the analytical approach. Need to look into that further.

Ran but to the extent I trust is limited.

Once the Rmd is fixed, would be great to run it since there are some parameters I am testing that I would like to evaluate prior to finalizing DML identification.

[yaaminiv]

Except, the code is looking for pH, not ph like the column is named.

No! Once ph is renamed to pH, it can find the pH column and assign information in pHtreatment accordingly.

[kubu4]

Ah! Duh! I was thinking the renaming was the other way around (pH -> ph). Then, the subsequent calls to pH would fail.

NEVERMIND!

[sr320]

Is there a new script ready to roll?

[yaaminiv]

I've fixed the pH issue in the original script: https://github.com/RobertsLab/project-oyster-oa/blob/master/code/Haws/04-methylKit-parameter-testing.Rmd

Main aspect was chunks would often not run because object being calling was
not created prior in Rmd.

Do you know specifically which objects were causing issues? I've combed through it again but am not sure what specifically was the issue the first time

[sr320]

First fail - line 85

Error in read.csv("../../data/Haws/sample_metadata.csv") %>% dplyr::select(-c(1)) %>% : could not find function "%>%"s

yet

session info confirms tidy verse?

R version 4.4.3 (2025-02-28)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 24.04.4 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8       
 [4] LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
[10] LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: Etc/UTC
tzcode source: system (glibc)

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

other attached packages:
 [1] pheatmap_1.0.12      RColorBrewer_1.1-3   scales_1.4.0         vegan_2.7-5         
 [5] permute_0.9-10       methylKit_1.32.1     GenomicRanges_1.58.0 GenomeInfoDb_1.42.3 
 [9] IRanges_2.40.1       S4Vectors_0.44.0     BiocGenerics_0.52.0  BiocManager_1.30.25 
[13] data.table_1.18.4    ggplot2_4.0.3        tidyverse_2.0.0     

loaded via a namespace (and not attached):
 [1] farver_2.1.2                R.utils_2.13.0              Biostrings_2.74.1          
 [4] S7_0.2.2                    bitops_1.0-9                fastmap_1.2.0              
 [7] RCurl_1.98-1.19             GenomicAlignments_1.42.0    XML_3.99-0.18              
[10] digest_0.6.39               lifecycle_1.0.5             cluster_2.1.8.1            
[13] statmod_1.5.2               magrittr_2.0.5              compiler_4.4.3             
[16] rlang_1.2.0                 tools_4.4.3                 yaml_2.3.12                
[19] rtracklayer_1.66.0          knitr_1.51                  S4Arrays_1.6.0             
[22] mclust_6.1.3                curl_7.1.0                  DelayedArray_0.32.0        
[25] plyr_1.8.9                  abind_1.4-8                 BiocParallel_1.40.2        
[28] withr_3.0.2                 numDeriv_2016.8-1.1         R.oo_1.27.1                
[31] grid_4.4.3                  gtools_3.9.5                MASS_7.3-65                
[34] SummarizedExperiment_1.36.0 bbmle_1.0.25.1              cli_3.6.6                  
[37] mvtnorm_1.4-1               rmarkdown_2.31              crayon_1.5.3               
[40] otel_0.2.0                  rstudioapi_0.19.0           httr_1.4.8                 
[43] reshape2_1.4.4              rjson_0.2.23                bdsmatrix_1.3-7            
[46] qvalue_2.38.0               stringr_1.6.0               zlibbioc_1.52.0            
[49] splines_4.4.3               parallel_4.4.3              XVector_0.46.0             
[52] restfulr_0.0.17             matrixStats_1.5.0           vctrs_0.6.5                
[55] Matrix_1.7-3                jsonlite_2.0.0              fastseg_1.52.0             
[58] limma_3.62.2                glue_1.8.1                  emdbook_1.3.14             
[61] codetools_0.2-20            stringi_1.8.7               gtable_0.3.6               
[64] BiocIO_1.16.0               UCSC.utils_1.2.0            tibble_3.3.1               
[67] pillar_1.11.1               htmltools_0.5.9             GenomeInfoDbData_1.2.13    
[70] R6_2.6.1                    evaluate_1.0.5              lattice_0.22-7             
[73] Biobase_2.66.0              R.methodsS3_1.8.2           Rsamtools_2.22.0           
[76] Rcpp_1.0.14                 nlme_3.1-168                SparseArray_1.6.2          
[79] coda_0.19-4.1               mgcv_1.9-3                  xfun_0.58                  
[82] MatrixGenerics_1.18.1       pkgconfig_2.0.3        

[yaaminiv]

huh...I don't know what to do about that. Wonder if loading dplyr specifically would help?