## ----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()