ProjectRNAseqCelegans/Week 8 code at main · olinm/ProjectRNAseqCelegans · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
#WEEK 8

---
TITLE: "Oliwia's RNASeq notebook: part 5"
---

Exploring different analysis methods for the project, merging and comparing acquired data with other sources.

---
RStudio info
---
platform       x86_64-w64-mingw32
arch           x86_64
os             mingw32
system         x86_64, mingw32
status
major          4
minor          1.3
year           2022
month          03
day            10
svn rev        81868
language       R
version.string R version 4.1.3 (2022-03-10)
nickname       One Push-Up

---
PART 1: POSITION COUNT MATRICES FOR MOTIF VISUALISATION
---
References:
http://jvanheld.github.io/cisreg_course/from_sites_to_matrices/from_TFBS_to_PSSM_solutions.html
https://bioconductor.org/packages/release/bioc/vignettes/universalmotif/inst/doc/IntroductionToSequenceMotifs.pdf
https://bioconductor.org/packages/release/bioc/vignettes/universalmotif/inst/doc/MotifManipulation.pdf
http://biocworkshops2019.bioconductor.org.s3-website-us-east-1.amazonaws.com/page/motifStackWorkshop__motifStackWorkshop/

#installing the required package
```{r}
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("motifStack")
```

```{r}
library(motifStack)
library(knitr)
```


#HIF1A
data preparation
```{r eval=FALSE, include=FALSE}
hif300hits <- head(hif300hits, - 3)
```

```{r}
sequenceshif <- hif300hits$matched_sequence
print(as.data.frame(sequenceshif))
```

#alignment matric - one row per site, one column per aligned position
```{r}
alignment <- t(as.data.frame(strsplit(sequenceshif, "")))
colnames(alignment) <- 1:ncol(alignment)
row.names(alignment) <- paste(sep="_", "site", 1:nrow(alignment))

kable(alignment, align="c")
```

#PCM (position count matrix)
```{r}
countshif <- data.frame()
residues <- c("A", "C", "G", "T")
for (r in residues) {
  countshif <- rbind(countshif, apply(alignment==r, 2, sum))
}
colnames(countshif) <- 1:ncol(countshif)
row.names(countshif) <- residues

kable(countshif, align="c")
```

#PPM (position probability matrix) - raw counts
```{r}
freq <- countshif / apply(countshif, 2, sum)
kable(freq, align="c")
```

#residue frequencies set to 0.25
```{r}
equal.priors <- vector()
equal.priors[residues] <- 1/length(residues)

kable(data.frame(equal.priors), align="c")
```

#adding pseudocounts
How to decide?
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2647310/

As suggested by "Pseudocounts for transcription factor binding sites', I will use a count of 0.8.

```{r}
k <- 0.8
corrected.freq <- (countshif + k/length(residues)) / apply(countshif + k*0.25, 2, sum)
```

#visualising the results
```{r}
plotMotifLogo(corrected.freq)
```

#can be modified using each TF data to construct all the visualisations

---
SESSION INFO
---
R version 4.1.3 (2022-03-10)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19043)

Matrix products: default

locale:
[1] LC_COLLATE=English_Europe.1252  LC_CTYPE=English_Europe.1252    LC_MONETARY=English_Europe.1252
[4] LC_NUMERIC=C                    LC_TIME=English_Europe.1252

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

loaded via a namespace (and not attached):
 [1] MatrixGenerics_1.6.0        Biobase_2.54.0              httr_1.4.2
 [4] bit64_4.0.5                 splines_4.1.3               assertthat_0.2.1
 [7] BiocManager_1.30.16         stats4_4.1.3                blob_1.2.2
[10] GenomeInfoDbData_1.2.7      yaml_2.3.5                  pillar_1.7.0
[13] RSQLite_2.2.10              lattice_0.20-45             glue_1.6.1
[16] digest_0.6.29               GenomicRanges_1.46.1        RColorBrewer_1.1-2
[19] XVector_0.34.0              colorspace_2.0-3            htmltools_0.5.2
[22] Matrix_1.4-0                DESeq2_1.34.0               XML_3.99-0.8
[25] pkgconfig_2.0.3             genefilter_1.76.0           zlibbioc_1.40.0
[28] purrr_0.3.4                 xtable_1.8-4                scales_1.1.1
[31] tzdb_0.2.0                  BiocParallel_1.28.3         tibble_3.1.6
[34] annotate_1.72.0             KEGGREST_1.34.0             generics_0.1.2
[37] IRanges_2.28.0              ggplot2_3.3.5               ellipsis_0.3.2
[40] cachem_1.0.6                SummarizedExperiment_1.24.0 BiocGenerics_0.40.0
[43] cli_3.2.0                   survival_3.2-13             magrittr_2.0.2
[46] crayon_1.5.0                memoise_2.0.1               evaluate_0.15
[49] fansi_1.0.2                 tools_4.1.3                 hms_1.1.1
[52] lifecycle_1.0.1             matrixStats_0.61.0          S4Vectors_0.32.3
[55] munsell_0.5.0               locfit_1.5-9.4              DelayedArray_0.20.0
[58] AnnotationDbi_1.56.2        Biostrings_2.62.0           compiler_4.1.3
[61] GenomeInfoDb_1.30.1         rlang_1.0.1                 grid_4.1.3
[64] RCurl_1.98-1.6              rstudioapi_0.13             bitops_1.0-7
[67] rmarkdown_2.13              gtable_0.3.0                DBI_1.1.2
[70] R6_2.5.1                    knitr_1.37                  dplyr_1.0.8
[73] fastmap_1.1.0               bit_4.0.4                   utf8_1.2.2
[76] readr_2.1.2                 parallel_4.1.3              Rcpp_1.0.8
[79] vctrs_0.3.8                 geneplotter_1.72.0          png_0.1-7
[82] tidyselect_1.1.2            xfun_0.29