## ----setup, include=FALSE----------------------------------------------------- knitr::opts_chunk$set(echo = TRUE) knitr::opts_chunk$set(collapse = TRUE, comment = "#>") library(devtools) load_all("./") ## ----eval = FALSE,message=FALSE----------------------------------------------- # #To install this package, start R (version "3.6" or later) and enter: # #if (!requireNamespace("BiocManager", quietly = TRUE)) # # install.packages("BiocManager") # # # #BiocManager::install("spiky") # # library(spiky) ## ----------------------------------------------------------------------------- data(spike) ## ----------------------------------------------------------------------------- sb <- system.file("extdata", "example.spike.bam", package="spiky", mustWork=TRUE) outFasta <- paste(system.file("extdata", package="spiky", mustWork=TRUE),"/spike_contigs.fa",sep="") show(generate_spike_fasta(sb, spike=spike,fa=outFasta)) ## ----------------------------------------------------------------------------- spikes <- system.file("extdata", "spikes.fa", package="spiky", mustWork=TRUE) spikemeth <- spike$methylated process_spikes(spikes, spikemeth) ## ----eval=TRUE---------------------------------------------------------------- genomic_bam_path <- system.file("extdata", "example_chr21.bam", package="spiky", mustWork=TRUE) genomic_coverage <- scan_genomic_contigs(genomic_bam_path,spike=spike) spike_bam_path <- system.file("extdata", "example.spike.bam", package="spiky", mustWork=TRUE) spikes_coverage <- scan_spike_contigs(spike_bam_path,spike=spike) ## ----------------------------------------------------------------------------- genomic_bedpe_path <- system.file("extdata", "example_chr21_bedpe.bed.gz", package="spiky", mustWork=TRUE) genomic_coverage <- scan_genomic_bedpe(genomic_bedpe_path,genome="hg38") spike_bedpe_path <- system.file("extdata", "example_spike_bedpe.bed.gz", package="spiky", mustWork=TRUE) spikes_coverage <- scan_spike_bedpe(spike_bedpe_path,spike=spike) ## ----eval=TRUE---------------------------------------------------------------- ##Calculate methylation specificity methyl_spec <- methylation_specificity(spikes_coverage,spike=spike) print(methyl_spec) ## ----eval=TRUE---------------------------------------------------------------- ## Build the Gaussian generalized linear model on the spike-in control data gaussian_glm <- model_glm_pmol(spikes_coverage,spike=spike) summary(gaussian_glm) ## ----eval=TRUE---------------------------------------------------------------- # Predict pmol concentration # To select a genome other than hg38, use BSgenome::available.packages() to find valid BSgenome name #library("BSgenome.Hsapiens.UCSC.hg38") sample_data_pmol <- predict_pmol(gaussian_glm, genomic_coverage,bsgenome="BSgenome.Hsapiens.UCSC.hg38",ret="df") head(sample_data_pmol,n=1) ## ----eval=TRUE---------------------------------------------------------------- sample_binned_data <- bin_pmol(sample_data_pmol) head(sample_binned_data,n=1) ## ----------------------------------------------------------------------------- sessionInfo()