-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathPertScore calculation.R
39 lines (38 loc) · 1.06 KB
/
PertScore calculation.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
# PertScore calculate after Differentially expressed analysis
ttest = function(exp, a, b, c, d){
Pvalue = c()
FC = c()
logFC = c()
for(i in 1:nrow(exp)){
# un-impute
y = try(t.test(as.numeric(exp[i, a:b]), as.numeric(exp[i, c:d ])),silent=FALSE)
if('try-error' %in% class(y))
{
Pvalue = c(Pvalue, NA)
}else{
y = t.test(as.numeric(exp[i, a:b]), as.numeric(exp[i, c:d ]))
Pvalue = c(Pvalue, y$p.value)
}
pre = rowMeans( exp[i, a:b] , na.rm = TRUE)
post = rowMeans( exp[i, c:d ] , na.rm = TRUE)
FC =c(FC, pre/post)
logFC = c(logFC, log2(pre/post))
}
FDR=p.adjust(Pvalue, "BH")
out<-cbind(exp, Pvalue, FDR, FC, logFC)
out
}
ttest_out = ttest(dat, a, b,c,d)
sig = c()
for (i in 1:nrow(ttest_out)){
if (is.na(ttest_out[i, "Pvalue"]) ){
sig = c(sig, 0)
}else if (ttest_out[i, "Pvalue"]<0.05 & ttest_out[i, "logFC"]>log2(1.2)){
sig = c(sig, 1)
}else if (ttest_out[i, "Pvalue"]<0.05 & ttest_out[i, "logFC"]< (-log2(1.2))){
sig = c(sig, -1)
}else{
sig = c(sig, 0)
}
}
ttest_out$UpDown = sig