-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathDGE_analysis_SARSCoV2_med.R
257 lines (219 loc) · 11.3 KB
/
DGE_analysis_SARSCoV2_med.R
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
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
#Hundreds of thousands of ribo counts
library(purrr)
library(dplyr)
library(edgeR)
library(limma)
library(DESeq2)
library(stringr)
library(Glimma)
library(sva)
setwd("D://Projects/kenezam/gene_counts/SARSCOV2/")
df1 <- read.csv("PRJNA704763_ribo_mock_novirus_24h_1_5s/ReadsPerGene.out.tab", sep="\t")[-c(1:3),c(1,3)]
df2 <- read.csv("PRJNA704763_ribo_mock_novirus_24h_2_5s/ReadsPerGene.out.tab", sep="\t")[-c(1:3),c(1,3)]
df3 <- read.csv("PRJNA704763_ribo_mock_novirus_24h_3_5s/ReadsPerGene.out.tab", sep="\t")[-c(1:3),c(1,3)]
df4 <- read.csv("PRJNA704763_ribo_virus_SARSCOV2_6h_1_5s/ReadsPerGene.out.tab", sep="\t")[-c(1:3),c(1,3)]
df5 <- read.csv("PRJNA704763_ribo_virus_SARSCOV2_6h_2_5s/ReadsPerGene.out.tab", sep="\t")[-c(1:3),c(1,3)]
df6 <- read.csv("PRJNA704763_ribo_virus_SARSCOV2_6h_3_5s/ReadsPerGene.out.tab", sep="\t")[-c(1:3),c(1,3)]
df7 <- read.csv("PRJNA704763_ribo_virus_SARSCOV2_24h_1_5s/ReadsPerGene.out.tab", sep="\t")[-c(1:3),c(1,3)]
df8 <- read.csv("PRJNA704763_ribo_virus_SARSCOV2_24h_2_5s/ReadsPerGene.out.tab", sep="\t")[-c(1:3),c(1,3)]
df9 <- read.csv("PRJNA704763_ribo_virus_SARSCOV2_24h_3_5s/ReadsPerGene.out.tab", sep="\t")[-c(1:3),c(1,3)]
dfr1 <- read.csv("PRJNA704763_rna_mock_novirus_24h_1_4s/ReadsPerGene.out.tab", sep="\t")[-c(1:3),c(1,3)]
dfr2 <- read.csv("PRJNA704763_rna_mock_novirus_24h_2_4s/ReadsPerGene.out.tab", sep="\t")[-c(1:3),c(1,3)]
dfr3 <- read.csv("PRJNA704763_rna_mock_novirus_24h_3_4s/ReadsPerGene.out.tab", sep="\t")[-c(1:3),c(1,3)]
dfr4 <- read.csv("PRJNA704763_rna_virus_SARSCOV2_6h_1_4s/ReadsPerGene.out.tab", sep="\t")[-c(1:3),c(1,3)]
dfr5 <- read.csv("PRJNA704763_rna_virus_SARSCOV2_6h_2_4s/ReadsPerGene.out.tab", sep="\t")[-c(1:3),c(1,3)]
dfr6 <- read.csv("PRJNA704763_rna_virus_SARSCOV2_6h_3_4s/ReadsPerGene.out.tab", sep="\t")[-c(1:3),c(1,3)]
dfr7 <- read.csv("PRJNA704763_rna_virus_SARSCOV2_24h_1_4s/ReadsPerGene.out.tab", sep="\t")[-c(1:3),c(1,3)]
dfr8 <- read.csv("PRJNA704763_rna_virus_SARSCOV2_24h_2_4s/ReadsPerGene.out.tab", sep="\t")[-c(1:3),c(1,3)]
dfr9 <- read.csv("PRJNA704763_rna_virus_SARSCOV2_24h_3_4s/ReadsPerGene.out.tab", sep="\t")[-c(1:3),c(1,3)]
cts <- purrr::reduce(
#list(df1, df2, df3, df4, df5, df6, df7, df8, df9, df10, dfr1, dfr2, dfr3, dfr4, dfr5, dfr6, dfr7, dfr8, dfr9, dfr10),
list(df1, df2, df3, df4, df5, df6, df7, df8, df9, dfr1, dfr2, dfr3, dfr4, dfr5, dfr6, dfr7, dfr8, dfr9),
full_join,
by = "N_unmapped")
#Not actual while using filterByExpr
#cts <- cts[rowSums(dplyr::select(cts,where(is.numeric)))>0,]
id = gsub("\\..*","",cts$N_unmapped)
cts <- cbind(id, cts)
library("AnnotationDbi")
library("org.Hs.eg.db")
cts$id <- AnnotationDbi::mapIds(x = org.Hs.eg.db,
keys = cts$id,
column = "SYMBOL",
keytype = "ENSEMBL",
multiVals = "list")
colnames(cts) <- c("GENES_SYMBOL",
"ENSEMBL_ID",
"ribo_mock_novirus_24h_1",
"ribo_mock_novirus_24h_2",
"ribo_mock_novirus_24h_3",
"ribo_virus_SARS-CoV-2_6h_1",
"ribo_virus_SARS-CoV-2_6h_2",
"ribo_virus_SARS-CoV-2_6h_3",
"ribo_virus_SARS-CoV-2_24h_1",
"ribo_virus_SARS-CoV-2_24h_2",
"ribo_virus_SARS-CoV-2_24h_3",
"rna_mock_novirus_24h_1",
"rna_mock_novirus_24h_2",
"rna_mock_novirus_24h_3",
"rna_virus_SARS-CoV-2_6h_1",
"rna_virus_SARS-CoV-2_6h_2",
"rna_virus_SARS-CoV-2_6h_3",
"rna_virus_SARS-CoV-2_24h_1",
"rna_virus_SARS-CoV-2_24h_2",
"rna_virus_SARS-CoV-2_24h_3")
cts$GENES_SYMBOL <- gsub('[c()"",]', "", cts$GENES_SYMBOL)
cts_rna <- cbind(cts[,c(1:2)],dplyr::select(cts, contains("rna_")))
cts_ribo <- cbind(cts[,c(1:2)],dplyr::select(cts, contains("ribo_")))
time <- c(rep(c(rep(24, 3), rep(6, 3), rep(24, 3)),2))
time <- factor(time)
type <- c(rep('mock', 3), rep('sarscov2', 6), rep('mock', 3), rep('sarscov2', 6))
type <- factor(type)
reg <- c(rep('ribo', 9), rep('rna', 9))
reg <- factor(reg)
rep <- c(rep(c('1', '2', '3'), 6))
rep <- factor(rep)
type <- relevel(type, ref = 'mock')
reg <- relevel(reg, ref = 'rna')
group = reg:type
ann <- data.frame(type=type,reg=reg,rep=rep,time=time)
x <- cts[,c(3:ncol(cts))]
rownames(x) <- paste(cts$GENES_SYMBOL, cts$ENSEMBL_ID, sep="__")
y <- DGEList(counts=x,group=group)
keep <- filterByExpr(y, group=group)
y <- y[keep,,keep.lib.sizes=FALSE]
#glMDSPlot(y[,1:6], groups=group, folder = "glimma-plots/")
glimmaMDS(y[,1:9], groups=ann[1:9,], gene.selection = "common", folder="mds")
#glMDSPlot(y[,7:12], groups=group, folder = "glimma-plots/")
glimmaMDS(y[,10:18], groups=ann[10:18,], gene.selection = "common", folder="mds")
y <- calcNormFactors(y)
plot(density(log2(y$counts[,1]*1000000/(y$samples$lib.size[1]+1))), from = -2, to = 20, lwd = 3)
for (k in 2:9){
lines(density(log2(y$counts[,k]*1000000/(y$samples$lib.size[k]+1))),col=k)
}
plot(density(log2(y$counts[,7]*1000000/(y$samples$norm.factors[1]*y$samples$lib.size[7]+1))), from = -2, to = 20, lwd = 3)
for (k in 10:18){
lines(density(log2(y$counts[,k]*1000000/(y$samples$norm.factors[k]*y$samples$lib.size[k]+1))),col=k)
}
design <- model.matrix(~0 + group + group:time)
design_c <- design[,-c(5,7,9,11,13,15)]
ncol(design_c) <= qr(design_c)$rank
y <- estimateDisp(y,design_c)
fit <- glmQLFit(y,design_c)
#grouprna:wt grouprna:ko groupribo:wt groupribo:ko
#1. 4
qlf_RO_sarscov2_15_4 <- glmQLFTest(fit, contrast = c(1, -1, -1, 1, rep(0,10)))
qlf_rna_sarscov2_15_4 <- glmQLFTest(fit, contrast = c(c(-1, 1), rep(0, 12)))
qlf_ribo_sarscov2_15_4 <- glmQLFTest(fit, contrast = c(rep(0,2), c(-1, 1), rep(0,10)))
res_RO_sarscov2_15_4 <- topTags(qlf_RO_sarscov2_4, n = Inf)[[1]]
res_rna_sarscov2_15_4 <- topTags(qlf_rna_sarscov2_4, n =Inf)[[1]]
res_ribo_sarscov2_15_4 <- topTags(qlf_ribo_sarscov2_4, n =Inf)[[1]]
#Another join
x$row_names <- rownames(x)
library("IHW")
res_RO_sarscov2_4$row_names <- rownames(res_RO_sarscov2_4)
res_RO_sarscov2_4 <- left_join(res_RO_sarscov2_4, x, by="row_names")
#IHW
ihwRes <- ihw(pvalues = res_RO_sarscov2_4$PValue, covariates = rowMeans(dplyr::select(res_RO_sarscov2_4, contains("virus"))), alpha = 0.1)
adj_pvalues(ihwRes)
res_RO_sarscov2_4$adj_palues <- adj_pvalues(ihwRes)
count(adj_pvalues(ihwRes) < 0.05)
#BH
FDR <- p.adjust(res_RO_sarscov2_4$PValue, method="BH")
sum(FDR < 0.1)
res_ribo_48$row_names <- rownames(res_ribo_48)
res_ribo_48 <- left_join(res_ribo_48, x, by="row_names")
#IHW
ihwRes <- ihw(pvalues = res_ribo_48$PValue, covariates = rowMeans(dplyr::select(res_ribo_48, contains("ribo_"))), alpha = 0.1)
adj_pvalues(ihwRes)
res_ribo_48$adj_palues <- adj_pvalues(ihwRes)
count(adj_pvalues(ihwRes) < 0.05)
#BH
FDR <- p.adjust(res_ribo_48$PValue, method="BH")
sum(FDR < 0.1)
#2. 24
qlf_RO_sarscov2_24 <- glmQLFTest(fit, contrast = c(1, -1, -1, 1, -1, 1, rep(0, 4), 0.217, -0.217, -0.217, 0.217))
qlf_rna_sarscov2_24 <- glmQLFTest(fit, contrast = c(c(-1, 1), rep(0, 2), 1, rep(0, 5), c(-0.217, 0.217), rep(0,2)))
qlf_ribo_sarscov2_24 <- glmQLFTest(fit, contrast = c(rep(0, 2), c(-1, 1), 0, 1, rep(0, 4), rep(0,2), c(-0.217, 0.217)))
res_RO_sarscov2_24 <- topTags(qlf_RO_sarscov2_24, n = Inf)[[1]]
res_rna_sarscov2_24 <- topTags(qlf_rna_sarscov2_24, n =Inf)[[1]]
res_ribo_sarscov2_24 <- topTags(qlf_ribo_sarscov2_24, n =Inf)[[1]]
#Another join
x$row_names <- rownames(x)
library("IHW")
res_RO_sarscov2_24$row_names <- rownames(res_RO_sarscov2_24)
res_RO_sarscov2_24 <- left_join(res_RO_sarscov2_24, x, by="row_names")
#IHW
ihwRes <- ihw(pvalues = res_RO_sarscov2_24$PValue, covariates = rowMeans(dplyr::select(res_RO_sarscov2_24, contains("virus"))), alpha = 0.1)
adj_pvalues(ihwRes)
res_RO_sarscov2_24$adj_palues <- adj_pvalues(ihwRes)
count(adj_pvalues(ihwRes) < 0.05)
#BH
FDR <- p.adjust(res_RO_sarscov2_24$PValue, method="BH")
sum(FDR < 0.1)
res_ribo_48$row_names <- rownames(res_ribo_48)
res_ribo_48 <- left_join(res_ribo_48, x, by="row_names")
#IHW
ihwRes <- ihw(pvalues = res_ribo_48$PValue, covariates = rowMeans(dplyr::select(res_ribo_48, contains("ribo_"))), alpha = 0.1)
adj_pvalues(ihwRes)
res_ribo_48$adj_palues <- adj_pvalues(ihwRes)
count(adj_pvalues(ihwRes) < 0.05)
#BH
FDR <- p.adjust(res_ribo_48$PValue, method="BH")
sum(FDR < 0.1)
res_rna_48$row_names <- rownames(res_rna_48)
res_rna_48 <- left_join(res_ribo_48, x, by="row_names")
#IHW
ihwRes <- ihw(pvalues = res_rna_48$PValue, covariates = rowMeans(dplyr::select(res_rna_48, contains("rna_"))), alpha = 0.1)
adj_pvalues(ihwRes)
res_rna_48$adj_palues <- adj_pvalues(ihwRes)
count(adj_pvalues(ihwRes) < 0.05)
#BH
FDR <- p.adjust(res_rna_48$PValue, method="BH")
sum(FDR < 0.1)
write.table(res_rna_48, str_c("./",'RNA_48.tsv'), sep="\t", quote = F, row.names = F)
write.table(res_ribo_48, str_c("./",'Ribo_48.tsv'), sep="\t", quote = F, row.names = F)
write.table(res_RO_48, str_c("./",'RO_48.tsv'), sep="\t", quote = F, row.names = F)
#3. 48
qlf_RO_sarscov2_48 <- glmQLFTest(fit, contrast = c(1, 0, -1, -1, 0, 1))
qlf_rna_sarscov2_48 <- glmQLFTest(fit, contrast = c(-1, 0, 1, 0, 0, 0))
qlf_ribo_sarscov2_48 <- glmQLFTest(fit, contrast = c(0, 0, 0, -1, 0, 1))
res_RO_sarscov2_48 <- topTags(qlf_RO_sarscov2_48, n = Inf)[[1]]
res_rna_sarscov2_48 <- topTags(qlf_rna_sarscov2_48, n =Inf)[[1]]
res_ribo_sarscov2_48 <- topTags(qlf_ribo_sarscov2_48, n =Inf)[[1]]
#Another join
x$row_names <- rownames(x)
library("IHW")
res_RO_sarscov2_48$row_names <- rownames(res_RO_sarscov2_48)
res_RO_sarscov2_48 <- left_join(res_RO_sarscov2_48, x, by="row_names")
#IHW
ihwRes <- ihw(pvalues = res_RO_sarscov2_48$PValue, covariates = rowMeans(dplyr::select(res_RO_sarscov2_48, contains("virus"))), alpha = 0.1)
adj_pvalues(ihwRes)
res_RO_sarscov2_48$adj_palues <- adj_pvalues(ihwRes)
count(adj_pvalues(ihwRes) < 0.05)
#BH
FDR <- p.adjust(res_RO_sarscov2_48$PValue, method="BH")
sum(FDR < 0.1)
res_ribo_sarscov2_48$row_names <- rownames(res_ribo_sarscov2_48)
res_ribo_sarscov2_48 <- left_join(res_ribo_sarscov2_48, x, by="row_names")
#IHW
ihwRes <- ihw(pvalues = res_ribo_sarscov2_48$PValue, covariates = rowMeans(dplyr::select(res_ribo_sarscov2_48, contains("ribo_"))), alpha = 0.1)
adj_pvalues(ihwRes)
res_ribo_sarscov2_48$adj_palues <- adj_pvalues(ihwRes)
count(adj_pvalues(ihwRes) < 0.05)
#BH
FDR <- p.adjust(res_ribo_sarscov2_48$PValue, method="BH")
sum(FDR < 0.1)
res_rna_sarscov2_48$row_names <- rownames(res_rna_sarscov2_48)
res_rna_sarscov2_48 <- left_join(res_ribo_sarscov2_48, x, by="row_names")
#IHW
ihwRes <- ihw(pvalues = res_rna_sarscov2_48$PValue, covariates = rowMeans(dplyr::select(res_rna_sarscov2_48, contains("rna_"))), alpha = 0.1)
adj_pvalues(ihwRes)
res_rna_sarscov2_48$adj_palues <- adj_pvalues(ihwRes)
count(adj_pvalues(ihwRes) < 0.05)
#BH
FDR <- p.adjust(res_rna_sarscov2_48$PValue, method="BH")
sum(FDR < 0.1)
write.table(res_rna_72, str_c("./",'RNA_72.tsv'), sep="\t", quote = F, row.names = F)
write.table(res_ribo_72, str_c("./",'Ribo_72.tsv'), sep="\t", quote = F, row.names = F)
write.table(res_RO_72, str_c("./",'RO_72.tsv'), sep="\t", quote = F, row.names = F)