## ----style, echo=FALSE, results='hide', message=FALSE-------------------------
library(BiocStyle)
library(knitr)
opts_chunk$set(error = FALSE, message = FALSE, warning = FALSE)
opts_chunk$set(fig.asp = 1)

## ----installation, echo=TRUE, eval=FALSE--------------------------------------
# ## try http:// if https:// URLs are not supported
# if (!requireNamespace("BiocManager", quietly=TRUE))
#     install.packages("BiocManager")
# BiocManager::install("Melissa")
# 
# ## Or download from Github repository
# # install.packages("devtools")
# devtools::install_github("andreaskapou/Melissa", build_vignettes = TRUE)

## ----bismark, eval=FALSE------------------------------------------------------
# # Requires Bismark
# bismark_methylation_extractor --comprehensive --merge_non_CpG \
#   --no_header --gzip --bedGraph input_file.bam

## ----binarise, eval=FALSE-----------------------------------------------------
# library(Melissa)
# # Binarise scBS-seq data
# binarise_files(indir = "path_to_met_files_dir")

## ----compress_files, eval=FALSE-----------------------------------------------
# gzip filenames

## ----melissa_data_obj, echo=TRUE, message=FALSE, eval=FALSE-------------------
# melissa_data <- create_melissa_data_obj(met_dir = "path_to_bin_met_dir",
#                                         anno_file = "anno_file", cov = 3)

## ----save_obj, eval=FALSE-----------------------------------------------------
# saveRDS(file = "melissa_data_obj.rds", melissa_data)

## ----filter_regions_by_coverage, eval=FALSE-----------------------------------
# melissa_data <- filter_by_cpg_coverage(melissa_data, min_cpgcov = 10)

## ----filter_regions_by_variability, eval=FALSE--------------------------------
# melissa_data <- filter_by_variability(melissa_data, min_var = 0.2)

## ----filter_by_coverage_across_cells, eval=FALSE------------------------------
# melissa_data <- filter_by_coverage_across_cells(melissa_data,
#                                                 min_cell_cov_prcg = 0.5)

## ----save_obj_filtered, eval=FALSE--------------------------------------------
# saveRDS(file = "melissa_data_obj_filtered.rds", melissa_data)

## ----eval=FALSE---------------------------------------------------------------
# #=================
# # 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

## ----eval=FALSE---------------------------------------------------------------
# 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

## ----eval=FALSE---------------------------------------------------------------
# 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

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

## ----eval=FALSE---------------------------------------------------------------
# bismark_genome_preparation hg19/

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

## ----session_info, echo=TRUE, message=FALSE-----------------------------------
sessionInfo()

