-
Notifications
You must be signed in to change notification settings - Fork 0
/
02_DifferentialAbundance.R
86 lines (70 loc) · 3.52 KB
/
02_DifferentialAbundance.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
###############################################################################################
##' CyTOF data analysis of ALS patient PBMCs
##' 29-03-2021 Johannes Schroth
###############################################################################################
## Setup --------------------------------------------------------------------------------------------------------------------------------
#rm(list = ls())
set.seed(1234)
setwd('~/Dropbox/PhD/01 Data/Other/Ozlem/CyTOF data/Johannes CyTOF Analysis/Combined/')
library(dplyr)
library(edgeR)
## Differential abundance ---------------------------------------------------------------------------------------------------------------
counts <- df %>%
dplyr::count(metaclusters, sample_id) %>%
spread(key = 'sample_id', value = 'n') %>%
filter(metaclusters != 'Excluded') %>%
column_to_rownames('metaclusters') %>%
as.data.frame.matrix() %>%
replace(is.na(.), 0)
counts <- counts[!rowMeans(counts) < 5,]
group <- as.factor(md[match(colnames(counts), md$sample_id), 'Group'])
progression <- as.factor(md[match(colnames(counts), md$sample_id), 'Progression'])
onset <- as.factor(md[match(colnames(counts), md$sample_id), 'Onset'])
riluzole <- as.factor(md[match(colnames(counts), md$sample_id), 'Riluzole'])
date <- as.factor(md[match(colnames(counts), md$sample_id), 'CyTOF_run_date'])
nrow(counts)
DA_res_group <- counts %>%
DGEList(lib.size=colSums(counts)) %>%
estimateDisp(design = model.matrix(~ 0 + group + date)) %>%
glmQLFit(robust = T) %>%
glmQLFTest(fit, contrast = makeContrasts(CTRLvsALS = groupCTRL - groupALS,
levels = model.matrix(~ 0 + group + date))) %>%
.$table %>%
mutate(FDR = p.adjust(PValue, 'fdr')) %>%
rownames_to_column('Celltype') %>%
arrange(FDR)
DA_res_progression <- counts %>%
DGEList(lib.size=colSums(counts)) %>%
estimateDisp(design = model.matrix(~ 0 + progression + date)) %>%
glmQLFit(robust=F) %>%
glmQLFTest(fit, contrast = makeContrasts(FvsS = progressionF - progressionS,
levels = model.matrix(~ 0 + progression + date))) %>%
.$table %>%
mutate(FDR = p.adjust(PValue, 'fdr')) %>%
rownames_to_column('Celltype') %>%
arrange(FDR)
DA_res_onset <- counts %>%
DGEList(lib.size=colSums(counts)) %>%
estimateDisp(design = model.matrix(~ 0 + onset + date)) %>%
glmQLFit(robust=TRUE) %>%
glmQLFTest(fit, contrast = makeContrasts(BvsL = onsetB - onsetL,
levels = model.matrix(~ 0 + onset + date))) %>%
.$table %>%
mutate(FDR = p.adjust(PValue, 'fdr')) %>%
rownames_to_column('Celltype') %>%
arrange(FDR)
DA_res_riluzole <- counts %>%
DGEList(lib.size=colSums(counts)) %>%
estimateDisp(design = model.matrix(~ 0 + riluzole)) %>%
glmQLFit(robust=TRUE) %>%
glmQLFTest(fit, contrast = makeContrasts(YesvNo = riluzoleyes - riluzoleno,
levels = model.matrix(~ 0 + riluzole))) %>%
.$table %>%
mutate(FDR = p.adjust(PValue, 'fdr')) %>%
rownames_to_column('Celltype') %>%
arrange(FDR)
# Write results to folder
write.csv(DA_res_group, 'Datatables/Results/Differential Abundance/Group CTRLvALS results.csv', row.names = F)
write.csv(DA_res_progression, 'Datatables/Results/Differential Abundance/Progression SvF results.csv', row.names = F)
write.csv(DA_res_onset, 'Datatables/Results/Differential Abundance/Onset BvL results.csv', row.names = F)
write.csv(DA_res_riluzole, 'Datatables/Results/Differential Abundance/Riluzole YvN results.csv', row.names = F)