Hi,
I have been pre-processing my RRBS data with Bismark for which I have obtained bam files and coverage files that could be loaded into MethylKit for downstream analysis. However, the methRead function ran more than half a day on my HPCC cluster where I have 28 samples with coverage files as input. I did not receive any error messages. Below is my code, may I ask whether this is normal, or is my process stalled in the background without any messages popping up?
library(methylKit)
library(data.table)
coverage_files <- list.files(path = "/home/lab/bismark_methylation", pattern = "\\.cov.gz$")
coverage_files <- paste0("/home/lab/bismark_methylation/", coverage_files)
coverage_fread <- lapply(coverage_files, fread)
meta <- read_excel("/home/lab/supplementary_table.xlsx")
meta <- meta[match(substr(coverage_files , start = 57 , stop = 66), meta$Sample.Name),]
myobj <- methRead(list(coverage_fread),sample.id=meta$Sample.Name,pipeline = "bismarkCoverage", assembly="hg18", treatment = list(ifelse(meta$type%in% "disease", 1, 0)),context = "CpG")
saveRDS(myobj, file = "/home/lab/methRead.RDS")
•
written
by
CTLong •
1 answer
You are passing a list of data.table objects (loaded via fread) into methRead():
coverage_fread <- lapply(coverage_files, fread)
myobj <- methRead(list(coverage_fread), ...)
This is not expected input. The methRead() function expects a list of file paths, not data tables.
Given that you are passing a list of lists of data.tables and assuming that methRead should fail with a warning indicting that length of 'location' and 'name' should be same, I would say that your process stalls during the reading of the coverage files.
Also, please note that the function list.files has the argument full.names which returns the file path with the path prepended.
Try this cleaned-up and correct version of your script:
library(methylKit)
coverage_files <- list.files(path = "/home/lab/bismark_methylation", pattern = "\\.cov.gz$", full.names = TRUE)
meta <- read_excel("/home/lab/supplementary_table.xlsx")
meta <- meta[match(substr(coverage_files , start = 57 , stop = 66), meta$Sample.Name),]
myobj <- methRead(list(coverage_fread),sample.id=meta$Sample.Name,pipeline = "bismarkCoverage", assembly="hg18", treatment = list(ifelse(meta$type%in% "disease", 1, 0)),context = "CpG")
saveRDS(myobj, file = "/home/lab/methRead.RDS")
Log in to answer this question.
