FISTULA.Rmd

---
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).**