## ----pqsfinder, message=FALSE, warning=FALSE---------------------------------- library(pqsfinder) ## ----basic_detection, results='hold'------------------------------------------ seq <- DNAString("TTTTGGGCGGGAGGAGTGGAGTTTTTAACCCCAAAAATTTGGGAGGGTGGGTGGGAGAA") pqs <- pqsfinder(seq, min_score = 20) pqs ## ----accessors---------------------------------------------------------------- elementMetadata(pqs) ## ----overlapping-------------------------------------------------------------- pqsfinder(seq, overlapping = TRUE, min_score = 30) ## ----density------------------------------------------------------------------ pqs <- pqsfinder(seq, deep = TRUE, min_score = 20) density(pqs) ## ----density_viz, fig.height=1.5, fig.width=8.0, message=FALSE---------------- library(Gviz) ss <- DNAStringSet(seq) names(ss) <- "chr1" dtrack <- DataTrack( start = 1:length(density(pqs)), width = 1, data = density(pqs), chromosome = "chr1", genome = "", name = "density") strack <- SequenceTrack(ss, chromosome = "chr1", name = "sequence") suppressWarnings(plotTracks(c(dtrack, strack), type = "h")) ## ----perfect_g4, results='hold'----------------------------------------------- pqsfinder(seq, max_defects = 0, min_score = 20) ## ----min_score, results='hold'------------------------------------------------ pqsfinder(seq, min_score = 70) ## ----granges_conversion------------------------------------------------------- gr <- as(pqs, "GRanges") gr ## ----granges_export----------------------------------------------------------- library(rtracklayer) export(gr, "test.gff", version = "3") ## ----gff_file, echo=FALSE, comment=NA----------------------------------------- text <- readLines("test.gff") cat(strwrap(text, width = 80, exdent = 3), sep = "\n") ## ----dnastringset_conversion-------------------------------------------------- dss <- as(pqs, "DNAStringSet") dss ## ----dnastringset_export------------------------------------------------------ writeXStringSet(dss, file = "test.fa", format = "fasta") ## ----fasta_file, echo=FALSE, comment=NA--------------------------------------- text <- readLines("test.fa") cat(text, sep = "\n") ## ----rwe_libraries, message=FALSE--------------------------------------------- library(pqsfinder) library(BSgenome.Hsapiens.UCSC.hg38) library(rtracklayer) library(ggplot2) library(Gviz) ## ----rwe_biomart, warning=FALSE----------------------------------------------- gnm <- "hg38" gene <- "AHNAK" # Load cached AHNAK region track: load(system.file("extdata", "gtrack_ahnak.RData", package="pqsfinder")) # Alternatively, query biomaRt API with the following commands: # library(biomaRt) # gtrack <- BiomartGeneRegionTrack(genome = gnm, symbol = gene, name = gene) ## ----rwe_seq------------------------------------------------------------------ extend <- 1000 seq_start <- min(start(gtrack)) - extend seq_end <- max(end(gtrack)) + extend chr <- chromosome(gtrack) seq <- Hsapiens[[chr]][seq_start:seq_end] ## ----rwe_pqsfinder, results='hide'-------------------------------------------- pqs <- pqsfinder(seq, deep = TRUE) ## ----rwe_show, message=FALSE-------------------------------------------------- pqs ## ----rwe_sort----------------------------------------------------------------- pqs_s <- pqs[order(score(pqs), decreasing = TRUE)] pqs_s ## ----rwe_gff------------------------------------------------------------------ export(as(pqs, "GRanges"), "test.gff", version = "3") ## ----rwe_gff_file, echo=FALSE, comment=NA------------------------------------- text <- readLines("test.gff", n = 5) cat(strwrap(text, width = 80, exdent = 3), sep= "\n") ## ----rwe_fasta---------------------------------------------------------------- writeXStringSet(as(pqs, "DNAStringSet"), file = "test.fa", format = "fasta") ## ----rwe_fasta_file, echo=FALSE, comment=NA----------------------------------- text <- readLines("test.fa", n = 6) cat(text, sep = "\n") ## ----rwe_hist, fig.height=4, fig.width=6-------------------------------------- sf <- data.frame(score = score(pqs)) ggplot(sf) + geom_histogram(mapping = aes(x = score), binwidth = 5) ## ----rwe_viz, fig.height=4.0, fig.width=8.0----------------------------------- strack <- DataTrack( start = start(pqs)+seq_start, end = end(pqs)+seq_start, data = score(pqs), chromosome = chr, genome = gnm, name = "score") dtrack <- DataTrack( start = (seq_start):(seq_start+length(density(pqs))-1), width = 1, data = density(pqs), chromosome = chr, genome = gnm, name = "density") atrack <- GenomeAxisTrack() suppressWarnings(plotTracks(c(gtrack, strack, dtrack, atrack), type = "h")) ## ----custom_scoring_fn-------------------------------------------------------- c_loop_bonus <- function(subject, score, start, width, loop_1, run_2, loop_2, run_3, loop_3, run_4) { l1 <- run_2 - loop_1 l2 <- run_3 - loop_2 l3 <- run_4 - loop_3 if (l1 == l2 && l1 == l3 && subject[loop_1] == DNAString("C") && subject[loop_1] == subject[loop_2] && subject[loop_1] == subject[loop_3]) { score <- score + 20 } return(score) } ## ----no_custom_scoring, results='hold'---------------------------------------- seq <- DNAString("GGGCGGGCGGGCGGGAAAAAAAAAAAAAGGGAGGGAGGGAGGG") pqsfinder(seq) ## ----custom_scoring, results='hold'------------------------------------------- pqsfinder(seq, custom_scoring_fn = c_loop_bonus) ## ----isg4_scoring_function---------------------------------------------------- isG4 <- function(subject, score, start, width, loop_1, run_2, loop_2, run_3, loop_3, run_4) { r1 <- loop_1 - start r2 <- loop_2 - run_2 r3 <- loop_3 - run_3 r4 <- start + width - run_4 if (!(r1 == r2 && r1 == r3 && r1 == r4)) return(0) run_1_s <- subject[start:start+r1-1] run_2_s <- subject[run_2:run_2+r2-1] run_3_s <- subject[run_3:run_3+r3-1] run_4_s <- subject[run_4:run_4+r4-1] if (length(grep("^G+$", run_1_s)) && length(grep("^C+$", run_2_s)) && length(grep("^G+$", run_3_s)) && length(grep("^C+$", run_4_s))) return(r1 * 20) else return(0) } ## ----isg4_pqsfinder, results='hold'------------------------------------------- pqsfinder(DNAString("AAAAGGGATCCCTAAGGGGTCCC"), strand = "+", use_default_scoring = FALSE, run_re = "G{3,6}|C{3,6}", custom_scoring_fn = isG4) ## ----sessionInfo, echo=FALSE-------------------------------------------------- sessionInfo()