## ----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 ## ----load_data------------------------------------------------------------- getCTSS(ce) ce ## -------------------------------------------------------------------------- CTSStagCountSE(ce) CTSScoordinatesGR(ce) CTSStagCountDF(ce) ## -------------------------------------------------------------------------- head(CTSScoordinates(ce)) head(CTSStagCountDf(ce)) head(CTSStagCount(ce)) ## -------------------------------------------------------------------------- sampleLabels(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") ## -------------------------------------------------------------------------- mergeSamples(ce, mergeIndex = c(3,2,4,4,1), mergedSampleLabels = c("zf_unfertilized_egg", "zf_high", "zf_30p_dome", "zf_prim6")) annotateCTSS(ce, exampleZv9_annot) ## -------------------------------------------------------------------------- librarySizes(ce) ## ----ReverseCumulatives, fig.cap="Reverse cumulative distribution of CAGE tags"---- plotReverseCumulatives(ce, fitInRange = c(5, 1000), onePlot = TRUE) ## -------------------------------------------------------------------------- normalizeTagCount(ce, method = "powerLaw", fitInRange = c(5, 1000), alpha = 1.2, T = 5*10^4) ce[["tagCountMatrix"]] ## -------------------------------------------------------------------------- clusterCTSS( object = ce , threshold = 1 , thresholdIsTpm = TRUE , nrPassThreshold = 1 , method = "distclu" , maxDist = 20 , removeSingletons = TRUE , keepSingletonsAbove = 5) ## -------------------------------------------------------------------------- head(tagClusters(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") ## -------------------------------------------------------------------------- cumulativeCTSSdistribution(ce, clusters = "tagClusters", useMulticore = T) 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) ## -------------------------------------------------------------------------- aggregateTagClusters(ce, tpmThreshold = 5, qLow = 0.1, qUp = 0.9, maxDist = 100) ce$outOfClusters / ce$librarySizes *100 ## -------------------------------------------------------------------------- consensusClustersGR(ce) ## -------------------------------------------------------------------------- annotateConsensusClusters(ce, exampleZv9_annot) cumulativeCTSSdistribution(ce, clusters = "consensusClusters", useMulticore = TRUE) 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) ## -------------------------------------------------------------------------- exportCTSStoBedGraph(ce, values = "normalized", format = "bedGraph", oneFile = TRUE) ## -------------------------------------------------------------------------- exportCTSStoBedGraph(ce, values = "normalized", format = "BigWig") ## ----CTSSbedGraph, echo=FALSE, fig.cap="CAGE data bedGraph track visualized in the UCSC Genome Browser"---- knitr::include_graphics("images/CTSSbedGraph.svg") ## -------------------------------------------------------------------------- exportToBed(object = ce, what = "tagClusters", qLow = 0.1, qUp = 0.9, oneFile = TRUE) ## ----tagClustersBed, echo=FALSE, fig.align="center", fig.cap="Tag clusters visualization in the genome browser"---- knitr::include_graphics("images/TagClustersBed.svg") ## -------------------------------------------------------------------------- # Not supported yet for CAGEexp objects, sorry. # getExpressionProfiles(ce, what = "consensusClusters", tpmThreshold = 10, # nrPassThreshold = 1, method = "som", xDim = 4, yDim = 2) ## -------------------------------------------------------------------------- # Not supported yet for CAGEexp objects, sorry. # plotExpressionProfiles(ce, what = "consensusClusters") ## ----ConsensusClustersExpressionProfiles, echo=FALSE, fig.align="center", fig.cap="Expression clusters"---- knitr::include_graphics("images/ConsensusClustersExpressionProfiles.svg") ## -------------------------------------------------------------------------- # Not supported yet for CAGEexp objects, sorry. # class3_1 <- extractExpressionClass(ce, # what = "consensusClusters", which = "3_1") # head(class3_1) ## -------------------------------------------------------------------------- # Not supported yet for CAGEexp objects, sorry. # exportToBed(ce, what = "consensusClusters", # colorByExpressionProfile = TRUE) ## ----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) ## -------------------------------------------------------------------------- cumulativeCTSSdistribution(ce, clusters = "consensusClusters") ## -------------------------------------------------------------------------- # Not supported yet for CAGEexp objects, sorry. # 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, # 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") ## ----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()