-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathNEUTROPHILS.Rmd
59 lines (45 loc) · 2.3 KB
/
NEUTROPHILS.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
---
title: "*NEUTROPHIL*"
author: "Rosamaria Tricarico"
date: "February 1, 2018"
output: html_notebook
---
***
**Average fold change of gene expression from baseline and p-value for each gene displayed in a Volcano plot.** Marked in black are the genes with p-value > .01 and fold change > 1.3 to identify those which might play a role during fistula remodeling and/or failure.
**Patient#:** all patients
**Parameters of interest:** expression of genes extracted from neutrophils
**Paramters units:** [] fold changes from baseline
**Time points:** 2 weeks post-fistula
**Description:** The number of genes with larger expression variation is equal to 42; the list of genes is attached below.
**Goal:** to identify the list of genes potentially involved in fistula adaptmation and/or failure mechanisms.
***
```{r}
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})
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))
Significant <- subset(Table, pvalue<.01 & abs(log2FoldChange)>a)
write.csv(Significant, file = "Significant_Neutrophils.csv")
Significant
```