[Bug]: cannot run aggregate_adjustment1.sh

[ks164]

Contact Details

No response

Scripts

03b-methylation_adjustment1.sh

What happened?

I ran 03b-methylation_adjustment in parallel on our cluster, generated the appropriate array of .Rdata files.

I am unable to use ./resources/methylation/aggregate_adjustement1.sh to combine the files, i get error:

Aggregating methylation_adjustment1 chunks There were 18 warnings (use warnings() to see them) Checking all files are present Error in main() : object 'norm.beta' not found Execution halted

N.B. the warnings are all related to R versions packages are installed under

How can the bug be reproduced?

No response

R version

4.3.3 (February, 2024)

Python version

None

Relevant log output

Aggregating methylation_adjustment1 chunks
There were 18 warnings (use warnings() to see them)
Checking all files are present
Error in main() : object 'norm.beta' not found
Execution halted

[ZXiaopu]

Hi,

Sorry for the late reply.

object norm.beta was supposed to be used in aggregate_adjustement2 rather than aggregate_adjustement1. one possibility is there is "pc" in your variable ${transformed_methylation_adjusted}

Can you go to your script path and then run

source config
echo ${transformed_methylation_adjusted} | grep pc

If this is the reason, can you change script [script path]/resources/methylation/aggregate_chunks.R line 30-32 to

if(grepl("methylation_adjusted_pcs", rootname)){
            betas <- norm.beta
}

please let me know whether it works.

[ks164]

thanks so much for this; I think the reason is that my path to the data contains the string "pc" since our organization's cluster naming contains "hpc".

I have changed lines 30-32 in aggregate_chunks.R to be explicit as suggested, and the code is currently running.

[ks164]

Unfortunately, the job failed with the same error as previously even after changing the code.

[ZXiaopu]

Can you show me your code in line 30-32 and also the results with below code

cd [Script path]
source config
${R_directory}Rscript ${scripts_directory}/resources/methylation/aggregate_chunks.R \
    ${transformed_methylation_adjusted} \
    ${meth_chunks} \
    ${methylation_array} \
    ${methylation_no_outliers} \
    ${subset}

[ks164]

thanks for looking into this!

This is my amended code in the file './resources/methylation/aggregate_chunks.R ', with surrounding lines of code for reference fs <- fs[index] load(fs[1])

    if(grepl("methylation_adjusted_pcs", rootname)){
            betas <- norm.beta
    }

    if(is.na(rootname1)){
        message("Reading in data")

This is the result of running the code you provided:

source config ${R_directory}Rscript ${scripts_directory}/resources/methylation/aggregate_chunks.R \ ${transformed_methylation_adjusted} \ ${meth_chunks} \ ${methylation_array} \ ${methylation_no_outliers} \ ${subset} There were 18 warnings (use warnings() to see them) Checking all files are present Reading in data Error in main() : object 'norm.beta' not found Execution halted

In case it helps at all, my slurm submission script is identical to the example 'godmc_03b_aggregate.sh' except I don't specify a partition.

I haven't progressed beyond step 03b, so I haven't generated any pcs yet, and grepping for 'pcs' instead of 'pc' doesn't return any matches.

[ZXiaopu]

Can you show me the output of

source config
echo ${scripts_directory}/resources/methylation/aggregate_chunks.R
${transformed_methylation_adjusted}
${meth_chunks}
${methylation_array}
${methylation_no_outliers}
${subset}

I would like to check the full input of the script submission which might be helpful

[ks164]

Thank you!

here's the output:

/hpc/group/sugdenlab/godmc_phase2/resources/methylation/aggregate_chunks.R -bash: /hpc/group/sugdenlab/godmc_phase2/processed_data/methylation_data/transformed_methylation_adjusted: No such file or directory -bash: 100: command not found -bash: epic: command not found -bash: /hpc/group/sugdenlab/godmc_phase2/processed_data/methylation_data/methylation_no_outliers.Robj: cannot execute binary file: Exec format error

The 'transformed_methylation_adjusted' .Rdata files do exist in /processed_data/ and I have checked them as instructed in the guidelines for running on a cluster.

[ZXiaopu]

Can you also change line 40 to 42 to

    if(grepl("methylation_adjusted_pcs", rootname)){
            betas <- norm.beta
    }

Hope it works

[ks164]

Thanks so much, this worked!