3 A workflow for PAR-CLIP data analysis
3.1 Reading sorted BAM files
A sorted and indexed BAM file can be loaded into the R session with readSortedBam. This function extracts aligned reads, sequences and the mismatch MD field, and returns a GRanges object.
library(wavClusteR)## Loading required package: GenomicRanges## Loading required package: stats4## Loading required package: BiocGenerics## Loading required package: parallel##
## Attaching package: 'BiocGenerics'## The following objects are masked from 'package:parallel':
##
## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
## clusterExport, clusterMap, parApply, parCapply, parLapply,
## parLapplyLB, parRapply, parSapply, parSapplyLB## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs## The following objects are masked from 'package:base':
##
## Filter, Find, Map, Position, Reduce, anyDuplicated, append,
## as.data.frame, cbind, colMeans, colSums, colnames, do.call,
## duplicated, eval, evalq, get, grep, grepl, intersect,
## is.unsorted, lapply, lengths, mapply, match, mget, order,
## paste, pmax, pmax.int, pmin, pmin.int, rank, rbind, rowMeans,
## rowSums, rownames, sapply, setdiff, sort, table, tapply,
## union, unique, unsplit, which, which.max, which.min## Loading required package: S4Vectors##
## Attaching package: 'S4Vectors'## The following object is masked from 'package:base':
##
## expand.grid## Loading required package: IRanges## Loading required package: GenomeInfoDb## Loading required package: Rsamtools## Loading required package: Biostrings## Loading required package: XVector##
## Attaching package: 'Biostrings'## The following object is masked from 'package:base':
##
## strsplitfilename <- system.file( "extdata", "example.bam", package = "wavClusteR" )
Bam <- readSortedBam(filename = filename)
Bam## GRanges object with 5358 ranges and 2 metadata columns:
## seqnames ranges strand | qseq
## <Rle> <IRanges> <Rle> | <DNAStringSet>
## [1] chrX [24001819, 24001844] - | CAGAGATAAA...TATTTTAAAG
## [2] chrX [24001819, 24001843] - | CAGAGATAAA...ATATTTTAAG
## [3] chrX [24001834, 24001863] - | ATATTTTAGA...ATTTTATTTA
## [4] chrX [24001836, 24001865] - | TTTTTAAAGA...TTTATTTAAA
## [5] chrX [24001841, 24001876] - | AAAGATTAAA...TTTTCTTCAT
## ... ... ... ... . ...
## [5354] chrX [24023018, 24023051] - | GTTTCACAGC...AAAAATATGT
## [5355] chrX [24023018, 24023051] - | GTTTCACAGC...AAAAATATGT
## [5356] chrX [24023019, 24023051] - | TTTCACAGCG...AAAAATATGT
## [5357] chrX [24023019, 24023051] - | TTTCACAGCG...AAAAATATGT
## [5358] chrX [24023067, 24023090] - | CAAAGGCGCG...GGTTTATTTT
## MD
## <character>
## [1] 26
## [2] 24A0
## [3] 8A21
## [4] 0A13A15
## [5] 24A11
## ... ...
## [5354] 10A23
## [5355] 10A23
## [5356] 9A23
## [5357] 9A23
## [5358] 9A14
## -------
## seqinfo: 25 sequences from an unspecified genome; no seqlengths3.2 Extracting informative genomic positions
Prior to model fitting, genome-wide substitutions need to be identified and filtered based on a coverage threshold. The getAllSub function identifies all genomic positions exhibiting at least one substitution and coverage above this threshold.
countTable <- getAllSub( Bam, minCov = 10 )## Loading required package: doMC## Considering substitutions, n = 497, processing in 1 chunks## Warning: executing %dopar% sequentially: no parallel backend registered## chunk #: 1## considering the + strand## Computing local coverage at substitutions...## considering the - strand## Computing local coverage at substitutions...head( countTable )## GRanges object with 6 ranges and 3 metadata columns:
## seqnames ranges strand | substitutions coverage
## <Rle> <IRanges> <Rle> | <character> <numeric>
## [1] chrX [24001959, 24001959] - | TC 17
## [2] chrX [24001973, 24001973] - | TC 17
## [3] chrX [24001977, 24001977] - | TC 13
## [4] chrX [24002046, 24002046] - | TC 10
## [5] chrX [24002057, 24002057] - | TC 10
## [6] chrX [24002147, 24002147] - | TC 22
## count
## <integer>
## [1] 2
## [2] 12
## [3] 1
## [4] 1
## [5] 6
## [6] 3
## -------
## seqinfo: 1 sequence from an unspecified genome; no seqlengthsThe function returns a GRanges object specifying genomic position, strand, substitution type (e.g. “TC”: T in the reference genome; C in the read), strand-specific coverage and number of observed substitutions at the position.
3.2.1 How to choose the coverage threshold?
The coverage threshold minCov defines the genomic positions to be used for parameter estimation, thus providing a means to tune the stringency of the analysis. Currently, wavClusteR does not allow to learn this threshold from the data. However, since the model is based on relative substitution frequencies, the value of minCov will influence the variance of the estimated parameters. Therefore, values smaller than default (minCov=10) are not recommended.
3.3 Inspecting the substitutions profile (diagnostic I)
Once all substitutions are computed, a diagnostic substitution profile can be plotted with plotSubstitutions.
plotSubstitutions( countTable, highlight = "TC" )
The function returns a barplot showing the total number of genomic positions that exhibit a given type of substitution and highlights the substitution type that is expected to be generated by the experimental procedure. The percentage of substitution of this type is also shown. This plot can be readily used to assess the quality of the sequenced libraries and can be used to compare different data sets generated under the same experimental condition.
3.4 Estimating the model
wavClusteR uses the identified genomic positions to learn a non-parametric mixture model discriminating PAR-CLIP-specific from extrinsic transitions. Model parameters are estimated by fitMixtureModel.
model <- fitMixtureModel(countTable, substitution = "TC")The function returns a list of:
- the two mixing coefficients (l1 and l2)
- the two model components (p1 and p2)
- the full density (p)
As the small size of our example data set does not to estimate the model reliably, the mixture model for the entire AGO2 dataset has been precomputed and is provided by the package.
data(model)
str(model)## List of 5
## $ l1: Named num 0.181
## ..- attr(*, "names")= chr "TC"
## $ l2: Named num 0.819
## ..- attr(*, "names")= chr "TC"
## $ p : num [1:999] 7.52 9.44 10.05 10.38 10.48 ...
## $ p1: num [1:999] 89.6 64.4 50.4 41.5 35.3 ...
## $ p2: num [1:999] 0 0 1.14 3.51 5 ...3.5 Inspecting the model fit (diagnostic II)
The model fit can be inspected using getExpInterval.
(support <- getExpInterval( model, bayes = TRUE ) )
## $supportStart
## [1] 0.007
##
## $supportEnd
## [1] 0.98The function returns two diagnostic plots. The first plot illustrates the estimated densities \(p\), \(p_1\) and \(p_2\), and log odds
\[ o=]
The second plot shows the resulting posterior class probability, i.e. the probability that a given relative substitution frequency (RSF, horizontal axis) has been induced by PAR-CLIP. The area under the curve corresponding to the returned RSF interval is colored, and the RSF interval indicated. By default, getExpInterval returns the RSF interval according to the Bayes classifier, i.e. a posterior probability cutoff of 0.5. However, the user can specify a custom RSF interval:
By means of the rightProb and leftProb parameters, e.g.
(support <- getExpInterval( model, bayes = FALSE, leftProb = 0.9, rightProb = 0.9 ) )
## $supportStart ## [1] 0.076 ## ## $supportEnd ## [1] 0.905By inspecting the posterior class probability density and passing the RSF interval boundaries when calling high-confidence substitutions (see point 6).
Finally, the model can be used to produce further diagnostic plots. Particularly, the total number of reads exhibitng a given substitution and an RSF-based partitioning of genomic positions with substitutions can be obtained by
plotSubstitutions( countTable, highlight = "TC", model )
3.6 Selecting high-confidence PAR-CLIP induced transitions
The RSF support is used to filter all observed transitions to define a set of high-confidence, PAR-CLIP induced transitions. These are identified by getHighConfSub. The function returns a GRanges object with genomic position, strand, strand-specific coverage, occurence (count), and relative substitution frequency (rsf) for each high-confidence substitution. In a call to getHighConfSub, the RSF interval returned by getExpInterval can be supplied as support argument directly
highConfSub <- getHighConfSub( countTable,
support = support,
substitution = "TC" )
head( highConfSub ) ## GRanges object with 6 ranges and 3 metadata columns:
## seqnames ranges strand | coverage count
## <Rle> <IRanges> <Rle> | <numeric> <integer>
## [1] chrX [24001959, 24001959] - | 17 2
## [2] chrX [24001973, 24001973] - | 17 12
## [3] chrX [24001977, 24001977] - | 13 1
## [4] chrX [24002046, 24002046] - | 10 1
## [5] chrX [24002057, 24002057] - | 10 6
## [6] chrX [24002147, 24002147] - | 22 3
## rsf
## <numeric>
## [1] 0.117647058823529
## [2] 0.705882352941177
## [3] 0.0769230769230769
## [4] 0.1
## [5] 0.6
## [6] 0.136363636363636
## -------
## seqinfo: 1 sequence from an unspecified genome; no seqlengthsor, alternatively, by specifying supportStart and supportEnd, which define the range of RSF of interest.
highConfSub <- getHighConfSub( countTable,
supportStart = 0.2,
supportEnd = 0.7,
substitution = "TC" )
head( highConfSub )3.7 Identifying protein binding sites (clusters)
Binding sites (clusters) are identified by getClusters. The function takes high-confidence substitution sites and the coverage function
coverage <- coverage( Bam )
coverage$chrX## integer-Rle of length 24023090 with 914 runs
## Lengths: 24001818 15 2 ... 1 15 24
## Values : 0 2 3 ... 21 0 1as an input. From package version 2.0, cluster boundaries are resolved using the Mini-Rank Norm (MRN) Comoglio et al., 2015, which is up to 10x faster and more sensitive than the previously adopted algorithm based on continuous wavelet transform of the coverage function Sievers et al., 2012. Briefly, the MRN algorithm finds an optimal cluster boundary for each high-confidence substitution by solving an optimization problem that integrates prior knowledge on the geometry of PAR-CLIP clusters. Two options are available:
Hard thresholding, i.e. the coverage function is denoised using a global threshold. Empirically,
minCov=1worked well on all tested datasets for whichminCov = 10. Alternatively, 10% of the mode of the coverage distribution at high-confidence substitutions represents a possible choice.clusters <- getClusters( highConfSub = highConfSub, coverage = coverage, sortedBam = Bam, method = "mrn", threshold = 1, cores = 1 )## Computing start/end read positions## Number of chromosomes exhibiting high confidence transitions: 1## Loading required package: doMC## ...Processing = chrXclusters## GRanges object with 184 ranges and 0 metadata columns: ## seqnames ranges strand ## <Rle> <IRanges> <Rle> ## [1] chrX [24001953, 24001975] - ## [2] chrX [24001953, 24001976] - ## [3] chrX [24001953, 24001977] - ## [4] chrX [24002044, 24002057] - ## [5] chrX [24002044, 24002057] - ## ... ... ... ... ## [180] chrX [24006559, 24006573] - ## [181] chrX [24006560, 24006565] - ## [182] chrX [24007061, 24007068] - ## [183] chrX [24007061, 24007083] - ## [184] chrX [24007061, 24007083] - ## ------- ## seqinfo: 1 sequence from an unspecified genome; no seqlengthsLocal thresholding, based on a global estimation of background levels via a Gaussian mixture model. Omitting the threshold parameter in the call to getClusters enables local thresholding
clusters <- getClusters( highConfSub = highConfSub, coverage = coverage, sortedBam = Bam, method = "mrn", cores = 1 )## Computing start/end read positions## Learning background threshold by fitting a GMM## Estimated threshold (% of maximum local coverage differences) from 184 sampled transitions: 5.26## Number of chromosomes exhibiting high confidence transitions: 1## Loading required package: doMC## ...Processing = chrXclusters## GRanges object with 184 ranges and 0 metadata columns: ## seqnames ranges strand ## <Rle> <IRanges> <Rle> ## [1] chrX [24001953, 24001959] - ## [2] chrX [24001963, 24001973] - ## [3] chrX [24001972, 24001977] - ## [4] chrX [24002044, 24002057] - ## [5] chrX [24002046, 24002057] - ## ... ... ... ... ## [180] chrX [24006569, 24006574] - ## [181] chrX [24006569, 24006574] - ## [182] chrX [24007062, 24007068] - ## [183] chrX [24007073, 24007076] - ## [184] chrX [24007074, 24007077] - ## ------- ## seqinfo: 1 sequence from an unspecified genome; no seqlengths
Note: the CWT-based method is maintained for reproducibility and can be called by specifying method = "cwt" in the function call above.
3.8 Merging clusters
Once the clusters are identified, the corresponding genomic regions can be merged in a strand-specific manner. Statistics for each merged cluster (a wavCluster), can be computed using filterClusters.
require(BSgenome.Hsapiens.UCSC.hg19)## Loading required package: BSgenome.Hsapiens.UCSC.hg19## Loading required package: BSgenome## Loading required package: rtracklayerwavclusters <- filterClusters( clusters = clusters,
highConfSub = highConfSub,
coverage = coverage,
model = model,
genome = Hsapiens,
refBase = "T",
minWidth = 12)## Computing log odds...## Refining cluster sizes...## Combining clusters...## Quantifying transitions within clusters...## Computing statistics...##
|
|= | 2%
|
|=== | 5%
|
|==== | 7%
|
|====== | 9%
|
|======= | 11%
|
|========= | 14%
|
|========== | 16%
|
|============ | 18%
|
|============= | 20%
|
|=============== | 23%
|
|================ | 25%
|
|================== | 27%
|
|=================== | 30%
|
|===================== | 32%
|
|====================== | 34%
|
|======================== | 36%
|
|========================= | 39%
|
|=========================== | 41%
|
|============================ | 43%
|
|============================== | 45%
|
|=============================== | 48%
|
|================================ | 50%
|
|================================== | 52%
|
|=================================== | 55%
|
|===================================== | 57%
|
|====================================== | 59%
|
|======================================== | 61%
|
|========================================= | 64%
|
|=========================================== | 66%
|
|============================================ | 68%
|
|============================================== | 70%
|
|=============================================== | 73%
|
|================================================= | 75%
|
|================================================== | 77%
|
|==================================================== | 80%
|
|===================================================== | 82%
|
|======================================================= | 84%
|
|======================================================== | 86%
|
|========================================================== | 89%
|
|=========================================================== | 91%
|
|============================================================= | 93%
|
|============================================================== | 95%
|
|================================================================ | 98%
|
|=================================================================| 100%## Consolidating results...wavclusters## GRanges object with 44 ranges and 7 metadata columns:
## seqnames ranges strand | Ntransitions MeanCov
## <Rle> <IRanges> <Rle> | <integer> <numeric>
## [1] chrX [24001950, 24001980] - | 3 14.870968
## [2] chrX [24002044, 24002057] - | 2 9.642857
## [3] chrX [24002139, 24002161] - | 3 21.956522
## [4] chrX [24002323, 24002352] - | 4 9.766667
## [5] chrX [24002670, 24002687] - | 3 14.333333
## ... ... ... ... . ... ...
## [40] chrX [24005918, 24005971] - | 5 15.14815
## [41] chrX [24006171, 24006191] - | 3 11.19048
## [42] chrX [24006537, 24006554] - | 3 32.00000
## [43] chrX [24006557, 24006577] - | 3 26.80952
## [44] chrX [24007059, 24007081] - | 3 16.30435
## NbasesInRef CrossLinkEff
## <integer> <numeric>
## [1] 12 0.25
## [2] 4 0.50
## [3] 8 0.38
## [4] 11 0.36
## [5] 7 0.43
## ... ... ...
## [40] 15 0.33
## [41] 7 0.43
## [42] 5 0.60
## [43] 7 0.43
## [44] 7 0.43
## Sequence SumLogOdds
## <factor> <numeric>
## [1] CGGCTTGGGAAAAATTACCTGACAAAAATGT 7.142001
## [2] CTAGGATTATTTGA 4.972714
## [3] TTAGGTAGAAAATAACCTCTCCC 7.457388
## [4] AATATTGAAGTTATACGGTGTACTGAAAGG 10.534644
## [5] TTAAATTATGAATTCTCA 7.795383
## ... ... ...
## [40] CAATGTTAGACCAATGGCTTTGATAGTAAAACTGGTGAGGTTTTTCTTTATATG 12.922715
## [41] AGGATGGAATCGCTGTAATGA 7.468965
## [42] GTGGAAGATGAGGTGATT 7.632600
## [43] ACCGGTGATAAACAATTTGTT 7.457649
## [44] TGCTGGTGAACATTCTGAAAGTA 8.298195
## RelLogOdds
## <numeric>
## [1] 0.5951667
## [2] 1.2431784
## [3] 0.9321735
## [4] 0.9576949
## [5] 1.1136261
## ... ...
## [40] 0.8615143
## [41] 1.0669950
## [42] 1.5265200
## [43] 1.0653784
## [44] 1.1854564
## -------
## seqinfo: 1 sequence from an unspecified genome; no seqlengthsThe function takes as input the following elements:
- clusters
- high-confidence substitution sites
- genome-wide coverage
- mixture model
- a BSgenome object containing the correct reference sequence
- reference base expected to be converted by the experimental procedure
- minimum required width of a wavCluster
and it returns a GRanges object with the following metadata:
- number of high-confidence transitions (Ntransitions)
- mean coverage (MeanCov)
- number of bases in the reference genome of the same type as the specified
refBase(NbasesInRef) - the estimated cross-linking efficiency (CrossLinkEff), i.e. the ratio between Ntransitions and NbasesInRef
- the genomic sequence (Sequence)
- the sum of the log odds (SumLogOdds), contributed by each high-confidence transition within the cluster
- the relative log odds (RelLogOdds), i.e. the ratio between SumLogOdds and Ntransitions
The relative log odds can be used to rank clusters according to statistical confidence.
3.9 Post-processing of the binding sites
wavClusteR provides a number of functions (summarized in the table below) for post-processing of the identified binding sites.
| Task | Function | Output format |
|---|---|---|
| Export all identified substitutions or high-confidence substitutions | exportHighConfSub | BED |
| Export clusters | exportClusters | BED |
| Export coverage function | exportCoverage | BigWig |
| Visualize the size distribution of wavClusters | plotSizeDistribution | histogram |
| Annotate clusters with respect to genomic features (e.g. CDS, introns, 3’-UTRs, 5’-UTRs) in a strand-specific manner | annotateClusters | dot chart, vector |
| Compute metagene profiles of wavClusters, where the density of wavClusters is represented as a function of a reference genomic coordinates | getMetaGene | line plot, vector |
| Compute metaTSS profiles based on all aligned reads in the input BAM file | line plot, vector | |
| Visualize wavClusteR statistics and meta data to learn pairwise relationships between variables | plotStatistics | pairs plot |
3.9.1 Exporting substitutions, wavClusters and coverage function
High-confidence substitutions can be exported with
exportHighConfSub( highConfSub = highConfSub,
filename = "hcTC.bed",
trackname = "hcTC",
description = "hcTC" )where trackname and description set the corresponding attributes in the UCSC BED file. By replacing highConfSub with another set of substitutions (e.g. all identified substitutions of a given type), those can be exported likewise.
wavClusters can be exported with
exportClusters( clusters = wavclusters,
filename = "wavClusters.bed",
trackname = "wavClusters",
description = "wavClusters" )and the coverage function can be exported with
exportCoverage( coverage = coverage, filename = "coverage.bigWig" )3.9.2 Annotating binding sites
wavClusters can be annotated with respect to genomic features using annotateClusters. This function generates a strand-specific dot chart representing wavClusters annotation. annotateClusters takes as an input the wavClusters and a transcriptDB object containing transcript annotations. The latter can be generated using makeTxDbFromUCSC (GenomicFeatures package)
txDB <- makeTxDbFromUCSC(genome = "hg19", tablename = "ensGene")and is automatically downloaded by annotateClusters if missing.
Then, the annotateClusters can be called as follows
annotateClusters( clusters = wavclusters,
txDB = txDB,
plot = TRUE,
verbose = TRUE)Four dot charts are returned by the function showing the percentage of clusters mapping to different transcript features localized on the same (sense) or on the opposite (antisense) strand, the relative sequence length of different compartments relative to the total transcriptome length and the normalized enrichment of clusters across functional compartments.
Note: multiple hits, i.e. wavClusters that overlap with more than one genomic feature, are reported as “multiple”, whereas wavClusters that map outside of the considered features are labeled as “other”. The latter are then annotated with respect to features on the antisense strand.
3.9.3 Computing cluster metagene profiles
A graphical representation of the density of wavClusters as a function of a binning of genomic coordinates across all annotated genes can be obtained using the getMetaGene function.
getMetaGene( clusters = wavclusters,
txDB = txDB,
upstream = 1e3,
downstream = 1e3,
nBins = 40,
nBinsUD = 10,
minLength = 1,
plot = TRUE,
verbose = TRUE ) In this example, genes were divided in 40 bins (nBins) and an upstream/downstream region spanning 1kb was considered. This region was subdivided in 10 bins (nBinsUD). No restriction on gene length was applied (minLength). getMetaGene returns a numeric vector of length nBins + 2*nBinsUD with normalized counts, which can be used, for instance, to compare the distribution of wavClusters across several PAR-CLIP samples.
In addition to metagene profiles, meta transcription start site (TSS) profiles based on all mapped reads can be generated using getMetaTSS.
getMetaTSS( sortedBam = Bam,
txDB = txDB,
upstream = 1e3,
downstream = 1e3,
nBins = 40,
unique = FALSE,
plot = TRUE,
verbose = TRUE ) Here, the upstream and downstream parameters control the width of the window to be considered, and nBins controls the resolution of the profile. If unique=TRUE, then overlapping windows are discarded. getMetaTSS returns a numeric vector of length nBins with normalized read counts.
3.9.4 Computing the size distribution of wavClusters
The size distribution of wavClusters can be visualized by
plotSizeDistribution( clusters = wavclusters, showCov = TRUE, col = "skyblue2" )## `geom_smooth()` using method = 'loess'
3.9.5 Visualizing wavCluster statistics and meta data
Cluster statistics can be plotted as pairs plot using
plotStatistics( clusters = wavclusters,
corMethod = "spearman",
lower = panel.smooth )