# NADfinder Guide

Jianhong Ou, Jun Yu, Paul Kaufman and Lihua Julie Zhu

#### 30 October 2017

#### Abstract

Nucleolus is an important structure inside the nucleus in eukaryotic cells. It is the site for transcribing rDNA into rRNA and for assembling ribosomes, aka ribosome biogenesis. In addition, nucleoli are dynamic hubs through which numerous proteins shuttle and contact specific non-rDNA genomic loci. Deep sequencing analyses of DNA associated with isolated nucleoli (NAD- seq) have shown that specific loci, termed nucleolar- associated domains (NADs) form frequent three- dimensional associations with nucleoli. NAD-seq has been used to study the biological functions of NAD and the dynamics of NAD distribution during embryonic stem cell (ESC) differentiation. Here, we developed a Bioconductor package NADfinder for bioinformatic analysis of the NAD-seq data, including normalization, smoothing, peak calling, peak trimming and annotation.

#### Package

NADfinder 1.2.0

# Contents

# 1 Introduction

Nucleoli serve as major organizing hubs for the three-dimensional structure of mammalian heterochromatin1. Specific loci, termed nucleolar-associated domains (NADs), form frequent three-dimensional associations with nucleoli2. Early mammalian development is a critical period to study NAD biological function, because interactions between pericentromeric chromatin and perinucleolar regions are particularly dynamic during embryonic development4. We therefore propose for the first time to map the NADs in the mouse genome, determine how these associations are altered during embryonic stem cell (ESC) differentiation, and develop tools for study of these higher-order chromosome interactions in fixed and live single cells (Fig. 1).

Fig. 1 Workflow for NAD peak identification. Read counts are summarized for each 50kb moving window with step size of 1kb for nucleolar and genomic samples. Log2 ratio between nucleolar and genomic samples was computed for each window, followed by local background correction, smoothing, peak calling, filtering and annotation.

# 2 Single sample

Here is an example to use NADfinder for peak calling.

## 2.1 Coverage calculation

Here is the code snippet for calculating coverage with a sliding window in a given step along the genome using a pair of bam files from genomic sample as background and purified nucleoli as target signal.

    ## load the library
    library(NADfinder)
    library(SummarizedExperiment)
    ## bam file path
    path <- "path/to/your/bam/files"
    f <- c(genome="genome.bam", nucleosome="nucleosome.bam")
    ## window size for tile of genome. Here we set it to a big window size,
    ## ie., 50k because of the following two reasons:
    ## 1. peaks of NADs are wide;
    ## 2. data will be smoothed better with bigger window size.
    ws <- 50000
    ## step for tile of genome
    step <- 1000
    ## Set the background.
    ## 0.25 means 25% of the lower ratios will be used for background training.
    backgroundPercentage <- 0.25
    ## Count the reads for each window with a given step in the genome.
    ## The output will be a GRanges object.
    library(BSgenome.Mmusculus.UCSC.mm10)

    se <- tileCount(reads=file.path(path, f),
                    genome=Mmusculus,
                    windowSize=ws,
                    step=step)

Here we load the pre-computed coverage data single.count to save build time.

    data(single.count)
    se <- single.count

For quality asessment, `cumulativePercentage` can be used to plot the cumulative sums of sorted read counts for nucleosome and genomic DNA. We expect the cumulative sum in the genomic DNA to be close to a straight line because the coverage for the genomic DNA sample should be uniformly distributed. Unlike ChIP-seq data, the cumulative sum in the nucleosome sample will not exhibit sharp changes because most of NADs are broad regions as wide as 100 kb. However, we should observe clear differences between the two curves.

    ## Calculate ratios for peak calling. We use nucleosome vs genomic DNA.
    dat <- log2se(se, nucleosomeCols = "nucleosome.bam", genomeCols="genome.bam")
    ## Smooth the ratios for each chromosome.
    ## We found that for each chromosome, the ratios are higher in 5'end than 3'end.
    datList <- smoothRatiosByChromosome(dat, N=100)
    ## check the difference between the cumulative percentage tag allocation
    ## in genome and nucleosome samples.
    cumulativePercentage(datList[["chr18"]])

## 2.2 Call peaks

Before peak calling, the function `smoothRatiosByChromosome` is used for log ratios calculation, normalization and smoothing.

The peaks will be called if the ratios are significantly higher than chromosome-specific background determined by `trimPeaks`. The following figure shows the peaks (black rectangles) called using normalized (green curve) and smoothed (red curve) log2 ratios.

    ## check the results of smooth function
    chr18 <- datList[["chr18"]] ## we only have reads in chr18 in test data.
    chr18subset <- subset(chr18, seq.int(floor(length(chr18)/10))*10)
    chr18 <- assays(chr18subset)
    ylim <- range(c(chr18$ratio[, 1],
                    chr18$bcRatio[, 1],
                    chr18$smoothedRatio[, 1]))
    plot(chr18$ratio[, 1],
         ylim=ylim*c(.9, 1.1),
         type="l", main="chr18")
    abline(h=0, col="cyan", lty=2)
    points(chr18$bcRatio[, 1], type="l", col="green")
    points(chr18$smoothedRatio[, 1], type="l", col="red")
    legend("topright",
           legend = c("raw_ratio", "background_corrected", "smoothed"),
           fill = c("black", "green", "red"))
    ## call peaks for each chromosome
    peaks <- lapply(datList, trimPeaks,
                    backgroundPercentage=backgroundPercentage,
                    cutoffPvalue=0.05, countFilter=1000)
    ## plot the peaks in "chr18"
    peaks.subset <- countOverlaps(rowRanges(chr18subset), peaks$chr18)>=1
    peaks.run <- rle(peaks.subset)
    run <- cumsum(peaks.run$lengths)
    run <- data.frame(x0=c(1, run[-length(run)]), x1=run)
    run <- run[peaks.run$values, , drop=FALSE]
    rect(xleft = run$x0, ybottom = ylim[2]*.75,
         xright = run$x1, ytop = ylim[2]*.8,
         col = "black")

    ## convert list to a GRanges object
    peaks.gr <- unlist(GRangesList(peaks))

The following shows how to save or export the called peaks for downstream analysis.

    ## output the peaks
    write.csv(as.data.frame(unname(peaks.gr)), "peaklist.csv", row.names=FALSE)
    ## export peaks to a bed file.
    library(rtracklayer)
    export(peaks.gr, "peaklist.bed", "BED")

# 3 Samples with duplicates

Data analysis with multiple biological replicates follows the same steps as that of a single paired samples, i.e., coverage calculation, normalization and smoothing, and peak calling. The only difference is that limma is used to determine the statistical significance in peak calling.

    library(NADfinder)
    ## bam file path
    path <- "path/to/your/bam/files"
    f <- c("G26.bam", "G28.bam", "G29.bam", "N26.bam", "N28.bam", "N29.bam")
    ws <- 50000
    step <- 1000
    library(BSgenome.Mmusculus.UCSC.mm10)

    se <- tileCount(reads=file.path(path, f),
                    genome=Mmusculus,
                    windowSize=ws,
                    step=step)

Load saved coverage.

    data(triplicates.counts)
    se <- triplicates.counts

    ## Calculate ratios for nucleoli vs genomic sample.
    gps <- c("26", "28", "29")
    se <- log2se(se,
                 nucleosomeCols = paste0("N", gps, ".bam"),
                 genomeCols = paste0("G", gps, ".bam"))
    getCorrelations(se, "chr18")

    ## $cor
    ##           N26.bam   N28.bam   N29.bam
    ## N26.bam 1.0000000 0.8297445 0.7848867
    ## N28.bam 0.8297445 1.0000000 0.8001389
    ## N29.bam 0.7848867 0.8001389 1.0000000
    ##
    ## $coe
    ##           N26.bam   N28.bam   N29.bam
    ## N26.bam 1.0000000 0.7584727 0.6677331
    ## N28.bam 0.7584727 1.0000000 0.7111284
    ## N29.bam 0.6677331 0.7111284 1.0000000

    seList<- smoothRatiosByChromosome(se, chr="chr18")
    cumulativePercentage(seList[["chr18"]])

    peaks <- lapply(seList, callPeaks,
                    cutoffAdjPvalue=0.05, countFilter=1000)
    peaks <- unlist(GRangesList(peaks))
    peaks

    ## GRanges object with 115 ranges and 3 metadata columns:
    ##             seqnames               ranges strand |            AveSig
    ##                <Rle>            <IRanges>  <Rle> |         <numeric>
    ##     chr18.1    chr18   [2980989, 3235987]      * |  1.79742085885062
    ##     chr18.2    chr18   [3408987, 4205983]      * |  1.61510181933636
    ##     chr18.3    chr18   [4390983, 4644981]      * |  1.08469616537963
    ##     chr18.4    chr18   [4709982, 4926980]      * | 0.802259931424113
    ##     chr18.5    chr18   [5187980, 5340978]      * | 0.404022193835048
    ##         ...      ...                  ...    ... .               ...
    ##   chr18.111    chr18 [87132654, 87642651]      * | 0.409729968809612
    ##   chr18.112    chr18 [87735651, 88736646]      * | 0.518711810111407
    ##   chr18.113    chr18 [89213646, 89776642]      * | 0.527758540125288
    ##   chr18.114    chr18 [89833643, 90447640]      * | 0.565759644363372
    ##   chr18.115    chr18 [90556640, 90605639]      * | 0.503765262908673
    ##                          P.value            adj.P.Val
    ##                        <numeric>            <numeric>
    ##     chr18.1 2.27231525229115e-09  1.0368184523373e-08
    ##     chr18.2  5.5110358933687e-09 2.03394311719841e-08
    ##     chr18.3 5.48463420484863e-07 9.01006631303913e-07
    ##     chr18.4  4.2553510530638e-06 6.08701529569145e-06
    ##     chr18.5 0.000411348600562326 0.000484469838566605
    ##         ...                  ...                  ...
    ##   chr18.111 0.000195241367274636 0.000234565984267849
    ##   chr18.112 2.61153158181348e-05 3.41977340447212e-05
    ##   chr18.113 2.23409068310891e-05 2.94844781505791e-05
    ##   chr18.114 1.17783744408087e-05 1.60430349147631e-05
    ##   chr18.115 3.36282562512596e-05 4.34630589026378e-05
    ##   -------
    ##   seqinfo: 66 sequences from an unspecified genome; no seqlengths

The peaks can be visualized along the ideogram using `plotSig`.

    ideo <- readRDS(system.file("extdata", "ideo.mm10.rds",
                                package = "NADfinder"))
    plotSig(ideo=ideo, grList=GRangesList(peaks), mcolName="AveSig",
            layout=list("chr18"),
            parameterList=list(types="heatmap"))

# 4 Session Info

    sessionInfo()

    ## R version 3.4.2 Patched (2017-10-07 r73498)
    ## 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] stats4    parallel  stats     graphics  grDevices utils     datasets
    ## [8] methods   base
    ##
    ## other attached packages:
    ##  [1] SummarizedExperiment_1.8.0         DelayedArray_0.4.0
    ##  [3] matrixStats_0.52.2                 Biobase_2.38.0
    ##  [5] BSgenome.Mmusculus.UCSC.mm10_1.4.0 BSgenome_1.46.0
    ##  [7] rtracklayer_1.38.0                 Biostrings_2.46.0
    ##  [9] XVector_0.18.0                     NADfinder_1.2.0
    ## [11] GenomicRanges_1.30.0               GenomeInfoDb_1.14.0
    ## [13] IRanges_2.12.0                     S4Vectors_0.16.0
    ## [15] BiocGenerics_0.24.0                BiocStyle_2.6.0
    ##
    ## loaded via a namespace (and not attached):
    ##  [1] ProtGenerics_1.10.0           bitops_1.0-6
    ##  [3] bit64_0.9-7                   RColorBrewer_1.1-2
    ##  [5] progress_1.1.2                httr_1.3.1
    ##  [7] rprojroot_1.2                 tools_3.4.2
    ##  [9] backports_1.1.1               R6_2.2.2
    ## [11] rpart_4.1-11                  Hmisc_4.0-3
    ## [13] DBI_0.7                       lazyeval_0.2.1
    ## [15] Gviz_1.22.0                   colorspace_1.3-2
    ## [17] nnet_7.3-12                   gridExtra_2.3
    ## [19] prettyunits_1.0.2             RMySQL_0.10.13
    ## [21] curl_3.0                      bit_1.1-12
    ## [23] compiler_3.4.2                htmlTable_1.9
    ## [25] SparseM_1.77                  grImport_0.9-0
    ## [27] bookdown_0.5                  scales_0.5.0
    ## [29] checkmate_1.8.5               pbapply_1.3-3
    ## [31] stringr_1.2.0                 digest_0.6.12
    ## [33] Rsamtools_1.30.0              foreign_0.8-69
    ## [35] rmarkdown_1.6                 base64enc_0.1-3
    ## [37] dichromat_2.0-0               htmltools_0.3.6
    ## [39] ensembldb_2.2.0               limma_3.34.0
    ## [41] htmlwidgets_0.9               rlang_0.1.2
    ## [43] RSQLite_2.0                   BiocInstaller_1.28.0
    ## [45] shiny_1.0.5                   BiocParallel_1.12.0
    ## [47] acepack_1.4.1                 VariantAnnotation_1.24.0
    ## [49] RCurl_1.95-4.8                magrittr_1.5
    ## [51] GenomeInfoDbData_0.99.1       Formula_1.2-2
    ## [53] Matrix_1.2-11                 Rcpp_0.12.13
    ## [55] munsell_0.4.3                 trackViewer_1.14.0
    ## [57] stringi_1.1.5                 yaml_2.1.14
    ## [59] MASS_7.3-47                   zlibbioc_1.24.0
    ## [61] plyr_1.8.4                    AnnotationHub_2.10.0
    ## [63] grid_3.4.2                    blob_1.1.0
    ## [65] lattice_0.20-35               splines_3.4.2
    ## [67] GenomicFeatures_1.30.0        knitr_1.17
    ## [69] biomaRt_2.34.0                XML_3.98-1.9
    ## [71] evaluate_0.10.1               biovizBase_1.26.0
    ## [73] latticeExtra_0.6-28           data.table_1.10.4-3
    ## [75] httpuv_1.3.5                  gtable_0.2.0
    ## [77] assertthat_0.2.0              ggplot2_2.2.1
    ## [79] mime_0.5                      xtable_1.8-2
    ## [81] AnnotationFilter_1.2.0        baseline_1.2-1
    ## [83] survival_2.41-3               signal_0.7-6
    ## [85] tibble_1.3.4                  GenomicAlignments_1.14.0
    ## [87] AnnotationDbi_1.40.0          memoise_1.1.0
    ## [89] cluster_2.0.6                 interactiveDisplayBase_1.16.0

# References

1. Politz, J. C. R., Scalzo, D. & Groudine, M. The redundancy of the mammalian heterochromatic compartment. Current opinion in genetics & development 37, 1–8 (2016).

2. Németh, A. et al. Initial genomics of the human nucleolus. PLoS Genet 6, e1000889 (2010).

3. Koningsbruggen, S. van et al. High-resolution whole-genome sequencing reveals that specific chromatin domains from most human chromosomes associate with nucleoli. Molecular biology of the cell 21, 3735–3748 (2010).

4. Aguirre-Lavin, T. et al. 3D-fish analysis of embryonic nuclei in mouse highlights several abrupt changes of nuclear organization during preimplantation development. BMC developmental biology 12, 1 (2012).

5. Popken, J. et al. Reprogramming of fibroblast nuclei in cloned bovine embryos involves major structural remodeling with both striking similarities and differences to nuclear phenotypes of in vitro fertilized embryos. Nucleus 5, 555–589 (2014).