## ----setup, echo = FALSE, results = "hide"------------------------------------ options(signif = 3, digits = 3) knitr::opts_chunk$set(tidy = FALSE, cache = TRUE, autodep = TRUE, fig.height = 5.5, message = FALSE, error = FALSE, warning = TRUE) set.seed(0xdada) ## ----CAGEprotocol, echo=FALSE, fig.align="center", fig.cap="Overview of CAGE experiment"---- knitr::include_graphics("images/CAGEprotocol.svg") ## ----load_CAGEr--------------------------------------------------------------- library(CAGEr) ## ----specify_input_files------------------------------------------------------ inputFiles <- list.files( system.file("extdata", package = "CAGEr") , "ctss$" , full.names = TRUE) ## ----create_CAGEexp----------------------------------------------------------- ce <- CAGEexp( genomeName = "BSgenome.Drerio.UCSC.danRer7" , inputFiles = inputFiles , inputFilesType = "ctss" , sampleLabels = sub( ".chr17.ctss", "", basename(inputFiles)) ) ## ----display_created_object--------------------------------------------------- ce ## ----show_colData------------------------------------------------------------- colData(ce) ## ----load_data---------------------------------------------------------------- ce <- getCTSS(ce) ce ## ----------------------------------------------------------------------------- CTSStagCountSE(ce) CTSScoordinatesGR(ce) CTSStagCountDF(ce) CTSStagCountGR(ce, 1) # GRanges for one sample with expression count. ## ----------------------------------------------------------------------------- sampleLabels(ce) ## ----------------------------------------------------------------------------- ce <- annotateCTSS(ce, exampleZv9_annot) ## ----------------------------------------------------------------------------- colData(ce)[,c("librarySizes", "promoter", "exon", "intron", "unknown")] CTSScoordinatesGR(ce) ## ----------------------------------------------------------------------------- plotAnnot(ce, "counts") ## ----CorrelationScatterPlots, fig.cap="Correlation of raw CAGE tag counts per TSS"---- corr.m <- plotCorrelation2( ce, samples = "all" , tagCountThreshold = 1, applyThresholdBoth = FALSE , method = "pearson") ## ----------------------------------------------------------------------------- ce <- mergeSamples(ce, mergeIndex = c(3,2,4,4,1), mergedSampleLabels = c("Zf.unfertilized.egg", "Zf.high", "Zf.30p.dome", "Zf.prim6")) ce <- annotateCTSS(ce, exampleZv9_annot) ## ----------------------------------------------------------------------------- librarySizes(ce) ## ----ReverseCumulatives, fig.cap="Reverse cumulative distribution of CAGE tags"---- plotReverseCumulatives(ce, fitInRange = c(5, 1000), onePlot = TRUE) ## ----------------------------------------------------------------------------- ce <- normalizeTagCount(ce, method = "powerLaw", fitInRange = c(5, 1000), alpha = 1.2, T = 5*10^4) ce[["tagCountMatrix"]] ## ----------------------------------------------------------------------------- ce <- clusterCTSS( ce , threshold = 1 , thresholdIsTpm = TRUE , nrPassThreshold = 1 , method = "distclu" , maxDist = 20 , removeSingletons = TRUE , keepSingletonsAbove = 5) ## ----------------------------------------------------------------------------- tagClustersGR(ce, sample = "Zf.unfertilized.egg") ## ----CumulativeDistribution, echo=FALSE, fig.align="center", fig.cap="Cumulative distribution of CAGE signal and definition of interquantile width"---- knitr::include_graphics("images/CumulativeDistributionAndQuantiles.svg") ## ----------------------------------------------------------------------------- ce <- cumulativeCTSSdistribution(ce, clusters = "tagClusters", useMulticore = T) ce <- quantilePositions(ce, clusters = "tagClusters", qLow = 0.1, qUp = 0.9) ## ----------------------------------------------------------------------------- tagClustersGR( ce, "Zf.unfertilized.egg" , returnInterquantileWidth = TRUE, qLow = 0.1, qUp = 0.9) ## ----------------------------------------------------------------------------- plotInterquantileWidth(ce, clusters = "tagClusters", tpmThreshold = 3, qLow = 0.1, qUp = 0.9) ## ----------------------------------------------------------------------------- ce <- aggregateTagClusters(ce, tpmThreshold = 5, qLow = 0.1, qUp = 0.9, maxDist = 100) ce$outOfClusters / ce$librarySizes *100 ## ----------------------------------------------------------------------------- consensusClustersGR(ce) ## ----------------------------------------------------------------------------- ce <- annotateConsensusClusters(ce, exampleZv9_annot) ce <- cumulativeCTSSdistribution(ce, clusters = "consensusClusters", useMulticore = TRUE) ce <- quantilePositions(ce, clusters = "consensusClusters", qLow = 0.1, qUp = 0.9, useMulticore = TRUE) ## ----------------------------------------------------------------------------- consensusClustersGR( ce, sample = "Zf.unfertilized.egg" , returnInterquantileWidth = TRUE, qLow = 0.1, qUp = 0.9) ## ----------------------------------------------------------------------------- trk <- exportToTrack(CTSSnormalizedTpmGR(ce, "Zf.30p.dome")) ce |> CTSSnormalizedTpmGR("all") |> exportToTrack(ce, oneTrack = FALSE) ## ----------------------------------------------------------------------------- split(trk, strand(trk), drop = TRUE) ## ----CTSSbedGraph, echo=FALSE, fig.cap="CAGE data bedGraph track visualized in the UCSC Genome Browser"---- knitr::include_graphics("images/CTSSbedGraph.svg") ## ----------------------------------------------------------------------------- iqtrack <- exportToTrack(ce, what = "tagClusters", qLow = 0.1, qUp = 0.9, oneTrack = FALSE) iqtrack #rtracklayer::export.bed(iqtrack, "outputFileName.bed") ## ----tagClustersBed, echo=FALSE, fig.align="center", fig.cap="Tag clusters visualization in the genome browser"---- knitr::include_graphics("images/TagClustersBed.svg") ## ----------------------------------------------------------------------------- ce <- getExpressionProfiles(ce, what = "consensusClusters", tpmThreshold = 10, nrPassThreshold = 1, method = "som", xDim = 4, yDim = 2) consensusClustersGR(ce)$exprClass |> table(useNA = "always") ## ----ConsensusClustersExpressionProfiles-------------------------------------- plotExpressionProfiles(ce, what = "consensusClusters") ## ----ConsensusClustersExpressionProfiles_svg, echo=FALSE, fig.align="center", fig.cap="Expression clusters"---- knitr::include_graphics("images/ConsensusClustersExpressionProfiles.svg") ## ----------------------------------------------------------------------------- consensusClustersGR(ce) |> subset(consensusClustersGR(ce)$exprClass == "0_1") ## ----------------------------------------------------------------------------- cc_iqtrack <- exportToTrack(ce, what = "consensusClusters", colorByExpressionProfile = TRUE) cc_iqtrack #rtracklayer::export.bed(cc_iqtrack, "outputFileName.bed") ## ----ConsensusClustersBed, echo=FALSE, fig.align="center", fig.cap="Consensus clusters colored by expression profile in the genome browser"---- knitr::include_graphics("images/ConsensusClustersBed.svg") ## ----------------------------------------------------------------------------- ce$group <- factor(c("a", "a", "b", "b")) dds <- consensusClustersDESeq2(ce, ~group) ## ----------------------------------------------------------------------------- ce <- cumulativeCTSSdistribution(ce, clusters = "consensusClusters") ## ----------------------------------------------------------------------------- ce <- scoreShift(ce, groupX = "Zf.unfertilized.egg", groupY = "Zf.prim6", testKS = TRUE, useTpmKS = FALSE) ## ----ShiftingScore, echo=FALSE, fig.cap="Calculation of shifting score"------- knitr::include_graphics("images/ShiftingScore.svg") ## ----------------------------------------------------------------------------- # Not supported yet for CAGEexp objects, sorry. shifting.promoters <- getShiftingPromoters(ce, groupX = "Zf.unfertilized.egg", groupY = "Zf.prim6", tpmThreshold = 5, scoreThreshold = 0.6, fdrThreshold = 0.01) head(shifting.promoters) ## ----ShiftingPromoter, echo=FALSE, fig.cap="Example of shifting promoter"----- knitr::include_graphics("images/ShiftingPromoter.svg") ## ----------------------------------------------------------------------------- ce <- quickEnhancers(ce) ce[["enhancers"]] ## ----create_df---------------------------------------------------------------- TSS.df <- read.table(system.file( "extdata/Zf.unfertilized.egg.chr17.ctss" , package = "CAGEr")) # make sure the column names are as required colnames(TSS.df) <- c("chr", "pos", "strand", "Zf.unfertilized.egg") # make sure the column classes are as required TSS.df$chr <- as.character(TSS.df$chr) TSS.df$pos <- as.integer(TSS.df$pos) TSS.df$strand <- as.character(TSS.df$strand) TSS.df$Zf.unfertilized.egg <- as.integer(TSS.df$Zf.unfertilized.egg) head(TSS.df) ## ----coerce_to_df------------------------------------------------------------- ce.coerced <- as(TSS.df, "CAGEexp") ce.coerced ## ----sessionInfo-------------------------------------------------------------- sessionInfo()