Convert Bismark coverage format

Melissa depends heavily on the BPRMeth package [2, 3] for reading and processing bisulfite sequencing data. It assumes that the data are first processed using Bismark [4], hence from fastq and BAM files we will obtain a coverage file by running the bismark_methylation_extractor command as shown below,

# Requires Bismark
bismark_methylation_extractor --comprehensive --merge_non_CpG \
  --no_header --gzip --bedGraph input_file.bam

The format of the coverage file is the following

<chr> <start> <end> <met_prcg> <met_reads> <unmet_reads>

where each row corresponds to an observed CpG (i.e. we have at least one read mapped to this location). Note that CpGs with no coverage are not included in this file. This format however contains redundant information, hence we bring the scBS-seq files in the format that Melissa (and BPRMeth) require, which is

<chr> <start> <met_level>

where met_level corresponds to the binary methylation state, either 0 or 1. We can do this by calling the binarise_files helper function, which only requires the input directory of the files and optionally the path to the output directory. Each file of the indir corresponds to a different cell and it can also be compressed, e.g. in .gz file format, to save storage space.

library(Melissa)
# Binarise scBS-seq data
binarise_files(indir = "path_to_met_files_dir")

Note that the new binarised files will not be compressed, due to system commands not being portable inside the R package. To save storage space, the user could compress the files for instance in .gz file format by running in command line :

gzip filenames

Create methylation regions

Now we are ready to process the binarised input files and create methylation regions using functions from the BPRMeth package. Briefly, the steps required to create this object are as follows.

1. First we require annotation data using the read_anno file. Note that the annotation file can contain any genomic context: from promoters and gene bodies to enhancers, Nanog regulatory regions and CTCF regions; hence Melissa can be used for a plethora of analyses that want to take spatial genomic correlations into account.
2. Next we need to read the methylation data using the read_met function. We will do this independently per cell.
3. Finally, the create_region_object will create the methylation regions object which is the main object for storing methylation data.

The create_melissa_data_obj is a wrapper function for doing all the above steps at once. Note that this step is important so read carefully the purpose of all the parameters to obtain the right object for downstream analysis.

melissa_data <- create_melissa_data_obj(met_dir = "path_to_bin_met_dir",
                                        anno_file = "anno_file", cov = 3)

The melissa_data$met contains the methylation data whose structure is a list of length N (number of cells), and each element of this list is another list of length M (number of genomic regions). Each entry in the inner list is an I × 2 matrix, where I are the number of CpGs, where the 1st column contains the (relative) CpG locations and the 2nd column contains the methylation state: methylated or unmethylated.

Store object

It is often useful to save this object to file with the saveRDS function. The object can then be restored using the readRDS function. This allows us to conduct downstream analysis without having to repeat the processing steps described above.

saveRDS(file = "melissa_data_obj.rds", melissa_data)

Filter by CpG coverage

Genomic regions with really sparse coverage of CpGs are not informative to infer methylation profiles. Hence, we only consider genomic regions with at least min_cpgcov CpG coverage in each region. Note that this step will not actually remove any genomic regions, it will only set to NA those regions that have coverage below the threshold.

melissa_data <- filter_by_cpg_coverage(melissa_data, min_cpgcov = 10)

Filter by mean methylation variability across cells

Genomic regions that have not heterogeneity across different cells are often of no interest, for example if we were to use them for identifying cell subpopulations. This way we will both keep only the informative genomic regions and reduce the number of genomic regions for downstream analysis for efficiency.

melissa_data <- filter_by_variability(melissa_data, min_var = 0.2)

Filter by genomic coverage across cells

Genomic regions that have coverage only on a handful of cells are not powerful for sharing information across cells. For example, a specific promoter that has observations in 5 out of the 100 cells, will not contain enough to perform sharing of information, either for imputation or clustering. Hence, regions that are are not covered in at least min_cell_cov_prcg of the cells are filtered out.

melissa_data <- filter_by_coverage_across_cells(melissa_data, 
                                                min_cell_cov_prcg = 0.5)

Store object

saveRDS(file = "melissa_data_obj_filtered.rds", melissa_data)

Smallwood et al. (2014)

The Smallwood et al. (2014) [5] dataset can be downloaded with accession number GSE56879. For this dataset we used the already processed coverage files from Bismark. The filtering of cells that do not pass quality control (QC) was done according to the original studies. See supplementary information of [5] for IDs of cells that passed filtering.

Angermueller et al. (2016)

The Angermueller et al. (2016) [6] dataset can be downloaded with accession number GSE74535. For this dataset we used the already processed coverage files from Bismark. The filtering of cells that do not pass quality control (QC) was done according to the original studies. See supplementary information of [5] for IDs of cells that passed filtering.

Bulk WGBS Encode

The analysis on the subsampled ENCODE WGBS data was performed on the bulk GM12878 (GEO GSE86765) and H1-hESC (GEO GSE80911) cell lines. The BAM files for these studies can be obtained directly from the ENCODE project portal.

#=================
# 1. Download BAM data
DATA_DIR="../encode/wgbs/"
# Download GM12878 cell line
wget -P ${DATA_DIR}GM12878/ https://www.encodeproject.org/files/ENCFF681ASN/@@download/ENCFF681ASN.bam
# Download H1-hESC cell line
wget -P ${DATA_DIR}H1hESC/ https://www.encodeproject.org/files/ENCFF546TLK/@@download/ENCFF546TLK.bam

Then we subsample the WGBS data from the BAM file, that is, we are going to remove individual reads instead of individual CpGs to take into account the nature of missing values of scBS-seq data. To do so we will run the samtools view command which subsamples random lines from a BAM file. This way, we can generate artificially 40 pseudo-single cells by keeping only 0.5% of the bulk reads for each single cell.

data_dir="encode/wgbs/GM12878/SRR4235788.bam"
out_dir="encode/wgbs/GM12878/subsampled/GM12878"
for (( i=1; i <= 40; ++i ))
do
    my_command="samtools view -s ${i}.005 -b $data_dir > ${out_dir}_${i}.bam"
    eval $my_command
done

Finally, we run the bismark_methylation_extractor command to obtain the methylation state of each covered CpG fomr the resulting BAM files. The following command will result in files of coverage output and bedGraph output.

data_dir="encode/wgbs/GM12878/subsampled/"
proc_dir="encode/wgbs/GM12878/processed/"
for (( i=1; i <= 40; ++i ))
do
    my_command="bismark_methylation_extractor --ignore 2 --comprehensive --merge_non_CpG --no_header --multicore 4 -o $proc_dir --gzip --bedGraph ${data_dir}GM12878_${i}.bam"
    eval $my_command
done

Bulk RRBS Encode

The analysis on the subsampled ENCODE RRBS data was performed again on the bulk GM12878 H1-hESC cells lines. We can download the the raw fastq files from.

http://genome.ucsc.edu/cgi-bin/hgFileUi?db=hg19&g=wgEncodeHaibMethylRrbs

and search for GM12878 or H1-hESC and download the fastq files only for the 2nd replicate.

Next we run Bismark. First run the bismark_genome_preparation command to create a genome indexing for hg19.

bismark_genome_preparation hg19/ 

After that we run the bismark command which will create alignment files in bam format.

#=================
# 3. Run bismark
bismark --genome hg19/ encode/wgEncodeHaibMethylRrbsGm12878HaibRawDataRep2.fastq.gz
bismark --genome hg19/ encode/wgEncodeHaibMethylRrbsH1hescHaibRawDataRep2.fastq.gz

After this step, we follow the same process as we did for the bulk ENCODE WGBS data above.