-
Notifications
You must be signed in to change notification settings - Fork 0
/
FISTULA.Rmd
350 lines (306 loc) · 19.3 KB
/
FISTULA.Rmd
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
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
---
title: "*MULTISCALE ANALYSES ON ARTERIOVENOUS FISTULA MATURATION*"
author: "Rosamaria Tricarico"
date: "February 1, 2018"
output: html_notebook
---
***
**PROJECT OUTLINE**
Arteriovenous fistula for hemodialysis access fails to clinically mature in [28-53%](https://cjasn.asnjournals.org/content/clinjasn/1/2/332.full.pdf) of the cases.
My goal is to identify inter-correlations among time-changes in fistula morphology, hemodynamics, gene expression, and smooth muscle cell and endothelial cell dysfunctions.
These results can play a fundamental role in the understating of fistula maturation and failure causes, influencing clinical techniques as well as innovative strategies for fistula maturation.
***
**BACKGROUND**
Over 650,000 Americans are affected by [end stage renal disease](https://www.niddk.nih.gov/health-information/health-statistics/kidney-disease), which requires [hemodialysis treatment](https://www.niddk.nih.gov/health-information/kidney-disease/kidney-failure/hemodialysis) until availability for kidney transplant.
[Arteriovenous fistula](https://surgery.ucsf.edu/conditions--procedures/vascular-access-for-hemodialysis.aspx), a surgical connection between artery and vein of the arm, is recognized as the best hemodialysis access, as it has been reported to have lower thrombosis and infection rates compared to other fistulas.
However, obtaining clinically mature fistula is challenging and the causes of fistula maturation failures are unknown. [Research investigations](https://pubs.rsna.org/doi/pdf/10.1148/radiol.2015150385) are being performed to gain an insight into the mechanisms of maturation failure, from vascular anatomy and hemodynamics, to cell dysfunction and genomics that might play a role in vascular maturation, yet among the major limitations are small data-sets and lack of interconnection among events at micro- and macro-scales.
The [Division of Vascular Surgery at the University of Florida Shands Hospital](https://vascular.surgery.med.ufl.edu/research/national-institutes-of-health-study-of-av-fistulas/) has dedicated over 8 years on multiscale analyses on arteriovenous fistula maturation placing itself among the top six internationally recognized centers specialized on the topic. This vaste research has generated large patient-specific data-sets, now available for data science exploration.
***
**DATA-SETS**
The data-sets are composed by:
1. **Gene expression of monocytes and neutrophils (53,662 total) on 50 patients at two time points (baseline, 2 weeks).**
Scientific reasoning: circulating leukocytes (such as monocytes and neutrophils) play a role in development of intimal hyperplasia. Monocytes and neutrophils both exhibit phenotype plasticity that is associated with either pro-inflammatory or non-inflammatory states. Dynamic gene expression would give insight into their role during arteriovenous fistula adaption; unique responses might occur between fistulas that successfully mature and those that don’t.
RNA was isolated and expression levels were determined using microarrays.
2. **Anatomic and hemodynamic data along the fistula length of 50 patients at 3 time points (baseline, 6 weeks, 6 months).**
Scientific reasoning: fistula maturation consists of vessel adaptation to larger blood flow. For instance, the venous and arterial diameter increase and the vessel undergo structural change to withstand higher blood flow. If vascular adaptation mechanism fails, larger wall-shear stresses (WSS), wall-shear stress gradients (WSSG), and oscillatory shear index (OSI) can be observed.
Anatomic data was extracted from ultrasound imaging; hemodynamic data was extracted from patient-specific computational fluid-dynamic analyses.
3. **Brachial artery dilation on 50 patients at two time points (pre- and post-fistula).**
Scientific reasoning: Also endothelial and smooth muscle cells play a role in vessel wall adaption to stimuli. Flow-mediated and nitroglycerin-mediated dilation (FMD and NMD, respectively) are utilized for non-invasive assessment of endothelium-dependent and endothelium-independent (smooth muscle cell) vessel functions.
Dilation is measured with ultrasound imaging.
*The results related to this data-set are not be presented in this notebook.*
***
**EXPLORATORY ANALYSES**
We will first have an overview of the full anatomic-hemodyamic and genomic data-sets to undestand their major characteristics and distributions.
***1. GENOMICS***
*MONOCYTES*:
I extracted the average fold change of gene expression from baseline, for monocytes. Data sample of the fist five genes is displayed below.
```{r message=FALSE, warning=FALSE}
# Read the gene expression data set
Monocyte <- read.csv(file="Monocyte.csv",head=TRUE,sep=",")
p <- dim(Monocyte);
Gene <- Monocyte[1:1];
Data <- Monocyte[2:(p[2])];
n <- dim(Data);
m <- ((n[2])/2);
# Calculate the average fold change for each gene from baseline to 2-weeks after fistula
B_means <- rowMeans(Data[1:m]);
W2_means <- rowMeans(Data[(m+1):(n[2])]);
log2FoldChange <- (W2_means - B_means); #for semplicity gene expression is reported as log2(gene expression)
# Calculate the p-value of the t-test between baseline and 2-weeks after fistula
Pvalue <- apply(Data, 1, function(Data) {
t.test(x = Data[1:m], y = Data[(m+1):(n[2])])$p.value});
# Save data
Results = cbind(Gene, log2FoldChange);
Results = cbind(Results, pvalue = Pvalue);
write.csv(Results, file = "Monocyte_FC-PV.csv")
write.table(Results, file = "Monocyte_FC-PV.txt")
head(Results)
```
Then I created a Volcano plot to display the distribution of gene-specific fold change from baseline and marked the genes with p-value > .01 and fold change > 1.3 to identify those which might play a role during fistula remodeling and adaption.
```{r}
# Make a basic volcano plot
Table <- read.table("Monocyte_FC-PV.txt", header=TRUE)
with(Table, plot(log2FoldChange, -log10(pvalue), pch=1, main="Monocyte", xlim=c(-1,1)))
# Add colored points green if pvalue<.01 (alias -log10(pvalue)<2) and FoldChange>1.3
a <- log2(1.3);
with(subset(Table, pvalue>.01 | abs(log2FoldChange)<a), points(log2FoldChange, -log10(pvalue), pch=1, col="gray"))
with(subset(Table, pvalue<.01 & abs(log2FoldChange)>a), points(log2FoldChange, -log10(pvalue), pch=20, col="black"))
abline(v=c(+a,-a), col=c("black", "black"), lty=c(2,2), lwd=c(1, 1))
abline(h=2, col=c("black", "black"), lty=c(2,2), lwd=c(1, 1))
```
The genes with larger expression variation are listed below. The total number of significant genes from monocytes is equal to 7.
```{r echo=FALSE}
Significant <- subset(Table, pvalue<.01 & abs(log2FoldChange)>a)
Significant
write.csv(Significant, file = "Significant_Monocytes.csv")
```
*NEUTROPHILS*:
I replicated the same approach for neutrophil gene expression.
```{r echo=FALSE}
Neutrophil <- read.csv(file="Neutrophil.csv",head=TRUE,sep=",")
p <- dim(Neutrophil);
Gene <- Neutrophil[1:1];
Data <- Neutrophil[2:(p[2])];
n <- dim(Data);
m <- ((n[2])/2);
# Calculate the average fold change for each gene from baseline to 2-weeks after fistula
B_means <- rowMeans(Data[1:m]);
W2_means <- rowMeans(Data[(m+1):(n[2])]);
log2FoldChange <- (W2_means - B_means);
# Calculate the p-value of the t-test between baseline and 2-weeks after fistula
Pvalue <- apply(Data, 1, function(Data) {
t.test(x = Data[1:m], y = Data[(m+1):(n[2])])$p.value})
# Save data
Results = cbind(Gene, log2FoldChange);
Results = cbind(Results, pvalue = Pvalue);
write.csv(Results, file = "Neutrophil_FC-PV.csv")
write.table(Results, file = "Neutrophil_FC-PV.txt")
# Make a basic volcano plot
Table <- read.table("Neutrophil_FC-PV.txt", header=TRUE)
with(Table, plot(log2FoldChange, -log10(pvalue), pch=1, main="Neutrophil", xlim=c(-1,1)))
# Add colored points green if pvalue<.01 (alias -log10(pvalue)<2) and FoldChange>1.3
a <- log2(1.3);
with(subset(Table, pvalue>.01 | abs(log2FoldChange)<a), points(log2FoldChange, -log10(pvalue), pch=1, col="gray"));
with(subset(Table, pvalue<.01 & abs(log2FoldChange)>a), points(log2FoldChange, -log10(pvalue), pch=20, col="black"));
abline(v=c(+a,-a), col=c("black", "black"), lty=c(2,2), lwd=c(1, 1));
abline(h=2, col=c("black", "black"), lty=c(2,2), lwd=c(1, 1))
```
The genes with larger expression variation are listed below. The total number of significant genes from neutrophils is equal to 42.
```{r echo=FALSE}
Significant <- subset(Table, pvalue<.01 & abs(log2FoldChange)>a);
Significant
write.csv(Significant, file = "Significant_Neutrophils.csv")
```
Let's move to the anatomic-hemodynamic data-set.
***2. VESSELS' MORPHOLOGY AND HEMODYNAMICS***
*VEIN*:
the full anatomic and hemodynamic data-set of the *venous* portion of the fistula is displayed below.
```{r echo=FALSE}
# Read the hemodynamic data set
Hemodynamics <- read.csv(file="Hemodynamics_Fistula-Vein.csv", head=TRUE, sep=",")
## Load
library(lattice)
par(cex = .6)
options(scipen=999) # turn-off scientific notation like 1e+48
#--Define plot titles:
lab.Scan.1 <- "Baseline"
lab.Scan.2 <- "6 weeks"
lab.Scan.3 <- "6 months"
#--Custom strip function:
my.strip <- function(which.given, which.panel, ...) {
strip.labels <- c(lab.Scan.1, lab.Scan.2, lab.Scan.3)
panel.rect(0, 0, 1, 1, col="#ffe5cc", border=1)
panel.text(x=0.5, y=0.5, adj=c(0.5, 0.55), cex=0.95,
lab=strip.labels[which.panel[which.given]])}
#Create plots
xyplot(area ~ d.mm | Scan, data = Hemodynamics, groups = Patient, type = "l", lwd=0.5, grid = TRUE,
strip=my.strip, layout=c(1, 3, 1), xlab="Vein Length [mm]", ylab="Area [mm^2]")
xyplot(WSSave ~ d.mm | Scan, data = Hemodynamics, groups = Patient, type = "l", lwd=0.5, grid = TRUE,
strip=my.strip, layout=c(1, 3, 1), xlab="Vein Length [mm]", ylab="Average of WSS [Pa]")
xyplot(WSSmax ~ d.mm | Scan, data = Hemodynamics, groups = Patient, type = "l", lwd=0.5, grid = TRUE,
strip=my.strip, layout=c(1, 3, 1), xlab="Vein Length [mm]", ylab="Maximum WSS [Pa]")
xyplot(WSSGave ~ d.mm | Scan, data = Hemodynamics, groups = Patient, type = "l", lwd=0.5, grid = TRUE,
strip=my.strip, layout=c(1, 3, 1), xlab="Vein Length [mm]", ylab="Average of WSS gradient [1/Pa]")
xyplot(OSI ~ d.mm | Scan, data = Hemodynamics, groups = Patient, type = "l", lwd=0.5, grid = TRUE,
strip=my.strip, layout=c(1, 3, 1), xlab="Vein Length [mm]", ylab="Oscillatory Shear Index")
```
*ARTERY*:
an overview of the full anatomic and hemodynamic data-set follows for *arterial* portion of the fistula.
```{r echo=FALSE}
# Read the hemodynamic data set
Hemodynamics <- read.csv(file="Hemodynamics_Fistula-Artery.csv", head=TRUE, sep=",")
## Load
library(lattice)
par(cex = .6)
options(scipen=999) # turn-off scientific notation like 1e+48
#--Define plot titles:
lab.Scan.1 <- "Baseline"
lab.Scan.2 <- "6 weeks"
lab.Scan.3 <- "6 months"
#--Custom strip function:
my.strip <- function(which.given, which.panel, ...) {
strip.labels <- c(lab.Scan.1, lab.Scan.2, lab.Scan.3)
panel.rect(0, 0, 1, 1, col="#ffe5cc", border=1)
panel.text(x=0.5, y=0.5, adj=c(0.5, 0.55), cex=0.95,
lab=strip.labels[which.panel[which.given]])}
#Create plots
xyplot(area ~ d.mm | Scan, data = Hemodynamics, groups = Patient, type = "l", lwd=0.5, grid = TRUE,
strip=my.strip, layout=c(1, 3, 1), xlab="Artery Length [mm]", ylab="Area [mm^2]")
xyplot(WSSave ~ d.mm | Scan, data = Hemodynamics, groups = Patient, type = "l", lwd=0.5, grid = TRUE,
strip=my.strip, layout=c(1, 3, 1), xlab="Artery Length [mm]", ylab="Average of WSS [Pa]")
xyplot(WSSmax ~ d.mm | Scan, data = Hemodynamics, groups = Patient, type = "l", lwd=0.5, grid = TRUE,
strip=my.strip, layout=c(1, 3, 1), xlab="Artery Length [mm]", ylab="Maximum WSS [Pa]")
xyplot(WSSGave ~ d.mm | Scan, data = Hemodynamics, groups = Patient, type = "l", lwd=0.5, grid = TRUE,
strip=my.strip, layout=c(1, 3, 1), xlab="Artery Length [mm]", ylab="Average of WSS gradient [1/Pa]")
xyplot(OSI ~ d.mm | Scan, data = Hemodynamics, groups = Patient, type = "l", lwd=0.5, grid = TRUE,
strip=my.strip, layout=c(1, 3, 1), xlab="Artery Length [mm]", ylab="Oscillatory Shear Index")
```
As the *vein anatomy and hemodynamics* undergo to the largest changes, most visible from the area graphs, for the purposes of this brief presentation, from now on we will focus *vein* investigation.
*Moving to the next phase of the epxloratory analysis, we will dive in our data-sets in more depth. *
*By investigating the anatomy and hemodynamics on a patient-specific basis and compare it with the patient-specific genomic results, we can obtain a list of genes that are potentially involved with fistula remodelling and/or maturation failure mechanics. *
*Aim: identify whether unique responses occur for fistulas that successfully mature and those that don’t.*
***3. PATIENT 1 - MORPHOLOGY AND HEMODYNAMICS***
Fold changes from baseline are displayed below for the five parameters of interest at 6 weeks and 6 months post-intervention.
```{r echo=FALSE}
Patient1 <- filter(Hemodynamics, Patient == "031143 HJ")
Baseline <- filter(Patient1, Scan == 1)
Week6 <- filter(Patient1, Scan == 2)
Month6 <- filter(Patient1, Scan == 3)
n <- dim(Baseline);
m <- (n[1]);
FoldChange_Week6 <- Week6[5:9]/Baseline[5:9]
FoldChange_Month6 <- Month6[5:9]/Baseline[5:9]
xdim <- head(Patient1[4:4],m)
Patient1_FoldChange_Week6 = cbind(xdim, FoldChange_Week6);
Patient1_FoldChange_Month6 = cbind(xdim, FoldChange_Month6);
#--Define plot titles:
lab.area <- "Area"
lab.WSSave <- "Average of WSS"
lab.WSSmax <- "Maximum WSS"
lab.WSSGave <- "Average of WSS gradient"
lab.OSI <- "Oscillatory Shear Index"
#--Custom strip function:
my.strip <- function(which.given, which.panel, ...) {
strip.labels <- c(lab.area, lab.WSSave, lab.WSSmax, lab.WSSGave, lab.OSI)
panel.rect(0, 0, 1, 1, col="#ffe5cc", border=1)
panel.text(x=0.5, y=0.5, adj=c(0.5, 0.55), cex=0.95,
lab=strip.labels[which.panel[which.given]])}
#settings for graph title style
my.settings <- list(
par.main.text = list(font = 2, # make it bold
just = "left",
x = grid::unit(5, "mm")),
par.sub.text = list(font = 0.5,
just = "left",
x = grid::unit(5, "mm")))
#Creating the combined graph
Patient1_FoldChange = 0;
Patient1_FoldChange = rbind(Patient1_FoldChange_Week6, Patient1_FoldChange_Month6);
#View(Patient1_FoldChange)
Scan <- tail(Patient1[3:3],(n[1]*2));
Patient1_FoldChange = cbind(Scan,Patient1_FoldChange);
xyplot(area + WSSave + WSSmax + WSSGave + OSI ~ d.mm, data = Patient1_FoldChange, groups = Scan, scales=list(y="free", rot=0),type = "l",
grid = TRUE, strip=my.strip, outer=TRUE, layout=c(1, 5, 1), xlab="Vein Length [mm]", ylab="",
key=list(text=list(c("6 weeks", "6 months")),
title="Patient 31143 - Fold change from baseline",
col=c("deepskyblue3", "deeppink"), lty=c(1, 1),
columns=2, cex=0.7,
lines=TRUE))
```
This figure shows that along some portions of the vein the cross-sectional area increased up to 2.5 fold at six months, suggesting vessel adaptation and remodeling.
However, some other portions of the vein present close to no area increase. At these locations, hemodynamic parameters such as maximum (and average) WSS as well as average WSS gradient reach the highest values, underlying risk of platelets activation and thrombus formation.
For this reason, next we analyze the genomic of this patient to identify genes that might be involved in fistula maturation failure.
***4. PATIENT 1 - GENOMICS***
*MONOCYTES*:
the Volcano plot is displayed below and the genes with p-value > .01 and fold change > 1.3 are marked in blue to identify set of genes that might play a role during fistula failure.
```{r echo=FALSE}
Gene <- Monocyte[1:1];
Data <- Monocyte[2:(p[2])];
n <- dim(Data);
m <- ((n[2])/2);
Baseline <- Data$X31143.B;
Week2 <- Data$X31143.2W;
Patient1 <- data.frame(Gene,Baseline,Week2)
# Calculate the average fold change for each gene from baseline to 2-weeks after fistula
log2FoldChange <- (Week2 - Baseline);
# Calculate the p-value of the t-test between baseline and 2-weeks after fistula
Pvalue <- apply(Data, 1, function(Data) {
t.test(x = Data[1:m], y = Data[(m+1):(n[2])])$p.value})
# Save data
Results = cbind(Gene, log2FoldChange);
Results = cbind(Results, pvalue = Pvalue);
write.csv(Results, file = "Patient31143_Monocyte_FC-PV.csv")
write.table(Results, file = "Patient31143_Monocyte_FC-PV.txt")
# Make a basic volcano plot
Table <- read.table("Patient31143_Monocyte_FC-PV.txt", header=TRUE)
with(Table, plot(log2FoldChange, -log10(pvalue), pch=1, main="Patient 31143 Monocyte", xlim=c(-5,5)))
# Add colored points green if pvalue<.01 (alias -log10(pvalue)<2) and FoldChange>1.3
a <- log2(1.3);
with(subset(Table, pvalue>.01 | abs(log2FoldChange)<a), points(log2FoldChange, -log10(pvalue), pch=1, col="gray"));
with(subset(Table, pvalue<.01 & abs(log2FoldChange)>a), points(log2FoldChange, -log10(pvalue), pch=20, col="black"));
abline(v=c(+a,-a), col=c("black", "black"), lty=c(2,2), lwd=c(1, 1));
abline(h=2, col=c("black", "black"), lty=c(2,2), lwd=c(1, 1))
```
The list of significant change in gene expression from *monocytes* is here reported. Total number 124.
```{r echo=FALSE}
Significant <- subset(Table, pvalue<.01 & abs(log2FoldChange)>a)
Significant
write.csv(Significant, file = "Patient31143_Significant_Monocytes.csv")
```
*NEUTROPHILS*:
the Volcano plot is displayed below and the genes with p-value > .01 and fold change > 1.3 are marked in blue to identify set of genes that might play a role during fistula failure.
```{r echo=FALSE}
# Read the gene expression data set
Baseline <- Data$X31143.B; #will be changed with each patient code
Week2 <- Data$X31143.2W; #will be changed with each patient code
Patient1 <- data.frame(Gene,Baseline,Week2)
# Calculate the average fold change for each gene from baseline to 2-weeks after fistula
log2FoldChange <- (Week2 - Baseline);
# Calculate the p-value of the t-test between baseline and 2-weeks after fistula
Pvalue <- apply(Data, 1, function(Data) {
t.test(x = Data[1:m], y = Data[(m+1):(n[2])])$p.value})
# Save data
Results = cbind(Gene, log2FoldChange);
Results = cbind(Results, pvalue = Pvalue);
write.csv(Results, file = "Patient31143_Neutrophil_FC-PV.csv")
write.table(Results, file = "Patient31143_Neutrophil_FC-PV.txt")
# Make a basic volcano plot
Table <- read.table("Patient31143_Neutrophil_FC-PV.txt", header=TRUE)
with(Table, plot(log2FoldChange, -log10(pvalue), pch=1, main="Patient 31143 Neutrophil"))
# Add colored points green if pvalue<.01 (alias -log10(pvalue)<2) and FoldChange>1.3
a <- log2(1.3);
with(subset(Table, pvalue>.01 | abs(log2FoldChange)<a), points(log2FoldChange, -log10(pvalue), pch=1, col="gray"));
with(subset(Table, pvalue<.01 & abs(log2FoldChange)>a), points(log2FoldChange, -log10(pvalue), pch=20, col="black"));
abline(v=c(+a,-a), col=c("black", "black"), lty=c(2,2), lwd=c(1, 1));
abline(h=2, col=c("black", "black"), lty=c(2,2), lwd=c(1, 1));
```
The list of significant change in gene expression from *neutrophils* is here reported. Total number 1,193.
```{r echo=FALSE}
Significant <- subset(Table, pvalue<.01 & abs(log2FoldChange)>a)
Significant
write.csv(Significant, file = "Patient31143_Significant_Neutrophils.csv")
```
*By investigating the genomics of all patients which fistula already failied or show hemodynamic signatures of failure, we are able to obtain a list of genes potenitally involved with fistula failure mechanicsms.*
***
**Full data-sets can be found [here](https://github.com/RosamariaTricarico/2019TDI/tree/master/DATA.SETS).**
**Main R scripts can be found [here](https://github.com/RosamariaTricarico/2019TDI/tree/master/R.SCRIPTS).**