forked from kapelsinas/biomed-data-analysis-jdk
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Task_2_JDK.Rmd
191 lines (167 loc) · 8.59 KB
/
Task_2_JDK.Rmd
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
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
---
title: "Practical No. 2"
author: "Jaroslav Rutkovskij, Danielė Stasiūnaitė, Karolis Augustauskas (JDK)"
date: "3/18/2022"
output: pdf_document
---
```{r Preparations, message=FALSE, warning=FALSE, include=FALSE}
options(scipen = 1, digits = 3)
library(minfi)
library(IlluminaHumanMethylationEPICmanifest)
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("IlluminaHumanMethylationEPICanno.ilm10b4.hg19")
```
```{r Preparations_Task, message=FALSE, warning=FALSE, include=FALSE}
# To load the data
# load("objects.Rdata")
# To save the data
# save.image(file = "objects.Rdata")
```
```{r Task_1, message=FALSE, warning=FALSE, include=FALSE}
################################################################################################
# 1. Pasiruoškite sau paskirtus sample-key failus (apvalydami reikšmes, panašiai kaip aptarėme
# paskaitų metu)
# - sample-key nuskaitymui R galite naudoti read.csv() funkciją.
# Importing data:
################################################################################################
if (file.exists("data.csv")) {
data <- read.csv("data.csv")
}
```
```{r Task_2, message=FALSE, warning=FALSE, include=FALSE}
################################################################################################
# 2. Nuskaitykite parsisiųstus idat failus į "RGChannelSet" tipo objektą
# - tam naudokite `read.metharray.exp()` funkcią iš R minfi bibliotekos
# - ši funkcija kaip argumentą naudoja "sample-key" lentelę, kurioje stulpelis Basename
# turi rodyti į failo vietą jūsų kompiuteryje (galite paskaityti šios funkcijos help()).
# - nuskaičius duomenis gautą objektą išsaugokite su saveRDS()
# - vėliau jį pasikrauti į atmintį galėsite su readRDS() funkcija
################################################################################################
rgSet <- read.metharray.exp(base = "IDATS")
save(rgSet, file = "rgSet0.RData")
```
```{r Task_3, message=FALSE, warning=FALSE, include=FALSE}
################################################################################################
# 3. Gaukite "detection p-value" kiekvienam DNR modifikacijos įverčiui
# - naudokite funkciją `detectionP()` iš minfi
# - ši funkcija grąžins įverčių matricą (vienas ivertis vienai modifikacijos reiksmei)
# parodanti ar matuojamas intensyvumas skiriasi nuo background
################################################################################################
detP <- detectionP(rgSet)
save(detP, file = "detP0.Rdata")
NROW(rgSet) #
# - visas vertes kuriu p-verte didesne uz 0.01 laikykite "blogomis"
badPValues <- colMeans(detP) > 0.01
keep <- colMeans(detP) < 0.01
# Nieko neišsimetėme
# - is RGChannelSet objekto ismeskite visus meginius (stulpelius) kurie turi daugiau
# nei 1% "blogu" detection p reiksmiu.
# - PAGALBA: is RGChannelSet objekto reiksmes ismesti galima taip, kaip ir is
# matricos: `myData[,-1]` - ismestu pirma megini (stulpeli)
rgSet <- rgSet[, keep]
data <- data[keep,]
detP <- detP[, keep]
NROW(rgSet) #
save(rgSet, file = "rgSet1.RData")
save(detP, file = "detP1.Rdata")
save(data, file ="data0.RData")
```
```{r Task_4, message=FALSE, warning=FALSE, include=FALSE}
####################################################################################################
# 4. Normalizuokite savo duomenis
# - naudokite vieną iš šių funkcijų: preprocessSWAN(), preprocessFunnorm() arba
# preproccessIllumina() iš "minfi" bibliotekos.
# - galite pasirinkti laisvai, apie tai kuo skiriasi šie normalizavimo būdai galite
# pasiskaityti šių funkcijų help() dokumentacijoje.
####################################################################################################
mSetSq <- preprocessFunnorm(rgSet) # Buvo naudotą mūsų straipsnyje
save(mSetSq, file = "mSetSq0.RData")
mSetRaw <- preprocessRaw(rgSet)
save(mSetRaw, file = "mSetRaw0.RData")
```
```{r Task_7, message=FALSE, warning=FALSE, include=FALSE}
####################################################################################################
# 7. Pašalinkite mėginius, kurių nurodyta lytis skiriasi nuo spėjamos lyties (spėjimus galite gauti
# su funkcija getSex()).
####################################################################################################
estSex <- getSex(mSetSq, cutoff = -2)
estSex$predictedSex == 'F'
# Vien moterų mėginiai
rm(estSex)
```
```{r Task_5, message=FALSE, warning=FALSE, include=FALSE}
####################################################################################################
# 5. Išmeskite visas genomines pozicijas (eilutes) kurios turi daugiau nei 1% "blogų" detection
# p reikšmių
# - detection p-reikšmės paskaičiuojamos prieš normalizaciją, taigi naudokite prieš tai
# paskaičiuotą p-reiksmių matricą.
# - jeigu iš nuskaityto objekto buvo pašalinta pavyzdžių - nepamirškite prieš atliekant šį
# žingsnį jų pašalinti ir iš detection p-reiksmių matricos.
# - PAGALBA: normalizavimo žingsnis gali pašalinti kai kurias genomines pozicijas (eilutes)
# iš nuskaityto objekto, todėl naudokite
# `rownames()`, kad sužinotumėte, kurios pozicijos turi būti pašalintos.
# ensure probes are in the same order in the mSetSq and detP objects
####################################################################################################
# ensure probes are in the same order in the mSetSq and detP objects
detP <- detP[match(featureNames(mSetSq), rownames(detP)), ]
save(detP, file = "detP1.Rdata")
# remove any probes that have failed in one or more samples
keep <- rowSums(detP < 0.01) == ncol(mSetSq)
NROW(mSetSq) #
table(keep)
## keep
## FALSE TRUE
## 9713 856146
mSetSq <- mSetSq[keep, ]
NROW(mSetSq) #
rm(detP)
```
```{r Task_6, message=FALSE, warning=FALSE, include=FALSE}
####################################################################################################
# 6. Išmeskite genomines pozicijas (eilutes) neturinčias "CG" nukleotidų poros (CH) arba esančias
# šalia DNR polimorfizmų
# - tam naudokite funkciją dropMethylationLoci(), galite perskaityti jos help().
####################################################################################################
NROW(mSetSq) # 856146
mSetSq <- dropMethylationLoci(mSetSq, dropCH = TRUE)
NROW(mSetSq) # 853368
```
```{r Task_8, message=FALSE, warning=FALSE, include=FALSE}
####################################################################################################
# 8. Po duomenų paruošimo iš gauto objekto pasidarykite tris atskirus objektus:
# - pagrindinė modifikacijos įverčių matrica: getBeta()
# - informacija apie pagrindinės matricos mėginius (stulpelius): pData()
# - informacija apie pagrindinės matricos pozicijas (eilutes): getAnnotation()
####################################################################################################
beta <- getBeta(mSetSq)
pData <- pData(mSetSq)
annotation <- getAnnotation(mSetSq)
save(beta, file = "beta.Rdata")
save(pData, file = "pData.Rdata")
save(annotation, file = "annotation.Rdata")
```
```{r Task_9, message=FALSE, warning=FALSE, include=FALSE}
####################################################################################################
# 9. Savarankiškai parašykite kodą, atliekantį IAC išskirčių išmetimą.
# - Detalus sios proceduros aprasymas, net su kodo pavyzdziais, gali buti rastas cia:
# https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/HumanBrainTranscriptome/Identification%20and%20Removal%20of%20Outlier%20Samples.pdf)
####################################################################################################
```
```{r Task_10, message=FALSE, warning=FALSE, include=FALSE}
####################################################################################################
# 10. Atlikite bent vieną kokybės kontrolės žingsnį aptartą paskaitų metu, arba susigalvotą savo
# paties (žinoma galite atlikti ir daugiau nei vieną).
####################################################################################################
load("rgSet1.RData")
qcReport(rgSet, pdf= "qcReport.pdf")
rm(rgSet)
```
```{r Task_11, message=FALSE, warning=FALSE, include=FALSE}
####################################################################################################
# 11. Išsisaugokite paruoštus duomenis (su saveRDS()), juos naudosite kitose užduotyse
####################################################################################################
# load("beta.Rdata")
# load("pData.Rdata")
# load("annotation.Rdata")
```