02-Import.Rmddata_dir = "~/Documents/GSE90553/" dir.create(path = data_dir, showWarnings = FALSE, recursive = TRUE) filePaths <- GEOquery::getGEOSuppFiles(c("GSM3039352", "GSM3039353", "GSM3039354", "GSM3039348", "GSM3039349", "GSM3039350"), baseDir = data_dir, makeDirectory = TRUE)
Above command should download files to data_dir.
bed_files <- list.files( path = getwd(), pattern = "*bed\\.gz$", full.names = TRUE, recursive = TRUE ) print(basename(bed_files))
As a first step, we need a list of CpG sites in the respective reference genome. The CpG sites are extracted using the respective Bsgenome annotation package. The read_bedgraph function is also able to extract CpG sites on it’s own, however, it might be beneficial to do it separately.
hg19_cpgs <- methrix::extract_CPGs(ref_genome = "BSgenome.Hsapiens.UCSC.hg19")
An annotation table is also necessary to perform analyses. The data will be added to the methrix object, as a colData slot.
sample_anno <- data.frame( row.names = gsub( pattern = paste0("GSM[[:digit:]]*_BiSeq_(WT_D[[:digit:]]{1,2}_R[[:digit:]]).bed.gz"), replacement = "\\1", x = basename(bed_files) ), Day = gsub( pattern = paste0("GSM[[:digit:]]*_BiSeq_WT_(D[[:digit:]]{1,2})_R[[:digit:]].bed.gz"), replacement = "\\1", x = basename(bed_files) ), Replicate = gsub( pattern = paste0("GSM[[:digit:]]*_BiSeq_WT_D[[:digit:]]{1,2}_(R[[:digit:]]).bed.gz"), replacement = "\\1", x = basename(bed_files) ), stringsAsFactors = FALSE ) print(sample_anno)
read_bedgraphs function is a versatile bedgraph reader intended to import bedgraph files generated virtually by any sort of methylation calling program. It requires user to provide indices for chromosome names, start position and other required fields. However, there are also presets available to import bedgraphs from most common programs such as
Bismark (.cov format)MethylDackelMethylcToolsBisSNPBSseeker2_CGmapIn case if the bedGraphs are generated using any of the above programs use pipiline argument in read_bedgraphs() to specify whihc causes the function to automatically infer the requored column indices.
For the example in hand, since we do not know the program used to generate the files, it might worth to take a look at one of the files, and manually specify the column indices.
data.table::fread(bed_files[1], nrows = 5)
The read_bedgraphs function adds CpGs missing from the reference set, and creates a methylation/coverage matrices. Once the process is complete - it returns an object of class methrix which in turn inherits SummarizedExperiment class. methrix object contains ‘methylation’ and ‘coverage’ matrices (either in-memory or as on-disk HDF5 arrays) along with pheno-data and other basic info. This object can be passed to all downstream functions for various analysis. For further details on the data structure, see the SummarizedExperiment package.
For computers with large memory, increase speed using vect=TRUE. The function will open multiple files, defined by the vect_batch_size argument. Keep in mind, that increasing the batch size will increase the memory need as well.
If vect=FALSE, the function will iterate along the list of files, keeping only one open at a time.
A more detailed description of the read_bedgraphs function is here: read_bedgraphs)
meth <- methrix::read_bedgraphs( files = bed_files, ref_cpgs = hg19_cpgs, chr_idx = 1, start_idx = 2, M_idx = 5, U_idx = 6, stranded = FALSE, zero_based = TRUE, collapse_strands = FALSE, coldata = sample_anno, vect = TRUE, n_threads = 3 )
Too many missing CpGs, if not expected, might indicate that something went wrong.
An obvious error is using the wrong reference genome.
A common mistake is falsely indicating strandedness and zero or one based positioning. The combination of the two might lead to only using C-s from one strand (due to a 1 bp shift), resulting lower coverage for the sites. A nice description of the zero- and one-based coordination system can be found here