diff --git a/chp08.md b/chp08.md
index 1459cd1..b8a03ca 100644
--- a/chp08.md
+++ b/chp08.md
@@ -15,31 +15,79 @@ coldata_file <- system.file("extdata/rna-seq/SRP029880.colData.tsv",
 
 **solution:**
 ```{r,echo=FALSE,eval=FALSE}
-#coming soon
- 
+# first compare total counts for each sample
+counts <- as.matrix(read.table(counts_file, header = T, sep = '\t'))
+colSums(counts) / 1E9 # display billions of reads for each sample 
+# calculate TPM (transcripts per million)
+#find gene length normalized values 
+geneLengths <- as.vector(subset(counts, select = c(width)))
+rpk <- apply(subset(counts, select = c(-width)), 
+             2, 
+             function(x) {x/(geneLengths/1000)}
+            )
+#normalize by the sample size using rpk values
+tpm <- apply(rpk, 
+             2, 
+             function(x) {x / sum(as.numeric(x)) * 10^6}
+             )
+colSums(tpm) # confirm that normalized values sum to 10^6
+head(tpm) # glimpse data
+
 ```
 
 2. Plot a heatmap of the top 500 most variable genes. Compare with the heatmap obtained using the 100 most variable genes. [Difficulty: **Beginner**]
 
 **solution:**
 ```{r,echo=FALSE,eval=FALSE}
-#coming soon
+#compute the variance of each gene across samples
+V <- apply(X = tpm, 
+           MARGIN = 1, 
+           FUN = var)
+
+# sort the results by variance in decreasing order 
+selectedGenes <- names(V[order(V, decreasing = T)])
+# make heatmaps clustering genes and samples
+library(pheatmap)
+# add annotations for samples 
+colData <- read.table(coldata_file, header = T, sep = '\t', 
+                      stringsAsFactors = TRUE)
+# heatmap with the top 100 genes 
+pheatmap(tpm[selectedGenes[1:100], ], scale = 'row', show_rownames = FALSE,
+         annotation_col = colData)
+# heatmap with the top 500 genes 
+pheatmap(tpm[selectedGenes[1:500], ], scale = 'row', show_rownames = FALSE)
  
 ```
 
+*Clustering with the TPM values of the top 100 genes or top 500 genes produces a tree that separately clusters the control and case samples. However, there is some difference in which samples cluster together between the two heatmaps. For example, Case 1 is most closely related to case 4 in the 100 gene tree, but Case 1 is most closely related to Case 3 in the 500 gene tree.*
+
+
 3. Re-do the heatmaps setting the `scale` argument to `none`, and `column`. Compare the results with `scale = 'row'`. [Difficulty: **Beginner**]
 
 **solution:**
 ```{r,echo=FALSE,eval=FALSE}
-#coming soon
+pheatmap(tpm[selectedGenes[1:500], ], show_rownames = FALSE)
+pheatmap(tpm[selectedGenes[1:500], ], scale = 'column', show_rownames = FALSE)
  
 ```
 
+*Scaling by rows allows you to compare relative expression levels of each gene across the samples. Without any scaling, only extreme differences in absolute expression levels are noticeable; in particular, the variation among lowly-level genes (TPM < 50000) is not noticeable.  When scaling by columns, the heatmap indicates that nearly all of these genes with similar absolute values show average levels (log = 0) of expression among the nearly 20,000 genes overall.*
+
+
 4. Draw a correlation plot for the samples depicting the sample differences as 'ellipses', drawing only the upper end of the matrix, and order samples by hierarchical clustering results based on `average` linkage clustering method. [Difficulty: **Beginner**]
 
 **solution:**
 ```{r,echo=FALSE,eval=FALSE}
-#coming soon
+correlationMatrix <- cor(tpm)
+# visualize
+library(corrplot)
+corrplot(correlationMatrix, 
+         method = "ellipse", type = "upper",
+         order = 'hclust', 
+         hclust.method = "average",
+         # addrect = 2, 
+         addCoef.col = 'black', 
+         number.cex = 0.7) 
  
 ```
 
@@ -47,7 +95,27 @@ coldata_file <- system.file("extdata/rna-seq/SRP029880.colData.tsv",
 
 **solution:**
 ```{r,echo=FALSE,eval=FALSE}
-#coming soon
+# Heatmaps
+# compute the sum of each gene across samples
+TE <- apply(X = tpm, 
+           MARGIN = 1, 
+           FUN = sum)
+# sort the results by total of expression values in decreasing order 
+selectedGenes_TE <- names(TE[order(TE, decreasing = T)])
+# heatmap with the top 100 genes 
+pheatmap(tpm[selectedGenes_TE[1:100], ], scale = 'row', show_rownames = FALSE)
+
+# PCA
+library(ggfortify)
+# transpose the matrix so genes become columns
+M <- t(tpm[selectedGenes,])
+# transform the counts to log2 scale 
+M <- log2(M + 1) # + 1 to avoid log(0)
+# compute PCA 
+pcaResults <- prcomp(M)
+
+# plot PCA results making use of ggplot2's autoplot function
+autoplot(pcaResults, data = colData, colour = 'group')
  
 ```
 
@@ -55,7 +123,10 @@ coldata_file <- system.file("extdata/rna-seq/SRP029880.colData.tsv",
 
 **solution:**
 ```{r,echo=FALSE,eval=FALSE}
-#coming soon
+colData$batch <- c(rep(letters[1:2],times = 5)) # simulate batch effect
+pheatmap(tpm[selectedGenes_TE[1:100],], scale = 'row', 
+         show_rownames = FALSE, 
+         annotation_col = colData)
  
 ```
 
@@ -99,15 +170,49 @@ Now, you are ready to do the following exercises:
 
 **solution:**
 ```{r,echo=FALSE,eval=FALSE}
-#coming soon
- 
+# First create needed objects (see list above)
+# 1. read count data: remove the 'width' column from previous matrix
+counts <- as.matrix(read.table(counts_file, header = T, sep = '\t'))
+countData <- as.matrix(subset(counts, select = c(-width)))
+# 2. get table with experimental setup 
+colData <- read.table(coldata_file, header = T, sep = '\t', 
+                      stringsAsFactors = TRUE)
+# 3. define the design formula, indicating variable of interest in colData
+designFormula <- "~ group"
+# set up DESeqDataSet object
+library(DESeq2)
+# Now create a DESeq dataset object from the three elements above 
+dds <- DESeqDataSetFromMatrix(countData = countData, 
+                              colData = colData, 
+                              design = as.formula(designFormula))
+# For each gene, we count the total number of reads for that gene in all samples 
+# and only keep those that have at least 2 reads
+dds <- dds[rowSums(DESeq2::counts(dds)) > 1, ]
+# run analysis on this filtered data set
+dds <- DESeq(dds)
+# compute the contrast (difference) for the 'group' variable 
+# where 'CTRL' samples are used as the control group. 
+DEresults = results(dds, contrast = c("group", 'CASE', 'CTRL')) 
+# sort rows by increasing p-value
+DEresults <- DEresults[order(DEresults$pvalue),]
+
+# show the top 10 genes to briefly check results
+head(DEresults, n = 10) 
+
+# make volcano plot
+ggplot(data = as.data.frame(DEresults), 
+       aes(x = log2FoldChange, y = -log10(pvalue))) + 
+  geom_point()
+
 ```
 
 2. Use DESeq2::plotDispEsts to make a dispersion plot and find out the meaning of this plot. (Hint: Type ?DESeq2::plotDispEsts) [Difficulty: **Beginner**]
 
 **solution:**
 ```{r,echo=FALSE,eval=FALSE}
-#coming soon
+plotDispEsts(dds, ymin = 1e-03, # removes empty plot space due to outlier
+             # finalcol = NULL # uncomment to see all original values in black
+             )
  
 ```
 
@@ -115,18 +220,23 @@ Now, you are ready to do the following exercises:
 
 **solution:**
 ```{r,echo=FALSE,eval=FALSE}
-#coming soon
- 
+?DESeq2::results 
+
 ```
 
+*Default value is 0. If the default value is changed to 1, then this performs a hypothesis test that the absolute values of log2 fold changes are less than or equal to 1 (assuming that altHypothesis is at the default setting of greaterAbs).*
+
+
 4. What is independent filtering? What happens if we don't use it? Google `independent filtering statquest` and watch the online video about independent filtering. [Difficulty: **Intermediate**]
 
 **solution:**
 ```{r,echo=FALSE,eval=FALSE}
-#coming soon
  
 ```
 
+*Independent filtering is a method to reduce the number of tested genes to improve the power of statistical testing, the ability to detect true positives. The method removes genes with low expression values, which are more likely to yield false positives due to higher dispersion (see above). By reducing the number of false positives, the proportion of true positives is increased.*
+
+
 5. Re-do the differential expression analysis using the `edgeR` package. Find out how much DESeq2 and edgeR agree on the list of differentially expressed genes. [Difficulty: **Advanced**] 
 
 **solution:**
@@ -150,32 +260,161 @@ Now, you are ready to do the following exercises:
 
 **solution:**
 ```{r,echo=FALSE,eval=FALSE}
-#coming soon
+library(gprofiler2) # updated version of gProfileR 
+library(knitr) # to make nice tables
+
+# filter for genes with significant changes in expression
+#remove genes with NA values 
+DE <- DEresults[!is.na(DEresults$padj),]
+#select genes with adjusted p-values below 0.05
+DE <- DE[DE$padj < 0.05,]
+#select genes with absolute log2 fold change above 1 (two-fold change)
+DE <- DE[abs(DE$log2FoldChange) > 1,]
+
+#get the list of genes of interest
+genesOfInterest <- rownames(DE)
+
+#calculate enriched terms amongst this gene list for multiple data sources
+multiResults <- gost(query = genesOfInterest, # new function with gprofiler2
+                       organism = 'hsapiens', 
+                       sources = c('GO', "KEGG", "REACTOME", "CORUM"))
+kable(multiResults$result) # nicely formats data in table
+
+#sort GO terms for precision, recall or p.value
+goResults <- gost(query = genesOfInterest, # new function with gprofiler2
+                       organism = 'hsapiens', 
+                       sources = 'GO')
+
+# sorted by increasing p_value 
+# based on sizes of query (gene list), term (gene set), and intersection (overlap)
+kable(head(goResults$result[order(goResults$result$p_value),])) # nicely formatted table of top lists
+
+# sorted by decreasing precision, intersection_size/query_size
+kable(head(goResults$result[order(goResults$result$precision),])) # nicely formatted table of top lists
+
+# sorted by decreasing recall, intersection_size/term_size
+kable(head(goResults$result[order(goResults$result$recall, decreasing = T),])) # nicely formatted table of top lists
  
 ```
 
+*With this large gene list (query_size), the results sorted by decreasing precision is similar to the results sorted by increasing p-value. In both cases, enriched terms with very low p-values and gene sets with many entries (large term sizes) are at the top of the list. These gene sets are higher level, such as "protein binding" or "cytoplasm". In contrast, when the results are sorted by decreasing recall, gene sets with fewer entries (small term sizes) are at the top of the list. The latter sort can be useful for focusing on more specific biological features, such as "blood coagulation, fibrin clot formation" and "complement activation, alternative pathway".*
+
+
 2. Repeat the gene set enrichment analysis by trying different options for the `compare` argument of the `GAGE:gage`
 function. How do the results differ? [Difficulty: **Beginner**]
 
 **solution:**
 ```{r,echo=FALSE,eval=FALSE}
-#coming soon
+#use the normalized counts to carry out a GSEA
+normalizedCounts <- DESeq2::counts(dds, normalized = TRUE)
+
+library(gage)
+# use the top term from the GO results to create 1st gene set 
+# but first restrict goResults to terms that have < 100 genes 
+go <- goResults$result[goResults$result$term_size < 100,] 
+geneSet1 <- gconvert(go$term_id[1], organism = "hsapiens") 
+
+# randomly select 25 genes from the counts table for 2nd gene set
+geneSet2 <- sample(rownames(normalizedCounts), 25)
+
+geneSets <- list('top_GO_term' = geneSet1$name,
+                 'random_set' = geneSet2)
+
+# perform GSEA with different options for compare argument
+# 'paired' compares each reference with its matching sample (like paired t-test) 
+# appropriate here since each reference and sample is from same person
+print("compare = 'paired'")
+gseaResults <- gage(exprs = log2(normalizedCounts+1), # avoids log(0)
+                    ref = match(rownames(colData[colData$group == 'CTRL',]), 
+                                colnames(normalizedCounts)), 
+                    samp = match(rownames(colData[colData$group == 'CASE',]), 
+                                 colnames(normalizedCounts)),
+                    gsets = geneSets, 
+                    compare = 'paired')
+gseaResults$greater
+# 'unpaired' compares all possible pairs of ref and samples 
+print("compare = 'unpaired'")
+gseaResults <- gage(exprs = log2(normalizedCounts+1), # avoids log(0)
+                    ref = match(rownames(colData[colData$group == 'CTRL',]), 
+                                colnames(normalizedCounts)), 
+                    samp = match(rownames(colData[colData$group == 'CASE',]), 
+                                 colnames(normalizedCounts)),
+                    gsets = geneSets, 
+                    compare = 'unpaired')
+gseaResults$greater
+# '1ongroup' compares each sample to average of all references
+print("compare = '1ongroup'")
+gseaResults <- gage(exprs = log2(normalizedCounts+1), # avoids log(0)
+                    ref = match(rownames(colData[colData$group == 'CTRL',]), 
+                                colnames(normalizedCounts)), 
+                    samp = match(rownames(colData[colData$group == 'CASE',]), 
+                                 colnames(normalizedCounts)),
+                    gsets = geneSets, 
+                    compare = '1ongroup')
+gseaResults$greater
+# 'as.group' compares the two groups (like two-sample t-test) 
+print("compare = 'as.group'")
+gseaResults <- gage(exprs = log2(normalizedCounts+1), # avoids log(0)
+                    ref = match(rownames(colData[colData$group == 'CTRL',]), 
+                                colnames(normalizedCounts)), 
+                    samp = match(rownames(colData[colData$group == 'CASE',]), 
+                                 colnames(normalizedCounts)),
+                    gsets = geneSets, 
+                    compare = 'as.group')
+gseaResults$greater
  
 ```
 
+*The paired comparison matches each case with its control sample, which accounts for the variation between individuals. The unpaired comparison matches each case with all control samples and produces a similar p-value. The 1ongroup comparison matches each case to the average of all controls.  Finally, the as.group comparison doesn't account for individual variation, and this is reflected in the higher p-value compared to the other comparisons.*
+
+
 3. Make a scatter plot of GO term sizes and obtained p-values by setting the `gProfiler::gprofiler` argument `significant = FALSE`. Is there a correlation of term sizes and p-values? (Hint: Take -log10 of p-values). If so, how can this bias be mitigated? [Difficulty: **Intermediate**]
 
 **solution:**
 ```{r,echo=FALSE,eval=FALSE}
-#coming soon
+goResults_all <- gost(query = genesOfInterest, # new function with gprofiler2
+                  organism = 'hsapiens', 
+                  significant = FALSE,
+                  sources = 'GO')
+ggplot(goResults_all$result, aes(x= sqrt(term_size), y= -log10(p_value))) + 
+  geom_point()
+cor(x = sqrt(goResults_all$result$term_size), # square-root transformation explained below
+    y = -log10(goResults_all$result$p_value))
  
 ```
 
+*The scatter plot (and Pearson correlation coefficient) indicates that there is a correlation between term sizes and p-values. After trying several mathematical transformations of the explanatory variable, this was most evident with a square root transformation of term_size (and a negative log transformation of p_value). To mitigate for this bias, -log10(p-value) could be divided by the square root of the term size.*
+
+
 4. Do a gene-set enrichment analysis using gene sets from top 10 GO terms. [Difficulty: **Intermediate**]
 
 **solution:**
 ```{r,echo=FALSE,eval=FALSE}
-#coming soon
+# retrieve the top 10 GO term identifiers
+go.term_ids <- go$term_id[1:10] 
+# construct gene sets by multiply applying gconvert() for each GO term_id
+geneSets_10 <- mapply(function(x) {
+    gconvert(x, organism = "hsapiens")
+  }, go.term_ids
+)
+geneSets_10
+# geneSets_10 <- go.term_ids %>% 
+library(dplyr) # used for pull() and %>% pipe
+library(purrr) # used for map, an alternative to mapply
+geneSets_10 <- go.term_ids %>% map(function(x) {
+    gconvert(x, organism = "hsapiens") %>%
+        pull("name")
+    }
+)
+names(geneSets_10) <- go.term_ids # name each list element with its GO term_id
+gseaResults_10 <- gage(exprs = log2(normalizedCounts+1), # avoids log(0)
+                       ref = match(rownames(colData[colData$group == 'CTRL',]), 
+                                   colnames(normalizedCounts)), 
+                       samp = match(rownames(colData[colData$group == 'CASE',]), 
+                                    colnames(normalizedCounts)),
+                       gsets = geneSets_10, 
+                       compare = 'paired')
+gseaResults_10$greater
  
 ```
 
@@ -183,10 +422,12 @@ function. How do the results differ? [Difficulty: **Beginner**]
 
 **solution:**
 ```{r,echo=FALSE,eval=FALSE}
-#coming soon
  
 ```
 
+*A general search for "gene set enrichment analysis" at bioconductor.org returned an overwhelming number of entries.  At the bottom of the bioconductor.org page on [gage](http://bioconductor.org/packages/release/bioc/html/gage.html), the biocViews Details includes a link to the (GeneSetEnrichment)(http://bioconductor.org/packages/release/BiocViews.html#___GeneSetEnrichment) category, listing a more manageable number (139) of packages. An internet search for "gene set enrichment analysis with R" identifies additional packages, including gprofiler2.* 
+
+
 6.  Use the topGO package (https://bioconductor.org/packages/release/bioc/html/topGO.html) to re-do the GO term analysis. Compare and contrast the results with what has been obtained using the `gProfileR` package. Which tool is faster, `gProfileR` or topGO? Why? [Difficulty: **Advanced**]
 
 **solution:**
diff --git a/chp09.md b/chp09.md
index 4f93156..4d84856 100644
--- a/chp09.md
+++ b/chp09.md
@@ -6,7 +6,47 @@
 
 **solution:**
 ```{r,echo=FALSE,eval=FALSE}
-#coming soon
+# Some objects must be prepared before starting this exercise. First, get the paths to the data files of interest:
+# get path to all chip-seq datasets and subsets by file format
+data_path <- system.file('extdata/chip-seq',package='compGenomRData') 
+chip_files <- list.files(data_path, full.names=TRUE)
+bam_files <- list.files(data_path, full.names=TRUE, pattern='bam$')
+bw_files <- list.files(data_path, full.names=TRUE, pattern='bw$')
+
+# Next, get human genome data, focusing on chr21, since all files are limited to this chromosome.
+library(GenomeInfoDb) # load the chromosome info package
+hg_chrs <- getChromInfoFromUCSC('hg38') # get human chromosome lengths
+hg_chrs <- subset(hg_chrs, grepl('chr21$',chrom)) # get length of chromosome 21
+seqlengths = with(hg_chrs, setNames(size, chrom)) # convert to named vector
+
+# Now follow the process described in subsection 9.5.4, Plus and minus strand cross-correlation. The maximum cross-correlation value when shifting the strands should correspond to the average DNA fragment length present in the library.* 
+library(GenomicAlignments)
+# use apply function to perform commands on each bam file
+reads_list <- lapply(bam_files, function(x){ # for each bam file, 
+  reads <-  readGAlignments(x)  # load the reads 
+  reads <- granges(reads)  # and convert to a GRanges object
+  reads <- resize(reads, width=1, fix='start') # then set each range to start position
+  reads <- keepSeqlevels(reads, 'chr21', pruning.mode='coarse') # remove extra levels
+})
+# calculate the coverage profile for plus and minus strand for each dataset
+wsize <- 1:400 # define shift range (for later use)
+jaccard = function(x,y)sum((x & y)) / sum((x | y)) # calculates jaccard similarity
+cc_list <- lapply(reads_list, function(x){ # for each GRanges object (dataset), 
+  reads <- split(x, strand(x)) # split each object based on strand
+  cov <- lapply(reads, function(y){
+    coverage(y, width = seqlengths)[[1]] > 0 # convert coverage vector to boolean vector
+    }) 
+  cov <- lapply(cov, as.vector)
+  cc  <- shiftApply(SHIFT = wsize, # shift the + vector by 1 - 400 nucleotides and after each shift
+                    X     = cov[['+']], 
+                    Y     = cov[['-']], 
+                    FUN   = jaccard) # calculate the similarity between strands 
+  })
+
+# convert the results into data frames 
+cc_list <- lapply(cc_list, function(cc){
+  data.frame(fragment_size = wsize, cross_correlation = cc)
+})
  
 ```
 
@@ -14,24 +54,53 @@
 
 **solution:**
 ```{r,echo=FALSE,eval=FALSE}
-#coming soon
- 
+library(ggplot2)
+# make list of experiment names to label plots
+expt_names <- sub('.chr21.bam','', basename(bam_files)) # remove path and common suffix
+expt_names = sub('GM12878_hg38_','', expt_names) # remove common prefix
+# add experiment name to each dataframe in cc_list
+names(cc_list) <- expt_names
+
+lapply(cc_list, function(cc){
+  ggplot(data = cc, aes(fragment_size, cross_correlation)) +
+    geom_point() +
+    geom_vline(xintercept = which.max(cc$cross_correlation), 
+               size=2, color='red', linetype=2) +
+    theme_bw() +
+    theme(
+        axis.text = element_text(size=10, face='bold'),
+        axis.title = element_text(size=14,face="bold"),
+        plot.title = element_text(hjust = 0.5)) +
+    labs(x = 'Shift in base pairs',
+         y = ('Jaccard similarity')
+         )
+}) 
 ```
 
 3. How does the Input sample distribution differ from the ChIP samples? [Difficulty: **Beginner**]
 
 **solution:**
-```{r,echo=FALSE,eval=FALSE}
-#coming soon
- 
-```
+
+*The Jaccard similarity plots for the Input samples (7-11) are more likely to be diffuse with no defined peaks; the major exception is Input_r5 (#11). In contrast, the Jaccard similarity plots for the ChIP samples generally have a discrete peak, indicative of the fragment size; the major exception is H3K4me1, which is quite diffuse. (#4).*
+
 
 4. Write a function which converts the bam files into bigWig files. [Difficulty: **Beginner**]
 
 **solution:**
 ```{r,echo=FALSE,eval=FALSE}
-#coming soon
- 
+library(rtracklayer)
+bam2bw <- function(bam_file){
+  reads <-  readGAlignments(bam_file) # read in bam file
+  reads = granges(reads) # convert to GRanges object
+  reads = resize(reads, width=200, fix='start') # extend each read toward 3' end
+  seqnames <- unique(seqnames(reads)) # get names (once) of each chr in object  
+  reads <- keepSeqlevels(reads, seqnames, pruning.mode='coarse') # remove unused levels
+  cov = coverage(reads, width = seqlengths) # convert reads into signal profile
+  bw_file <- sub('.bam','.bigWig', bam_file) # change filename extension to .bigWig
+  work_dir <- getwd()
+  output_file <- file.path(work_dir, basename(bw_file))
+  export.bw(cov, output_file) # export as bigWig file
+} 
 ```
 
 5. Apply the function to all files, and visualize them in the genome browser.
@@ -39,33 +108,90 @@ Observe the signal profiles. What can you notice, about the similarity of the sa
   
 **solution:**
 ```{r,echo=FALSE,eval=FALSE}
-#coming soon
- 
+# convert bam reads to bigWig coverage profile
+lapply(bam_files, bam2bw)
+# visualize profiles in separate tracks
+library(Gviz)
+axis   = GenomeAxisTrack(
+    range = GRanges('chr21', IRanges(1, width=seqlengths))
+)
+# convert each signal into genomic ranges and define each signal track
+# read in .bigWig files with coverage profiles
+bw_files <- sub('.bam','.bigWig', bam_files) 
+work_dir <- getwd()
+bw_files <- file.path(work_dir, basename(bw_files)) # set path to working dir
+bw_list <- lapply(bw_files, import.bw)
+# name each bigWig object in the list with its experiment name
+expt_names <- sub('.chr21.bam','', basename(bam_files)) # remove path and common suffix
+expt_names = sub('GM12878_hg38_','', expt_names) # remove common prefix
+names(bw_list) <- expt_names
+# visualize each profile in a genome browser
+# dtrack_list <- lapply(bw_list, function(cov){
+#   gcov <- as(cov, 'GRanges')
+#   dtrack <- DataTrack(gcov, name = names(cov), type='l')
+#   return(dtrack)
+# })
+# alternative; b/c when using lapply to produce dtrack_list, plotTracks threw error
+i=1
+for(i in 1:length(bw_list)) {
+  gcov <- as(bw_list[[i]], 'GRanges')
+  dtrack <- DataTrack(gcov, name = names(bw_list[i]), type='l')
+  plotTracks(trackList = list(axis, dtrack), 
+             sizes            = c(.1,1), 
+             background.title = "black")
+  } 
 ```
 
+*All of the Input samples (#7-11) have a low level of reads (peaks rarely exceed 20 counts) that are distributed rather uniformly; again, Input_r5 is the exception. The ChIP samples have more peaks with high levels of reads (often hundreds per peak); SMC3_r1 is an exception.*
+
 6. Use `GViz` to visualize the profiles for CTCF, SMC3 and ZNF143. [Difficulty: **Beginner/Intermediate**]
 
 **solution:**
 ```{r,echo=FALSE,eval=FALSE}
-#coming soon
+# Visualization of profiles are generated by code in previous exercise.
  
 ```
 
-7. Calculate the cross correlation for both CTCF replicates, and
-the input samples. How does the profile look for the control samples? [Difficulty: **Intermediate**]
+7. Calculate the cross correlation for both CTCF replicates, and the input samples. How does the profile look for the control samples? [Difficulty: **Intermediate**]
   
 **solution:**
 ```{r,echo=FALSE,eval=FALSE}
-#coming soon
+# following the example in 9.6.2
+# use tileGenome() to return list of GRanges of given width, spanning whole chromosome
+tilling_window = tileGenome(seqlengths, tilewidth=1000)
+tilling_window = unlist(tilling_window) # convert list to one GRanges object
+so = summarizeOverlaps(tilling_window, bam_files) # count reads in each window
+counts = assays(so)[[1]] # extract counts from SummarizedExperiment
+cpm = t(t(counts)*(1000000/colSums(counts))) # calculate cpm from counts matrix
+cpm = cpm[rowSums(cpm) > 0,] # remove tiles with no reads
+colnames(cpm) = sub('.chr21.bam','', colnames(cpm)) # shorten column names
+colnames(cpm) = sub('GM12878_hg38_','', colnames(cpm))
+# calculate the pearson correlation coefficient between CTCF and input samples
+cor(cpm[,c(1,2,7:11)], method='pearson') # subset for specified columns
  
 ```
 
+*The correlation between the two CTCF replicates is very high (r = 0.94) indicating that the enriched fragments were consistent. As expected, the CTCF replicates do not exhibit strong correlation with any input samples, which do not contain enriched fragments. The input samples exhibit weak to moderate cross correlation (r = 0.39-0.64), consistent with similar composition but weak signals.*
+
 8. Calculate the cross correlation coefficients for all samples and 
 visualize them as a heatmap. [Difficulty: **Intermediate**]
 
 **solution:**
 ```{r,echo=FALSE,eval=FALSE}
-#coming soon
+correlation_matrix <- cor(cpm, method='pearson') # use entire matrix/all samples
+
+library(ComplexHeatmap)
+library(circlize)
+heatmap_col = circlize::colorRamp2( # define color palette where 0 is white
+    breaks = c(-1,0,1),
+    colors = c('blue','white','red')
+)
+
+# plot the heatmap using the Heatmap function
+Heatmap(
+    matrix = correlation_matrix, 
+    col    = heatmap_col
+)
  
 ```
 
@@ -75,7 +201,25 @@ visualize them as a heatmap. [Difficulty: **Intermediate**]
 
 **solution:**
 ```{r,echo=FALSE,eval=FALSE}
-#coming soon
+# rather than arbitrarily select a single Input file (with fewer reads), 
+# let's merge all 5 Input files into one
+library(Rsamtools)
+mergeBam(files = bam_files[7:11],
+         overwrite = TRUE,
+         destination = file.path(work_dir, "GM12878_hg38_Inputs_merged.chr21.bam"),
+         indexDestination = TRUE)
+
+# using code in 9.6.2 as template:
+library(normr)
+# peak calling using CTCF (1,2), SMC3 (12,13), and ZNF143 (14,15) samples with merged Input
+peaks_list <- lapply(bam_files[c(1,2,12:15)], function(chip_file){
+  ctcf_fit = enrichR(treatment = chip_file,
+                     control   = file.path(work_dir, "GM12878_hg38_Inputs_merged.chr21.bam"),
+                     genome    = "hg38",
+                     verbose   = FALSE)
+  })
+# show summary for all samples
+lapply(peaks_list, summary)
  
 ```
 
@@ -83,7 +227,36 @@ visualize them as a heatmap. [Difficulty: **Intermediate**]
 
 **solution:**
 ```{r,echo=FALSE,eval=FALSE}
-#coming soon
+# get regions with peaks for CTCF replicates (1,2)
+CTCF_peaks <- lapply(peaks_list[1:2], function(fit){
+  ctcf_peaks = getRanges(fit) # extract all ranges for an experiment
+  ctcf_peaks$qvalue = getQvalues(fit) # annotate ranges with adjusted p values
+  ctcf_peaks$enrichment = getEnrichment(fit) # annotate ranges with calculated enrichment
+  ctcf_peaks = subset(ctcf_peaks, !is.na(component)) # selects ranges corresponding to enriched class
+  ctcf_peaks = subset(ctcf_peaks, qvalue < 0.01) # filter by stringent q value threshold
+  ctcf_peaks = ctcf_peaks[order(ctcf_peaks$qvalue)] # sort peaks based on q values
+  ctcf_peaks = GenomicRanges::reduce(ctcf_peaks) # collapse neighboring regions
+})
+# clean up levels to avoid problems later
+CTCF_peaks <- lapply(CTCF_peaks, function(sample){
+  keepSeqlevels(sample, 'chr21', pruning.mode='coarse') 
+  })
+
+# find which reads overlap peak regions for each replicate
+# CTCF_r1 sample
+CTCF_r1_reads <- resize(reads_list[[1]], 200) # resize each read assuming 200 bp fragments
+# calculate percentage of reads overlapping a peak
+100*sum(countOverlaps(CTCF_r1_reads, CTCF_peaks[[1]]))/length(CTCF_r1_reads)
+
+# CTCF_r2 sample
+CTCF_r2_reads <- resize(reads_list[[2]], 200) # resize each read assuming 200 bp fragments
+# calculate percentage of reads overlapping a peak
+100*sum(countOverlaps(CTCF_r2_reads, CTCF_peaks[[2]]))/length(CTCF_r2_reads)
+
+# combined for both replicates
+100*(sum(countOverlaps(CTCF_r1_reads, CTCF_peaks[[1]])) + sum(countOverlaps(CTCF_r2_reads, CTCF_peaks[[2]])))/(length(CTCF_r1_reads) + length(CTCF_r2_reads))
+
+*Roughly one-third of the reads in the CTCF_r1 (26%) and CTCF_r2 (35%) map to CTCF peaks.*
  
 ```
 
@@ -101,19 +274,36 @@ How many peaks are specific to each biological replicate, and how many peaks ove
 
 **solution:**
 ```{r,echo=FALSE,eval=FALSE}
-#coming soon
+CTCF_peaks_unified <- CTCF_peaks[[1]][countOverlaps(CTCF_peaks[[1]], CTCF_peaks[[2]]) > 0]
+length(CTCF_peaks_unified)
+length(CTCF_peaks[[1]]) - length(CTCF_peaks_unified)
+length(CTCF_peaks[[2]]) - length(CTCF_peaks_unified)
  
 ```
 
-5. Plot a scatter plot of signal strengths for biological replicates. Do intersecting
-peaks have equal signal strength in both samples? [Difficulty: **Intermediate**]
+*Altogether 429 peaks overlap between the two CTCF replicates, while 47 peaks are specific to r1 and 286 are specific to r2.*
+
+5. Plot a scatter plot of signal strengths for biological replicates. Do intersecting peaks have equal signal strength in both samples? [Difficulty: **Intermediate**]
 
 **solution:**
 ```{r,echo=FALSE,eval=FALSE}
-#coming soon
+# for the 429 peaks shared by the CTCF replicates, count # reads per peak in each replicate
+reads_CTCFpeak <- tibble(CTCF_r1 = countOverlaps(CTCF_peaks_unified, CTCF_r1_reads),
+                         CTCF_r2 = countOverlaps(CTCF_peaks_unified, CTCF_r2_reads))
+# normalize to account for difference in total number of reads in each sample
+reads_CTCFpeak <- reads_CTCFpeak %>% 
+  mutate(cpm_r1 = 10^6*CTCF_r1/length(CTCF_r1_reads),
+         cpm_r2 = 10^6*CTCF_r2/length(CTCF_r2_reads))
+# plot the normalized count for each CTCF peak
+reads_CTCFpeak %>% 
+  ggplot(aes(x=cpm_r1, y=cpm_r2)) +
+  geom_point() +
+  geom_abline(slope = 1, color = "red")
  
 ```
 
+*After normalizing the counts, the signal strengths of the intersecting peaks are generally comparable in the two samples, as indicated by the close fit with the red line for y=x. However, there does seems to be some excess of normalized read counts in r2 vs. r1.*
+
 6. Quantify the combinatorial binding of all three proteins. Find
 the number of places which are bound by all three proteins, by 
 a combination of two proteins, and exclusively by one protein.
@@ -152,10 +342,39 @@ How many motifs do you observe? How do the motifs look (visualize the motif logs
 
 **solution:**
 ```{r,echo=FALSE,eval=FALSE}
-#coming soon
+ZNF143_peaks <- lapply(peaks_list[5:6], function(fit){
+  peaks = getRanges(fit) # extract all ranges for an experiment
+  peaks$qvalue = getQvalues(fit) # annotate ranges with adjusted p values
+  peaks$enrichment = getEnrichment(fit) # annotate ranges with calculated enrichment
+  peaks = subset(peaks, !is.na(component)) # selects ranges corresponding to enriched class
+  peaks = subset(peaks, qvalue < 0.01) # filter by stringent q value threshold
+  peaks = peaks[order(peaks$qvalue)] # sort peaks based on q values
+  peaks = GenomicRanges::reduce(peaks) # collapse neighboring regions
+})
+# clean up levels to avoid problems later
+ZNF143_peaks <- lapply(ZNF143_peaks, function(sample){
+  keepSeqlevels(sample, 'chr21', pruning.mode='coarse') 
+  })
+# focus each peak at its center
+ZNF143_peaks_resized <- lapply(ZNF143_peaks, function(peak){
+  resize(peak, width = 50, fix='center')
+  })
+library(rGADEM)
+# load the human genome sequence
+library(BSgenome.Hsapiens.UCSC.hg38)
+# run GADEM() for motif discovery 
+novel_motifs <- GADEM(unlist(as(ZNF143_peaks_resized, "GRangesList")), # use peaks of both replicates
+                    verbose=1, # print results to screen
+                    genome=Hsapiens) # need genome to get base composition
+consensus(novel_motifs) # view the consensus sequence of each motif
+nOccurrences(novel_motifs) # count occurrences of each motif
+getPWM(novel_motifs) # get position weight matrix for each motif
+plot(novel_motifs) # visualize each motif
  
 ```
 
+*Three motifs are observed.*
+
 2. Scan the ZNF143 peaks with the top motifs found in the previous exercise. 
 Where are the motifs located? [Difficulty: **Advanced**]