2.1 Input files required

Before you start, you will need:

1. Five restriction map information files (“design files”):

• Restriction map file (.rmap) - a bed file containing coordinates of the restriction fragments. By default, 4 columns: chr, start, end, fragmentID.
• Bait map file (.baitmap) - a bed file containing coordinates of the baited restriction fragments, and their associated annotations. By default, 5 columns: chr, start, end, fragmentID, baitAnnotation. The regions specified in this file, including their fragmentIDs, must be an exact subset of those in the .rmap file. The baitAnnotation is a text field that is used only to annotate the output and plots.
• nperbin file (.npb), nbaitsperbin file (.nbpb), proxOE file (.poe) - Precompute these tables from the .rmap and .baitmap files, using the Python script makeDesignFiles.py from chicagoTools at our Bitbucket repository. Refer to the chicagoTools README file for more details.

We recommend that you put all five of these files into the same directory (that we refer to as designDir). An example of a valid design folder, for a two-chromosome sample of the GM12878 data used in this vignette, is provided in the PCHiCdata package, as follows.

dataPath <- system.file("extdata", package="PCHiCdata")
testDesignDir <- file.path(dataPath, "hg19TestDesign")
dir(testDesignDir)

## [1] "h19_chr20and21.baitmap" "h19_chr20and21.nbpb"   
## [3] "h19_chr20and21.npb"     "h19_chr20and21.poe"    
## [5] "h19_chr20and21.rmap"

2. You will also need input data files, which should be in CHiCAGO input format, .chinput. You can obtain .chinput files from aligned Capture Hi-C BAM files by running bam2chicago.sh, available as part of chicagoTools. (To obtain BAM files from raw fastq files, use a Hi-C alignment & QC pipeline such as HiCUP.

Example .chinput files are provided in the PCHiCdata package, as follows:

testDataPath <- file.path(dataPath, "GMchinputFiles")
dir(testDataPath)

## [1] "GM_rep1.chinput" "GM_rep2.chinput" "GM_rep3.chinput"

files <- c(
    file.path(testDataPath, "GM_rep1.chinput"),
    file.path(testDataPath, "GM_rep2.chinput"),
    file.path(testDataPath, "GM_rep3.chinput")
  )

OPTIONAL: The data set in this vignette requires some additional custom settings, both to ensure that the vignette compiles in a reasonable time and to deal with the artificially reduced size of the data set:

settingsFile <- file.path(system.file("extdata", package="PCHiCdata"),
                          "sGM12878Settings", "sGM12878.settingsFile")

Omit this step unless you know which settings you wish to alter. If you are using non-human samples, or a very unusual cell type, one set of options that you might want to change is the weighting parameters: see Using different weights.

2.2 Example workflow

We run CHiCAGO on the test data as follows. First, we create a blank chicagoData object, and we tell it where the design files are. For this example, we also provide the optional settings file - in your analysis, you can omit the settingsFile argument.

library(Chicago)

cd <- setExperiment(designDir = testDesignDir, settingsFile = settingsFile)

The properties of chicagoData objects are discussed more in The chicagoData object.

Next, we read in the input data files:

cd <- readAndMerge(files=files, cd=cd)

## 
Read 0.0% of 295740 rows
Read 295740 rows and 5 (of 5) columns from 0.008 GB file in 00:00:04
## 
Read 0.0% of 1233748 rows
Read 1233748 rows and 5 (of 5) columns from 0.034 GB file in 00:00:08
## 
Read 0.0% of 873077 rows
Read 873077 rows and 5 (of 5) columns from 0.024 GB file in 00:00:05

Finally, we run the pipeline with chicagoPipeline():

cd <- chicagoPipeline(cd)

2.3 Output plots

chicagoPipeline() produces a number of plots. You can save these to disk by setting the outprefix argument in chicagoPipeline().

The plots are as follows (an explanation follows):

## 
Read 0.0% of 25794 rows
Read 25794 rows and 76 (of 76) columns from 0.004 GB file in 00:00:03

## 
Read 0.0% of 585348 rows
Read 585348 rows and 3 (of 3) columns from 0.012 GB file in 00:00:07

2.4 Output files

Two main output methods are provided. Here, we discuss the first: exporting to disk. However, it is also possible to export to a GenomeInteractions object, as described in Further downstream analysis.

You can export the results to disk, using exportResults(). (If you use runChicago.R, the files appear in ./<results-folder>/data). By default, the function outputs three different output file formats:

outputDirectory <- tempdir()
exportResults(cd, file.path(outputDirectory, "vignetteOutput"))

## Reading the restriction map file...

## Reading the bait map file...

## Preparing the output table...

## Writing out for seqMonk...

## Writing out interBed...

## Preprocessing for WashU outputs...

## Writing out text file for WashU browser upload...

Each called interaction is assigned a score that represents how strong CHiCAGO believes the interaction is (formally, it is based on -log(adjusted P-value)). Thus, a larger score represents a stronger interaction. In each case, the score threshold of 5 is applied.

Summary of output files:

##   bait_chr bait_start bait_end                  bait_name otherEnd_chr
## 1       20     119103   138049                    DEFB126           20
## 2       20     119103   138049                    DEFB126           20
## 3       20     161620   170741                    DEFB128           20
## 4       20     233983   239479                    DEFB132           20
## 5       20     268639   284501 AL034548.1;C20orf96;ZCCHC3           20
## 6       20     268639   284501 AL034548.1;C20orf96;ZCCHC3           20
##   otherEnd_start otherEnd_end otherEnd_name N_reads score
## 1         161620       170741       DEFB128      11  5.07
## 2         523682       536237       CSNK2A1       6  6.78
## 3          73978        76092             .      16  5.11
## 4         206075       209203       DEFB129      33  5.96
## 5         293143       304037             .      34  7.36
## 6         304038       305698             .      34  8.95

##   V1     V2     V3      V4 V5   V6
## 1 20 119103 138049 DEFB126 11 5.07
## 2 20 161620 170741 DEFB128 11 5.07
## 3 20 119103 138049 DEFB126  6 6.78
## 4 20 523682 536237 CSNK2A1  6 6.78
## 5 20 161620 170741 DEFB128 16 5.11
## 6 20  73978  76092       . 16 5.11

##                    V1                  V2   V3
## 1 chr20,119103,138049 chr20,161620,170741 5.07
## 2 chr20,119103,138049 chr20,523682,536237 6.78
## 3 chr20,161620,170741   chr20,73978,76092 5.11
## 4 chr20,233983,239479 chr20,206075,209203 5.96
## 5 chr20,268639,284501 chr20,293143,304037 7.36
## 6 chr20,268639,284501 chr20,304038,305698 8.95

2.5 Visualising interactions

The plotBaits() function can be used to plot the raw read counts versus linear distance from bait for either specific or random baits, labelling significant interactions in a different colour. By default, 16 random baits are plotted, with interactions within 1 Mb from bait passing the threshold of 5 shown in red and those passing the more lenient threshold of 3 shown in blue.

plottedBaitIDs <- plotBaits(cd, n=6)

2.6 Peak enrichment for features

peakEnrichment4Features() tests the hypothesis that other ends in the CHiCAGO output are enriched for genomic features of interest - for example, histone marks associated with enhancers. We find out how many overlaps are expected under the null hypothesis (i.e. that there is no enrichment) by shuffling the other ends around in the genome, while preserving the overall distribution of distances over which interactions span.

You will need additional files to perform this analysis - namely, a .bed file for each feature. We include ChIP-seq data from the ENCODE consortium (The ENCODE Project Consortium 2012), also restricted to chr20 and chr21. (Data accession numbers: Bernstein lab GSM733752, GSM733772, GSM733708, GSM733664, GSM733771, GSM733758)

First, we find the folder that contains the features, and construct a list of the features to use:

featuresFolder <- file.path(dataPath, "GMfeatures")
dir(featuresFolder)

## [1] "featuresGM.txt"                                                 
## [2] "spp.wgEncodeBroadHistoneGm12878CtcfStdAln_chr20and21.narrowPeak"
## [3] "wgEncodeBroadHistoneGm12878H3k27acStdAln_chr20and21.narrowPeak" 
## [4] "wgEncodeBroadHistoneGm12878H3k27me3StdAln_chr20and21.narrowPeak"
## [5] "wgEncodeBroadHistoneGm12878H3k4me1StdAln_chr20and21.narrowPeak" 
## [6] "wgEncodeBroadHistoneGm12878H3k4me3StdAln_chr20and21.narrowPeak" 
## [7] "wgEncodeBroadHistoneGm12878H3k9me3StdAln_chr20and21.narrowPeak"

featuresFile <- file.path(featuresFolder, "featuresGM.txt")
featuresTable <- read.delim(featuresFile, header=FALSE, as.is=TRUE)
featuresList <- as.list(featuresTable$V2)
names(featuresList) <- featuresTable$V1
featuresList

## $CTCF
## [1] "spp.wgEncodeBroadHistoneGm12878CtcfStdAln_chr20and21.narrowPeak"
## 
## $H3K4me1
## [1] "wgEncodeBroadHistoneGm12878H3k4me1StdAln_chr20and21.narrowPeak"
## 
## $H3K4me3
## [1] "wgEncodeBroadHistoneGm12878H3k4me3StdAln_chr20and21.narrowPeak"
## 
## $H3k27ac
## [1] "wgEncodeBroadHistoneGm12878H3k27acStdAln_chr20and21.narrowPeak"
## 
## $H3K27me3
## [1] "wgEncodeBroadHistoneGm12878H3k27me3StdAln_chr20and21.narrowPeak"
## 
## $H3K9me3
## [1] "wgEncodeBroadHistoneGm12878H3k9me3StdAln_chr20and21.narrowPeak"

Next, we feed this information into the peakEnrichment4Features() function.

As part of the analysis, peakEnrichment4Features() takes a distance range (by default, the full distance range over which interactions are observed), and divides it into some number of bins. We must select the number of bins; here, we choose that number to ensure that the bin size is approximately 10kb. If the defaults are changed, a different number of bins is more appropriate. See ?peakEnrichment4Features for more information.

no_bins <- ceiling(max(abs(intData(cd)$distSign), na.rm = TRUE)/1e4)

enrichmentResults <- peakEnrichment4Features(cd, folder=featuresFolder,
              list_frag=featuresList, no_bins=no_bins, sample_number=100)

Note the plot produced by this function. For each feature type, the yellow bar represents the number of features that overlap with interaction other ends. The blue bar represents what would be expected by chance, with a 95% confidence interval for the mean number of overlaps plotted. If the yellow bar lies outside of this interval, we reject the null hypothesis, thus concluding that there is enrichment/depletion of that feature.

The information displayed in the plot is also returned in tabular form (OL = Overlap, SI = Significant Interactions, SD = Standard Deviation, CI = Confidence Interval):

enrichmentResults

##          OLwithSI MeanOLwithSamples SDOLwithSample   LowerCI  HigherCI
## CTCF          349            123.17      10.234676 103.11004 143.22996
## H3K4me1       680            273.56      15.709689 242.76901 304.35099
## H3K4me3       369            136.05       9.541377 117.34890 154.75110
## H3k27ac       464            164.68      11.712792 141.72293 187.63707
## H3K27me3       70             72.62       7.397488  58.12092  87.11908
## H3K9me3       229            120.57      11.384248  98.25687 142.88313

2.7 Further downstream analysis

We can perform further downstream analysis in R or Bioconductor, using functionality from the GenomicInteractions package. First, we export the significant interactions into a GenomicInteractions object:

library(GenomicInteractions)
library(GenomicRanges)
gi <- exportToGI(cd)

From here, we can pass the CHiCAGO results through to other Bioconductor functionality. In the following example, we find out which other ends overlap with the H3K4me1 enhancer mark, using ENCODE data. We use AnnotationHub to fetch a relevant enhancer mark track from the ENCODE project:

library(AnnotationHub)
ah <- AnnotationHub()
hs <- query(ah, c("GRanges", "EncodeDCC", "Homo sapiens", "H3k4me1"))
enhancerTrack <- hs[["AH23254"]]

Next, we use the anchorTwo() function to extract the other end locations from the GenomicInteractions object (anchorOne() would give us the bait locations instead). Note that in this particular instance, the seqlevels() also need to be changed before performing the comparison, adding “chr” to make them match those of the annotation.

otherEnds <- anchorTwo(gi)
otherEnds <- renameSeqlevels(otherEnds, c("chr20","chr21"))

Finally, we look at which other ends overlap the enhancer marks:

findOverlaps(otherEnds, enhancerTrack)

## Hits object with 4534 hits and 0 metadata columns:
##          queryHits subjectHits
##          <integer>   <integer>
##      [1]         2       40732
##      [2]         2       40735
##      [3]         8       40717
##      [4]         8       40718
##      [5]         9       40721
##      ...       ...         ...
##   [4530]      4162       17634
##   [4531]      4162       17635
##   [4532]      4165       17660
##   [4533]      4165       17661
##   [4534]      4165       17662
##   -------
##   queryLength: 4166 / subjectLength: 109612

Further note that the annotation’s genome version should match that of the promoter capture data, namely hg19:

hs["AH23254"]$genome

## [1] "hg19"

2.8 The chicagoData object

In the above workflows, cd is a chicagoData object. It contains three elements:

• intData(cd) is a data.table (note: not a data.frame) that contains information about fragment pairs.
• settings(cd) is a list of settings, usually set with the setExperiment() function.
• params(cd) is a list of parameters. This list is populated as the pipeline runs, and CHiCAGO estimates them in turn.

A closer look at intData(cd):

head(intData(cd), 2)

##    baitID otherEndID distbin       s_j otherEndLen distSign isBait2bait N.1
## 1: 403463     403833      NA 0.2368791        2579  1652804       FALSE   0
## 2: 403463     403843      NA 0.2368791        6302  1690808       FALSE   0
##    N.2 N.3 N refBinMean       s_i NNb NNboe   tlb      tblb        Tmean
## 1:   1   0 1         NA 0.9014785   4     4 [0,1] [  2, 46) 0.0004199299
## 2:   1   0 1         NA 0.9987130   4     4 (8,9] [  2, 46) 0.0038708177
##         Bmean     log.p    log.w     log.q score
## 1: 0.08516784 -2.503705 1.406451 -3.910156     0
## 2: 0.09212130 -2.394217 1.350609 -3.744826     0

Columns:

• baitID: ID of baited fragment
• otherEndID: ID of other end fragment
• s_j: bait-specific scaling factor (Brownian component)
• otherEndLen: The length of the other end fragment
• distSign: The distance from the baited fragment to the other end fragment. Positive and negative values indicate that the other end is, respectively, \(5'\) and \(3'\) of the baited fragment. NA indicates a trans interactions.
• isBait2Bait: TRUE if the other end fragment is also a baited fragment
• N.1, N.2, …: Raw read counts per replicate (see ?mergeSamples).
• N: Merged count (see ?mergeSamples) or raw count in the case of single-replicate interaction calling.
• refBinMean: Can be ignored. (see ?normaliseBaits)
• s_i: other end-specific scaling factor (Brownian component)
• NNb: “N normalised for baits”, a count scaled up by accounting for s_j. May be useful for visualization.
• NNboe: “N normalised for baits and other ends”; may be useful for visualization.
• tlb: Class of other end, based on the number of fragments on other chromosomes that have read pairs.
• tblb: As tlb, for the bait fragment.
• Tmean: Expected count from technical noise.
• Bmean: Expected count from Brownian component. (Thus, the expected count under the null hypothesis is Tmean + Bmean.)
• log.p: p-value associated with fragment pair, on log-scale.
• log.w: p-value weight, on log-scale.
• log.q: weighted p-value, on log-scale.
• score: Final CHiCAGO score.

WARNING: Many functions in CHiCAGO update intData(cd) by reference, which means that intData(cd) can change even when you do not explicitly assign to it. To avoid this behaviour, copy the chicagoData object first:

newCd = copyCD(cd)

##Using different weights

CHiCAGO uses a p-value weighting procedure to upweight proximal interactions and downweight distal interactions. This procedure has four tuning parameters.

The default values of these tuning parameters were calibrated on calls from seven human Macrophage data sets. Provided that your cell type is not too dissimilar to these calibration data, it should be fine to leave the parameters at their default settings. However, if your data set is from a different species or an unusual cell type, you may wish to recalibrate these parameters using data from cell types similar to yours. You can do this with the fitDistCurve.R script in chicagoTools, which we demonstrate in this section.

First, run all of the samples through chicagoPipeline(), saving each chicagoData object in individual .rds files (see saveRDS()). Alternatively, if you are using the runChicago.R wrapper, .rds files should be generated automatically.

Second, run the fitDistCurve.R script. As an example, if we had three biological replicates of mESC cells, we might run the following script at the Unix command prompt:

Rscript fitDistCurve.R mESC --inputs mESC1.rds,mESC2.rds,mESC3.rds 

This script produces the file mESC.settingsFile, which you can read in to modifySettings() as usual - see the Input files required section.

Additionally, the script produces a plot (in this case, called mESC_mediancurveFit.pdf) that can be used to diagnose unreliable estimates. By default, five coloured lines are shown, each representing a parameter estimate from a subset of the data. An unreliable fit is typically diagnosed when the coloured lines are highly dissimilar to each other, and thus the black median line is not representative of them. (Some dissimilarity may be OK, since the median confers robustness.)

For the user’s convenience, a set of precomputed weights are also provided in the package:

weightsPath <- file.path(system.file("extdata", package="Chicago"),
                         "weights")
dir(weightsPath)

## [1] "GM12878-2reps.settings"         "humanMacrophage-7reps.settings"
## [3] "mESC-2reps.settings"

For example, to use the GM12878 weights, supply the appropriate settings file to setExperiment() as per the following:

weightSettings <- file.path(weightsPath, "GM12878-2reps.settings")
cd <- setExperiment(designDir = testDesignDir, settingsFile = weightSettings)

##Session info

sessionInfo()

## R version 3.5.0 (2018-04-23)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows Server 2012 R2 x64 (build 9600)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=C                          
## [2] LC_CTYPE=English_United States.1252   
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.1252    
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] BSgenome.Hsapiens.UCSC.hg19_1.4.0 BSgenome_1.48.0                  
##  [3] Biostrings_2.48.0                 XVector_0.20.0                   
##  [5] rtracklayer_1.40.0                AnnotationHub_2.12.0             
##  [7] GenomicInteractions_1.14.0        InteractionSet_1.8.0             
##  [9] SummarizedExperiment_1.10.0       DelayedArray_0.6.0               
## [11] BiocParallel_1.14.0               matrixStats_0.53.1               
## [13] Biobase_2.40.0                    GenomicRanges_1.32.0             
## [15] GenomeInfoDb_1.16.0               IRanges_2.14.0                   
## [17] S4Vectors_0.18.0                  BiocGenerics_0.26.0              
## [19] PCHiCdata_1.7.0                   Chicago_1.8.0                    
## [21] data.table_1.10.4-3               BiocStyle_2.8.0                  
## 
## loaded via a namespace (and not attached):
##  [1] colorspace_1.3-2              rprojroot_1.3-2              
##  [3] biovizBase_1.28.0             htmlTable_1.11.2             
##  [5] base64enc_0.1-3               dichromat_2.0-0              
##  [7] rstudioapi_0.7                bit64_0.9-7                  
##  [9] interactiveDisplayBase_1.18.0 AnnotationDbi_1.42.0         
## [11] splines_3.5.0                 knitr_1.20                   
## [13] Formula_1.2-2                 Rsamtools_1.32.0             
## [15] cluster_2.0.7-1               shiny_1.0.5                  
## [17] compiler_3.5.0                httr_1.3.1                   
## [19] backports_1.1.2               assertthat_0.2.0             
## [21] Matrix_1.2-14                 lazyeval_0.2.1               
## [23] later_0.7.1                   acepack_1.4.1                
## [25] htmltools_0.3.6               prettyunits_1.0.2            
## [27] tools_3.5.0                   bindrcpp_0.2.2               
## [29] igraph_1.2.1                  gtable_0.2.0                 
## [31] glue_1.2.0                    GenomeInfoDbData_1.1.0       
## [33] dplyr_0.7.4                   Rcpp_0.12.16                 
## [35] xfun_0.1                      stringr_1.3.0                
## [37] mime_0.5                      ensembldb_2.4.0              
## [39] XML_3.98-1.11                 zlibbioc_1.26.0              
## [41] MASS_7.3-50                   scales_0.5.0                 
## [43] VariantAnnotation_1.26.0      BiocInstaller_1.30.0         
## [45] promises_1.0.1                ProtGenerics_1.12.0          
## [47] AnnotationFilter_1.4.0        RColorBrewer_1.1-2           
## [49] yaml_2.1.18                   curl_3.2                     
## [51] memoise_1.1.0                 gridExtra_2.3                
## [53] ggplot2_2.2.1                 biomaRt_2.36.0               
## [55] rpart_4.1-13                  latticeExtra_0.6-28          
## [57] stringi_1.1.7                 RSQLite_2.1.0                
## [59] checkmate_1.8.5               GenomicFeatures_1.32.0       
## [61] rlang_0.2.0                   pkgconfig_2.0.1              
## [63] bitops_1.0-6                  evaluate_0.10.1              
## [65] lattice_0.20-35               bindr_0.1.1                  
## [67] GenomicAlignments_1.16.0      htmlwidgets_1.2              
## [69] bit_1.1-12                    plyr_1.8.4                   
## [71] magrittr_1.5                  bookdown_0.7                 
## [73] R6_2.2.2                      Hmisc_4.1-1                  
## [75] Delaporte_6.1.0               DBI_0.8                      
## [77] pillar_1.2.2                  foreign_0.8-70               
## [79] survival_2.42-3               RCurl_1.95-4.10              
## [81] nnet_7.3-12                   tibble_1.4.2                 
## [83] rmarkdown_1.9                 progress_1.1.2               
## [85] grid_3.5.0                    blob_1.1.1                   
## [87] digest_0.6.15                 xtable_1.8-2                 
## [89] httpuv_1.4.1                  munsell_0.4.3                
## [91] Gviz_1.24.0

##References