quiz-1.Rmd

---
title: "Quiz 1"
author: "Muhammad Abdullah Nabeel"
date: "7/6/2024"
output: html_document
---

```{r setup, include=FALSE}
#knitr::opts_chunk$set(echo = TRUE)
```

# Quiz 1

Useful links for initial questions before we start:

-   [https://www.coursera.org/learn/bioconductor/discussions/forums/ux-4Iij3Eea8jw6UvTi2Tw/threads/YLZyg9f_EeW17wrgkMLDIQhttps://www.coursera.org/learn/bioconductor/discussions/forums/ux-4Iij3Eea8jw6UvTi2Tw/threads/YLZyg9f_EeW17wrgkMLDIQ](https://www.coursera.org/learn/bioconductor/discussions/forums/ux-4Iij3Eea8jw6UvTi2Tw/threads/YLZyg9f_EeW17wrgkMLDIQ){.uri}

-   [Bioconductor for Genomic Data Science - Discussions \| Coursera](https://www.coursera.org/learn/bioconductor/discussions/forums/ux-4Iij3Eea8jw6UvTi2Tw/threads/JW10d3MTRd6tdHdzE_Xe0Q)

session direcory

```{r}
setwd("D:/CB/biocond/coursera-gds-bc")
```

## Q1.

Use the AnnotationHub package to obtain data on "CpG Islands" in the human genome.

**Question:** How many islands exists on the autosomes?

1 point

**26641**

**25557**

**22215**

**13810**

```{r}
# BiocManager::version()
# BiocManager::install("AnnotationHub")
# BiocManager::install("GenomicRanges")
library("AnnotationHub")
library(GenomicRanges)

ah = AnnotationHub()
ah

ha = subset(ah, species == "Homo sapiens")

ha

q1_data = query(ha, "CpG Islands")

q1_data

q1_datum = q1_data[[1]]
q1_datum

# we have right data
# now check how many ranges are on autosomes
# autosomes are chromosomes other than sex chromosome
# removing non-standad chromosomes/contigs

q1_std_chr = keepStandardChromosomes(q1_datum, pruning.mode = "coarse")


q1_std_chr
seqlevels(q1_std_chr)

q1_autosomes = dropSeqlevels(q1_std_chr, c("chrX", "chrY", "chrM"), pruning.mode="coarse")

q1_autosomes

#verify you only have autosomes
seqlevels(q1_autosomes)
q1_autosomes

ranges(q1_autosomes)
length(q1_autosomes) # this gives you the answer
#26641
```

## Q2.

Question: How many CpG Islands exists on chromosome 4.

1 point

1163

1688

1031

801

```{r}
# using standard chromosome ranes from q1
q1_std_chr
q2_chr4 = keepSeqlevels(q1_std_chr, "chr4", pruning.mode = "coarse")

q2_chr4

# answer
length(q2_chr4)
# 1031
```

## Q3.

#### Hint: Since choosing the right dataset is sifficult due to unclear instructions. Must see the first two links at top to know what we are going to search.

3.  Question 3 Obtain the data for the H3K4me3 histone modification for the H1 cell line from Epigenomics Roadmap, using AnnotationHub. Subset these regions to only keep regions mapped to the autosomes (chromosomes 1 to 22).

Question: How many bases does these regions cover?

1 point

37029137

40252794

37923727

41135164

```{r}
ah = AnnotationHub()
ah = subset(ah, species == "Homo sapiens")
ah3 = query(ah, c("H3K4me3", "EpigenomeRoadMap", "E003"))

ah3
ah3$species
mcols(ah3)$genome

ah3[1]
ah3[2]

gr1 = ah3[[1]]
gr2 = ah3[[2]]

summary(width(gr1))
summary(width(gr2))
# it seems like there is more dispersion in gr2
# dispersion is almost alike form minimum to medium/mean
# after mean gr1 is more dispersed and gr2 is less
# but this is still my guess, as I see the summary, I didnt run any # statistical tests

# choosing less dispersed gr2
seqlevels(gr2)
ah3_std = keepStandardChromosomes(gr2, pruning.mode = "coarse")
seqlevels(ah3_std)
ah3_auto = dropSeqlevels(ah3_std, c("chrX", "chrY", "chrM"), pruning.mode = "coarse")
seqlevels(ah3_auto)

seqlengths(ah3_auto)

# following two statements give same and they answer to the question
sum(width(reduce(ah3_auto)))
sum(width(ah3_auto))

# 41135164
```

## Q4.

Question 4

Obtain the data for the H3K27me3 histone modification for the H1 cell line from Epigenomics Roadmap, using the AnnotationHub package. Subset these regions to only keep regions mapped to the autosomes. In the return data, each region has an associated "signalValue".

**Question:** What is the mean signalValue across all regions on the standard chromosomes?

1 point

**4.770728**

**4.799766**

**4.372933**

**5.140497**

```{r}
ah = AnnotationHub()
ah = subset(ah, species == "Homo sapiens")
ah4 = query(ah, c("H3K27me3", "EpigenomeRoadMap", "E003"))

ah4
ah4$species
mcols(ah4)$genome

ah4[1]
ah4[2]

gr1 = ah4[[1]]
gr2 = ah4[[2]]

summary(width(gr1))
summary(width(gr2))

qqnorm(width(gr1))
qqnorm(width(gr2))

# answers
mean(gr1$signalValue)
mean(gr2$signalValue)

# for use in next question
seqlevels(gr2)
ah4_std = keepStandardChromosomes(gr2, pruning.mode = "coarse")
seqlevels(ah4_std)
ah4_auto = dropSeqlevels(ah4_std, c("chrX", "chrY", "chrM"), pruning.mode = "coarse")
seqlevels(ah4_auto)

seqlengths(ah4_auto)
# 4.771152
```

## Q5.

Question 5

Bivalent regions are bound by both H3K4me3 and H3K27me3.

**Question:** Using the regions we have obtained above, how many bases on the standard chromosomes are bivalently marked?

1 point

**10207246**

**10289096**

**10984729**

**9926503**

```{r}
ah3_auto
ah4_auto
seqlevels(ah4_auto) == seqlevels(ah3_auto)

bivalent = intersect(ah3_auto, ah4_auto)
bivalent

# bivalent = subsetByOverlaps(ah3_auto, ah4_auto)
# bivalent

sum(width(reduce(bivalent)))
# answer is 10289096
```

## Q6-7.

We will examine the extent to which bivalent regions overlap CpG Islands.

**Question:** how big a fraction (expressed as a number between 0 and 1) of the bivalent regions, overlap one or more CpG Islands?

1 point

**0.5140894**

**0.5839835**

**0.5383644**

**0.5059474**

**Question:** How big a fraction (expressed as a number between 0 and 1) of the bases which are part of CpG Islands, are also bivalent marked

1 point

**0.241688**

**0.2924248**

**0.1860021**

**0.2750512**

```{r}
q1_autosomes
bivalent
ov6 = findOverlaps(bivalent, q1_autosomes)
ov6
length(unique(queryHits(ov6)))

# Q6
length(unique(queryHits(ov6))) / length(bivalent)
# 0.5383644

# Q7
length(unique(subjectHits(ov6))) / length(q1_autosomes)
# 0.2915806 Choose nearest answer
```

## Q8.

**Question:** How many bases are bivalently marked within 10kb of CpG Islands?

**Tip:** consider using the "resize()"" function.

1 point

**11928130**

**11104114**

**9782086**

**7920203**

```{r}
# resize q1_autosomes
q1_autosomes
big_islands <- resize(q1_autosomes, width = (20000 + width(q1_autosomes)), fix = "center")
x = intersect(bivalent, big_islands)
x
sum(width(reduce(x)))
# 9782086
```

## Q9.

Question 9

**Question:** How big a fraction (expressed as a number between 0 and 1) of the human genome is contained in a CpG Island?

**Tip 1:** the object returned by AnnotationHub contains "seqlengths".

**Tip 2:** you may encounter an integer overflow. As described in the session on R Basic Types, you can address this by converting integers to numeric before summing them, "as.numeric()".

1 point

**0.006823921**

**0.007288502**

**0.005769540**

**0.007047481**

```{r}
# denominator is the length of human genome
genome = q1_autosomes
seqlengths(genome)
as.numeric(seqlengths(genome))

genome.length = sum(as.numeric(seqlengths(genome)))
genome.length

sum(as.numeric(width(q1_autosomes)))

sum(as.numeric(width(q1_autosomes))) / genome.length
# 0.007047481 is the answer
```

## Q10.

Question 10

**Question:** Compute an odds-ratio for the overlap of bivalent marks with CpG islands.

1 point

**193.1999**

**168.2223**

**169.0962**

**211.6532**

```{r}
inOut = matrix(0, ncol = 2, nrow = 2)
colnames(inOut) = c("in", "out")
rownames(inOut) = c("in", "out")

#inOut

inOut[1,1] = sum(width(intersect(
  bivalent, q1_autosomes, ignore.strand=TRUE)))
inOut[1,2] = sum(width(setdiff(
  bivalent, q1_autosomes, ignore.strand=TRUE)))
inOut[2,1] = sum(width(setdiff(
  q1_autosomes, bivalent, ignore.strand=TRUE)))
inOut[2,2] = genome.length - sum(inOut)

inOut

odd_ratio <- inOut[1,1]*inOut[2,2]/(inOut[1,2]*inOut[2,1])

odd_ratio
```