-
Notifications
You must be signed in to change notification settings - Fork 0
/
lipidAcylBondsAnalysis.r
170 lines (162 loc) · 7.28 KB
/
lipidAcylBondsAnalysis.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
#this script investigates pd specifically
library(RMySQL)
setwd("D:/project/lipidomics/data/lipid analysis")
source("D:/project/lipidomics/data/lipidomicLib.r")
detectionLimits = read.table(file= "D:/project/lipidomics/data/detectionLimit.txt", header = F, sep = "\t")
lipidArray = detectionLimits[,1]
tData = getPatientSpecieMatrix(TRUE)
nData = getPatientSpecieMatrix(FALSE)
tData = getMolarData(tData,TRUE)
nData = getMolarData(nData,TRUE)
for(k in 1:length(lipidArray)){
#tPS = filterClassPSM(applyDetectionLimit(tData),lipidArray[k])
#nPS = filterClassPSM(applyDetectionLimit(nData),lipidArray[k])
tPS = filterClassPSM(tData,lipidArray[k])
nPS = filterClassPSM(nData,lipidArray[k])
carbonAcylList_N = list()
doubleBondList_N = list()
carbonAcylList_T = list()
doubleBondList_T = list()
listOfLipid = colnames(tPS)
listOfCarbonLength = vector()
listOfDoubleBond = vector()
if(length(listOfLipid) > 0){
for(i in 1:length(listOfLipid)){
acyl = unlist(strsplit(unlist(strsplit(listOfLipid[i], split ="-"))[1], split = ":"))
listOfCarbonLength[i] = acyl[1]
listOfDoubleBond[i] = acyl[2]
}
listOfCarbonLength = unique(listOfCarbonLength)
listOfDoubleBond = unique(listOfDoubleBond)
listOfCarbonLength = as.character(sort(as.numeric(listOfCarbonLength)))
listOfDoubleBond = as.character(sort(as.numeric(listOfDoubleBond)))
# Collect all lipid in normal
carbonLengthToRemove = vector()
doubleBondToRemove = vector()
tTestCarbonLengthToRemove = vector()
tTestDoubleBondToRemove = vector()
for(i in 1:length(listOfCarbonLength)){
carbonLengthToRemove = vector()
r = which(grepl(paste(listOfCarbonLength[i],":", sep=""),listOfLipid))
# sum all the lipid of the same axyl carbon length
if(length(r) > 1){
lc1 = rowSums(tPS[,r])
lc2 = rowSums(nPS[,r])
}
else {
lc1 = tPS[,r]
lc2 = nPS[,r]
}
ind = which(lc1 == 0)
if(length(ind) > 0)
for(i1 in 1:length(ind))
carbonLengthToRemove[length(carbonLengthToRemove)+1] = ind[i1]
carbonAcylList_T[[listOfCarbonLength[i]]] = lc1
ind = which(lc2 == 0)
if(length(ind) > 0)
for(i1 in 1:length(ind))
carbonLengthToRemove[length(carbonLengthToRemove)+1] = ind[i1]
carbonAcylList_N[[listOfCarbonLength[i]]] = lc2
if(length(carbonLengthToRemove) > 0){
carbonAcylList_T[[listOfCarbonLength[i]]] = carbonAcylList_T[[listOfCarbonLength[i]]][-carbonLengthToRemove]
carbonAcylList_N[[listOfCarbonLength[i]]] = carbonAcylList_N[[listOfCarbonLength[i]]][-carbonLengthToRemove]
}
if(length(carbonAcylList_T[[listOfCarbonLength[i]]]) < 2 | length(carbonAcylList_N[[listOfCarbonLength[i]]]) < 2){
tTestCarbonLengthToRemove[length(tTestCarbonLengthToRemove)+1] = i
}
}
if(length(tTestCarbonLengthToRemove) > 0)
listOfCarbonLength = listOfCarbonLength[-tTestCarbonLengthToRemove]
for(i in 1:length(listOfDoubleBond)){
doubleBondToRemove = vector()
r = which(grepl(paste(":",listOfDoubleBond[i], sep = ""),listOfLipid))
# sum all the lipid of the same double bond number
if(length(r) > 1){
lb1 = rowSums(tPS[,r])
lb2 = rowSums(nPS[,r])
}
else{
lb1 = tPS[,r]
lb2 = nPS[,r]
}
ind = which(lb1 == 0)
if(length(ind) > 0)
for(i1 in 1:length(ind))
doubleBondToRemove[length(doubleBondToRemove)+1] = ind[i1]
doubleBondList_T[[listOfDoubleBond[i]]] = lb1
ind = which(lb2 == 0)
if(length(ind) > 0)
for(i1 in 1:length(ind))
doubleBondToRemove[length(doubleBondToRemove)+1] = ind[i1]
doubleBondList_N[[listOfDoubleBond[i]]] = lb2
if(length(doubleBondToRemove) > 0){
doubleBondList_T[[listOfDoubleBond[i]]] = doubleBondList_T[[listOfDoubleBond[i]]][-doubleBondToRemove]
doubleBondList_N[[listOfDoubleBond[i]]] = doubleBondList_N[[listOfDoubleBond[i]]][-doubleBondToRemove]
}
if(length(doubleBondList_T[[listOfDoubleBond[i]]]) < 2 | length(doubleBondList_N[[listOfDoubleBond[i]]]) < 2){
tTestDoubleBondToRemove[length(tTestDoubleBondToRemove)+1] = i
}
}
if(length(tTestDoubleBondToRemove) > 0)
listOfDoubleBond = listOfDoubleBond[-tTestDoubleBondToRemove]
carbonAcylList = list()
carbonAcylTestPValue = vector()
for(i in 1:length(listOfCarbonLength)){
t = t.test(carbonAcylList_T[[listOfCarbonLength[i]]],carbonAcylList_N[[listOfCarbonLength[i]]], paired = T)
carbonAcylList[[i]] = log10(carbonAcylList_T[[listOfCarbonLength[i]]]/carbonAcylList_N[[listOfCarbonLength[i]]])
carbonAcylTestPValue[i] = t$p.value
}
tCarbonAcylTestColVec = vector(length = length(carbonAcylTestPValue))
for(i in 1:length(carbonAcylTestPValue)){
if(!is.nan(carbonAcylTestPValue[i]) && carbonAcylTestPValue[i] < 0.05){
tCarbonAcylTestColVec[i] = "red"
}else{
tCarbonAcylTestColVec[i] = "white"
}
}
doubleBondList = list()
doubleBondTestPValue = vector()
for(i in 1:length(listOfDoubleBond)){
t = t.test(doubleBondList_T[[listOfDoubleBond[i]]],doubleBondList_N[[listOfDoubleBond[i]]], paired = T)
doubleBondList[[i]] = log10(doubleBondList_T[[listOfDoubleBond[i]]]/doubleBondList_N[[listOfDoubleBond[i]]])
doubleBondTestPValue[i] = t$p.value
}
tDoubleBondTestColVec = vector(length = length(doubleBondTestPValue))
for(i in 1:length(doubleBondTestPValue)){
if(!is.nan(doubleBondTestPValue[i]) && doubleBondTestPValue[i] < 0.05){
tDoubleBondTestColVec[i] = "red"
}else{
tDoubleBondTestColVec[i] = "white"
}
}
#axis(1, labels = rownames(lipidP)[o], las = 2,cex = 0.5, at = 1:length(orVecL))
windows(width=10, height=8)
par(mfrow=c(1,2))
pdf(file = paste("Lipid_Acyl_Bonds_Analysis_", paste(lipidArray[k], ".jpg")), width = 9, height = 8, onefile = T)
boxplot(carbonAcylList,main = lipidArray[k], las = 2, names = listOfCarbonLength , ylab = "Lipid Ratio (T:N)", col = tCarbonAcylTestColVec)
abline(h=0, col = "blue")
color = unique(tCarbonAcylTestColVec)
if(length(color) == 1){
if(color == "white"){
legend("bottomright", legend = c("non signif."), fill =c("white"))
}else{
legend("bottomright", legend = c("signif."), fill =c("red"))
}
}else
legend("bottomright", legend = c("non signif.", "signif."), fill =c("white","red"))
#axis(1, labels = listOfCarbonLength, at = 1:length(listOfCarbonLength))
boxplot(doubleBondList,main = lipidArray[k], las = 2, names = listOfDoubleBond , ylab = "Lipid Ratio (T:N)", col = tDoubleBondTestColVec)
abline(h=0, col = "blue")
color = unique(tDoubleBondTestColVec)
if(length(color) == 1){
if(color == "white"){
legend("topleft", legend = c("non signif."), fill =c("white"))
}else{
legend("topleft", legend = c("signif."), fill =c("red"))
}
}else
legend("topleft", legend = c("non signif.", "signif."), fill = c("white","red"))
dev.off()
#axis(1, labels = listOfDoubleBond, at = 1:length(listOfDoubleBond))
}
}