Repliscope is an R package for creating, normalising, comparing and plotting DNA replication timing profiles. The analysis pipeline starts with BED-formatted read count files (output of localMapper) obtained by high-throughput sequencing of DNA from replicating and non-replicating cells. There are three methods of measuring DNA replication dynamics using relative copy number (Fig.1): sort-seq, sync-seq and marker frequency analysis (MFA-seq). Sort-seq uses fluorescence-activated cell sorting (FACS) to enrich for non-replicating and replicating cells from an asynchronous population. Sync-seq requires cells to be arrested in non-replicating cell cycle phase (i.e. G1), followed by release into S phase. Samples are then taken throughout S phase when cells synchronously synthesise DNA according to the replication timing programme. In the case of MFA-seq, rapidly dividing cells in exponential growth phase are directly used as the replicating sample, while a saturated culture serves as a non-replicating control sample. While the latter approach of obtaining cells is the simplest, it also requires deeper sequencing due to decreased dynamic range and, thus, is more suitable for organisms with small genomes (typically, bacteria).
<img src=“flavours.png” alt=“Fig.1: Replication time profiling family”, width = “95%”>
For best experience, use the Repliscope in interactive mode. To do so, simply run the runGUI() function.
The typical command line analysis using Repliscope starts with loading BED-formatted read count files using the loadBed function. Various functions allow removal of genomic bins containing low quality data (rmChr,rmOutliers). To aid read count analysis, two visualisation functions are used: plotBed and plotCoverage. Next, read-depth adjusted ratio of reads from the replicating sample to reads in the non-replicating sample is calculated using the makeRatio function; the resulting ratio values are distributed around one. The normaliseRatio function is then used to transpose the ratio values to a biologically relevant relative DNA copy number scale (typically, from 1 to 2). The normalised ratio values are, essentially, replication timing profiles, where low values indicate late replication and high values - early replication. The plotRatio function helps to visualise the ratio values as a histogram. Genomic bins containing unusually high or low ratio values may be removed using the trimRatio function. smoothRatio uses cubic spline to smooth replication profile data. compareRatios can be used to calculate difference between two replication profiles using z-score statistics. Finally, replication profiles are plotted using the plotGenome function, which also allows for various genome annotations.
The runGUI function launches the package in interactive mode as a Shiny app. This is the recommended mode for the majority of applications.
runGUI()runGUI()The loadBed function reads BED-formatted files into an R session. The BED format specifications can be found at UCSC. When read using this function, The localMapper output BED files typically look like this:
| chrom | chromStart | chromEnd | name | score |
|---|---|---|---|---|
| chrI | 0 | 1000 | W303_G2 | 481 |
| chrI | 1000 | 2000 | W303_G2 | 311 |
| chrI | 2000 | 3000 | W303_G2 | 471 |
| chrI | 3000 | 4000 | W303_G2 | 1298 |
| chrI | 4000 | 5000 | W303_G2 | 1124 |
The BED file fields may be separated by tabs or spaces. If the BED file contains a header, it will be ignored. If a genome dataframe is provided, only data for the chromosomes/regions within it will be retained (this may be used to split genomes containing many chromosomes into several plots). The resulting data is ordered by “chromStart” columns, retaining original chromosome order, unless a genome dataframe is provided.
loadBed(file,genome,name)W303_G2 <- loadBed(system.file("extdata/W303_G2.bed",package="Repliscope"), name='W303_G2')
W303_G2_chrI <- loadBed(system.file("extdata/W303_G2.bed",package="Repliscope"), name='W303_G2',genome=sacCer3[["genome"]])A function to remove single chromosome data from a bed dataframe
rmChr(bed,chr)W303_S <- rmChr(W303_S,"chrM")A function to remove outliers from the “score” column of a supplied bed dataframe. There are three methods available: max, IQR and median. Max is used to remove 1 or more maximum values; IQR uses the InterQuartile Range to detect outliers, while the median method can be used to remove data based on the genome-wide median.
rmOutliers(bed,method,n=1,range=3,loLim=0.25,hiLim=NULL)bedDF <- rmOutliers(W303_S,method="max",n=2)
bedDF <- rmOutliers(W303_S,method="IQR",range=3)
bedDF <- rmOutliers(W303_S,method="median",loLim=0.25,hiLim=2)A function to plot values in the ‘score’ column of the supplied bed dataframe, a boxplot per unique chromosome name in the ‘chrom’ column. The resulting plot also highlights outliers based on the interquartile range. The genome wide median is plotted as a pink line through the boxplots.
plotBed(W303_G2)<img src=“plotBed.png” alt=“Fig.2: Output of the plotBed function”, width = “95%”>
plotBed(bed,plotting=TRUE)plotBed(W303_S)
plotObject <- plotBed(W303_G2,plotting=FALSE)A function to create a scatterplot of values in the ‘score’ column of the supplied bed dataframe, along chromosome coordinates. The genome wide median is plotted as a pink line.
plotCoverage(W303_S,region='chrII:0-813136')<img src=“plotCoverage.png” alt=“Fig.3: Output of the plotCoverage function”, width = “95%”>
plotCoverage(bed,region=FALSE,plotting=TRUE)plotCoverage(W303_G2)
plotObject <- plotCoverage(W303_S,plotting=FALSE)A function to calculate the ‘score’ ratio between two bed dataframes. It merges two supplied bed dataframes, calculates the ratios of their “score” values and normalises the ratios by the ‘score’ total sums.
| chrom | chromStart | chromEnd | name.rep | name.nonRep | ratio | ratioFactor |
|---|---|---|---|---|---|---|
| chrI | 0 | 1000 | W303_S | W303_G2 | 0.7682 | 1.0 |
| chrI | 1000 | 2000 | W303_S | W303_G2 | 0.6678 | 1.0 |
| chrI | 2000 | 3000 | W303_S | W303_G2 | 0.7277 | 1.0 |
| chrI | 3000 | 4000 | W303_S | W303_G2 | 0.7411 | 1.0 |
| chrI | 4000 | 5000 | W303_S | W303_G2 | 0.6994 | 1.0 |
makeRatio(bedRep,bedNonRep)Both bed dataframes must contain “chrom”,“chromStart”, “chromEnd” and “score” columns.
ratioDF <- makeRatio(W303_S,W303_G2)A function to normalise ratio values to fit a biologically relevant scale. It scales values from the ‘ratio’ column of the provided dataframe either using a supplied ‘rFactor’ value or automatically to best fit the scale from one to two.
normaliseRatio(ratioDF,rFactor=NULL,replace=TRUE)ratioDF <- normaliseRatio(W303)
ratioDF <- normaliseRatio(W303,rFactor=1.41,replace=FALSE)A function that plots a histogram of values in a supplied vector using ggplot2 and highlights the interval between 1 and 2 in green.
plotRatio(W303norm$ratio)<img src=“plotRatio.png” alt=“Fig.4: Output of the plotRatio function”, width = “95%”> ### Usage
plotRatio(ratio,plotting=TRUE)plotRatio(W303$ratio)
plotObject <- plotRatio(W303$ratio,plotting=FALSE)A function to remove outliers from the “ratio” column of a supplied ratio dataframe.
trimRatio(ratioDF,loLim,hiLim)W303 <- trimRatio(W303,0.5,1.5)A function to smooth ratio values using a cubic smoothing spline. The smoothRatio function splits values from the ‘ratio’ column by chromosome and based the supplied groupMin and split parameters and then applies the smooth.spline function from R stats package. The supplied dataframe may contain multiple ratios, i.e. ratios produced using multiple replicating samples and/or multiple non-replicating samples. This must be reflected in the ‘name.rep’ and ‘name.nonRep’ columns. In other words, different ratio dataframes may be merged using the rbind function before calling the smoothRatio function.
smoothRatio(ratioDF,groupMin=5,split=5)ratioDF <- smoothRatio(W303norm)
ratioDFs <- smoothRatios(rbind(Dbf4myc,W303norm))A function to compare two replication profiles. compareRatios takes two ratio dataframes that were binned the same way and uses z-score statistics to find p-values of their differences. The function outputs a combined dataframe containing the two input ratio dataframes in a long format with an added ‘p.value’ column.
compareRatios(ratio1,ratio2)ratioDFs <- compareRatios(W303norm,Dbf4myc)plotGenome plots scatterplot/barplot/polygon of the ‘score’ and/or the ‘splineSmooth’ column values by genomic coordinates, either for the whole genome or a specified region. It also allows annotation with additional BED-formatted data.
plotGenome(sortSeq,region='chrXII:0-410000',lines=sacCer3[["cen"]],circles=sacCer3[["ori"]])<img src=“plotGenome.png” alt=“Fig.5: Output of the plotGenome function”, width = “95%”>
plotGenome(ratioDFs, geom="geom_point", ylims=c(1,2), plotting=TRUE,
genome=NULL, region=FALSE, guide=NULL,
lines=NULL, circles=NULL, rectangles=NULL, pointers=NULL,
colourLines='#00FF00', colourCircles='#FFFFFF', colourRectangles='#FF0000', colourPointers='#FF7F00'
)plotGenome(sortSeq,genome=sacCer3[["genome"]],lines=sacCer3[["cen"]],circles=sacCer3[["ori"]])
plotGenome(sortSeq, geom="geom_ribbon", guide=guide,region="chrIX:250000-439885", lines=sacCer3[["cen"]], circles=sacCer3[["ori"]])
plotGenome(syncSeq[["data"]], geom="geom_segment", guide=syncSeq[["guide"]], genome=sacCer3[["genome"]], lines=sacCer3[["cen"]], circles=sacCer3[["ori"]], colourLines="black")