-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathedgeR-glm-analysis.R
158 lines (80 loc) · 3.16 KB
/
edgeR-glm-analysis.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
#1
genecounts <- read.table("18H18L.txt", row.names = 1)
colnames(genecounts) <- c(paste("18H",1:10,sep = "_"), paste("18L",1:10,sep="_"))
condition <-c(rep("H", 10), rep("L", 10))
dge <- DGEList(counts = genecounts, group = condition)
countsPerMillion <- cpm(dge)
countCheck <- countsPerMillion > 1
keep <- which(rowSums(countCheck)>=2)
dge <- dge[keep,]
dge <- calcNormFactors(dge)
design <- model.matrix(~ condition + 0, data = dge$samples)
colnames(design) = c("High", "Low")
disp <- estimateGLMCommonDisp(dge, design)
disp <- estimateGLMTrendedDisp(disp, design)
disp <- estimateGLMTagwiseDisp(disp, design)
fit <- glmFit(disp, design)
lrt <- glmLRT(fit)
res <- as.data.frame(lrt$table)
res <- cbind(res, FDR = p.adjust(res$PValue, method = "BH"))
fc1 <- res[res[,1] >= 1 | res[,1] <= -1,]
fc1p <- fc1[fc1$PValue <= 0.05,]
fc1q <- fc1[fc1$FDR <= 0.05,]
#2
genecounts <- read.table("18H18L.txt", row.names = 1)
colnames(genecounts) <- c(paste("18H",1:10,sep = "_"),paste("18L",1:10,sep="_"))
condition <-c(rep("H", 10), rep("L", 10))
dge <- DGEList(counts = genecounts, group = condition)
countsPerMillion <- cpm(dge)
countCheck <- countsPerMillion > 0.5
keep <- which(rowSums(countCheck)>=2)
dge <- dge[keep,]
dge <- calcNormFactors(dge)
design <- model.matrix(~ condition + 0, data = dge$samples)
colnames(design) = c("High", "Low")
disp <- estimateGLMCommonDisp(dge, design)
disp <- estimateGLMTrendedDisp(disp, design)
disp <- estimateGLMTagwiseDisp(disp, design)
fit <- glmFit(disp, design)
highlow <- makeContrasts(High-Low,levels = design)
lrt <- glmLRT(fit, contrast =highlow)
res <- as.data.frame(lrt$table)
res <- cbind(res, FDR = p.adjust(res$PValue, method = "BH"))
fc1 <- res[res[,1] >= 1 | res[,1] <= -1,]
fc1p <- fc1[fc1$PValue <= 0.05,]
fc1q <- fc1[fc1$FDR <= 0.05,]
write.table(fc1p,"fc1p_18H18L_1.txt", quote = F, sep = "\t", row.names = T, col.names = T)
#3
genecounts <- read.table("18H18L.txt", row.names = 1)
colnames(genecounts) <- c(paste("18H",1:10,sep = "_"),paste("18L",1:10,sep="_"))
condition <-c(rep("H", 10), rep("L", 10))
dge <- DGEList(counts = genecounts, group = condition)
countsPerMillion <- cpm(dge)
countCheck <- countsPerMillion > 0.5
keep <- which(rowSums(countCheck)>=2)
dge <- dge[keep,]
dge <- calcNormFactors(dge)
design <- model.matrix(~ condition + 0, data = dge$samples)
colnames(design) = c("High", "Low")
disp <- estimateGLMCommonDisp(dge, design)
disp <- estimateGLMTrendedDisp(disp, design)
disp <- estimateGLMTagwiseDisp(disp, design)
fit <- glmFit(disp, design)
lrt <- glmLRT(fit)
res <- as.data.frame(lrt$table)
res <- cbind(res, FDR = p.adjust(res$PValue, method = "BH"))
fc1 <- res[res[,1] >= 1 | res[,1] <= -1,]
fc1p <- fc1[fc1$PValue <= 0.05,]
fc1q <- fc1[fc1$FDR <= 0.05,]
#2-1
library(gplots)
data1 <- read.table("fc1p_18H18L_1.txt", header = T)
data2 <- read.table("fc1p_18H30H_1.txt", header = T)
data3 <- read.table("fc1p_18L30L_1.txt", header = T)
data4 <- read.table("fc1p_30H30L_1.txt", header = T)
test1 <- data1[,1]
test2 <- data2[,1]
test3 <- data3[,1]
test4 <- data4[,1]
venn(list(A = test1, B = test2, C = test3, D = test4))
A group = 105, B group = 1778, C group = 1665, D group = 47