01_IFPA_abstract_analysis.Rmd

---
title: "piRNA processing"
author: "Magda Price"
date: "Thursday, April 07, 2016"
output: html_document
---

```{r load libraries, message = F}
#source("http://www.bioconductor.org/biocLite.R")
#biocLite("DESeq")
library(DESeq)
library(ggplot2)
library(plyr)
library(dplyr)
library(reshape)
library(gplots)
library("RColorBrewer")
options("scipen"=100, "digits"=4)
library(matrixStats)
library(limma)
library(gridExtra)
```

```{r read data}
piRNA <- read.delim("Y:\\RobinsonLab\\Magda\\Project---microRNA\\01 - Input data\\01C_Placenta-piRNA-expression_edited.txt", row.names = 1)
head(piRNA)
colnames(piRNA)
piRNA <- as.matrix(piRNA)

desMat <- read.delim("Y:\\RobinsonLab\\Magda\\Project---microRNA\\01 - Input data\\01A_desMat_placenta_controls.txt", row.names = 1)
head(desMat)
rownames(desMat)

piRNA <- piRNA[,rownames(desMat)]
all(colnames(piRNA) == rownames(desMat))
dim(piRNA) # 6,029 piRNAs, 26 samples
```

```{r}
corMat <- desMat[,c(1,10,11,12)] # desMat input file changed since I did this. I think, gibven below, that originally I was interested in LibraryID (Tech1), trimester, sex and gestational age (since this was already a number?)
corMat$LibraryID <- as.numeric(corMat$LibraryID)
corMat$trimester <- as.numeric(corMat$trimester)
corMat$Sex <- as.numeric(corMat$Sex)

cor(corMat)
```

```{r Filtering of piRNAs}
# Q: How many piRNAs are not expressed in any sample? (N=2,979)
notExp <- piRNA[rowSums(piRNA) == 0,]
dim(notExp) 
nrow(notExp)/nrow(piRNA) 
# A: 49% of piRNAs not expressed

# Q: How many piRNAs are expressed in at least one sample? (N = 3,049)
piRNA.Exp <- piRNA[rowSums(piRNA) != 0,]
dim(piRNA.Exp) 
nrow(piRNA.Exp)/nrow(piRNA) 
# A: 51% of piRNAs expressed 

# Q: How many piRNAs expressed ONLY in one sample? (N = 2,134)
piRNA.count <- rowSums(piRNA.Exp > 0) 
head(piRNA.count)

Exp1 <- piRNA.Exp[names(piRNA.count[piRNA.count == 1]),]
dim(Exp1)
head(Exp1)
nrow(Exp1)/nrow(piRNA.Exp) 
# A: 70% of expressed piRNAs are only expressed in one sample

# Q: How highly expressed are the "one sample only" piRNAs?
hist(rowSums(Exp1), breaks = 100, main = "piRNAs expressed in ONLY one sample",
		 xlab = "level of expression (RPKM)")
# A: Most "one sample only" piRNAs have RPKM <20

# Q: How many piRNAs are expressed in MORE THAN ONE sample? (N = 915)
Exp <- piRNA.Exp[!rownames(piRNA.Exp) %in% rownames(Exp1),]
dim(Exp)
nrow(Exp)/nrow(piRNA.Exp) #30%

hist(Exp, breaks = 100, main = "piRNAs expressed in >1 sample",
		 xlab = "level of expression RPKM")

hist(log10(Exp+1), breaks = 100, main = "piRNAs expressed in >1 sample",
		 xlab = "level of expression log10(RPKM)")

# A: 30% of expressed piRNAs are expressed in more than one sample

Exp.NA <- Exp
Exp.NA[Exp.NA == 0] <- NA
```

```{r}
Exp26 <- piRNA.Exp[names(piRNA.count[piRNA.count == 26]),]
Exp26 <- Exp26[rowSums(Exp26) > 26,]
dim(Exp26)
head(Exp26)
nrow(Exp26)/nrow(piRNA.Exp) * 100

hist(rowSds(log10(Exp26 +1 )), breaks = 100)
Exp26.sd <- rowSds(log10(Exp26 +1))
```

```{r}
grey<-colorRampPalette(brewer.pal(n=5,"GnBu"))

plot.new()
heatmap.2(log10(Exp +1),    col = grey, dendrogram = "both",
				scale="none", margins=c(5,10),trace="none",
				labRow = "")
dev.off()

plot.new()
heatmap.2(cor(log10(Exp +1)),
          main="Correlation of all data",
          trace="none", col=grey,dendrogram="row",
          cexRow=1.0,
          cexCol=1.0,
          keysize=1)
dev.off()

source("C:/450k analyses/MTHFR project/heatmap.3.R")

col.trimester=gsub("3", "#EDCA0C", gsub("2", "#BC8F8F", gsub("1","#5f9EA0",desMat$trimester)))
col.sex=gsub("FEMALE", "#FF7352", desMat$Sex)
col.sex=gsub("MALE","#72CBEB", col.sex)
row.tech1=gsub("MX1303", "#D9D9D9",gsub("MX1307", "#F7F7F7", gsub("MX1310", "#BDBDBD", gsub("MX1355", "#737373",
																				gsub("MX1356", "#525252", gsub("MX1356", "#969696", gsub("MX1357", "#252525", desMat$LibraryID)))))))

row.tech2=gsub("C5JC1ACXX", "#F7F7F7", gsub("C5JC4ACXX", "#252525", gsub("C6RGTANXX", "#D9D9D9", desMat$Flowcell)))
#from this, I know that Tech1 = LibraryID, Tech2=Flowcell

clab=cbind(row.tech1, row.tech2, col.trimester, col.sex)
clab
colnames(clab)=c("Tech1", "Tech2", "Trimester", "Sex")


#Define custom dist and hclust functions for use with heatmaps
mydist=function(c) {dist(c,method="euclidian")}
myclust=function(c) {hclust(c,method="average")}

#Create heatmap using custom heatmap.3 source code loaded above
#pdf(file="C:\\450k analyses\\MTHFR project\\CHGC 2015\\heatmap.pdf", width=10, height=10)
main_title="" #Inter-sample correlation post processing (442,348 CpGs)
par(cex.main=1)
grey<-colorRampPalette(brewer.pal(n=9,"Greys"))

heatmap.3(cor(log10(Exp + 1)),hclustfun=myclust, distfun=mydist,na.rm = TRUE, scale="none", 
          dendrogram="both", margins=c(6,12), Rowv=TRUE, Colv=TRUE,  
          symbreaks=FALSE, key=TRUE, symkey=FALSE, density.info="none", 
          trace="none", ColSideColors=clab, 
          main=main_title, labCol=rownames(pData), 
          labRow=rownames(pData), cexRow=1, col=grey,ColSideColorsSize=5,
          KeyValueName="Pearson's r",keysize=1.2,
          cexCol=1)
legend("topright",legend=c("Female","Male",
                     " ","1st T","2nd T","Term" 
                     ),
       fill=c("#FF7352","#72CBEB",
              "white","#5f9EA0","#BC8F8F","#EDCA0C"
              ), 
       border=FALSE, bty="n", y.intersp = 0.9, cex=0.9)
```

piRNAs that are expressed in more than one sample (N=915)
```{r Descriptive stats}
# 1A. Features by trimester
with(desMat, table(trimester))

features.T <- 
	t(aggregate(t(log10(Exp + 1)),
								by = list(as.factor(desMat$trimester)), 
													 FUN = function(x) {
													 	c(count = sum(x >=1), 
      	mean = mean(x, na.rm = T), 
      	median = median(x, na.rm = T), 
      	max = max(x, na.rm = T),
      	sd = sd(x, na.rm = T))
													 	}))
head(features.T)
features.T <- as.data.frame(features.T[-1,])
str(features.T)
features.T$V1 <- as.numeric(features.T$V1)
features.T$V2 <- as.numeric(features.T$V2)
features.T$V3 <- as.numeric(features.T$V3)
str(features.T)

category <- do.call(rbind, strsplit(rownames(features.T), "\\."))
features.T <- cbind(features.T, category[,2])
head(features.T)
str(features.T)
features.T.melt <- melt(features.T)
colnames(features.T.melt) <- c("stat", "variable", "value")

ggplot(features.T.melt, aes(x = value, group = variable, fill = variable)) +
	geom_histogram(bins = 30) +
	theme_minimal() +
	facet_wrap( ~ stat, nrow = 2) +
	labs(x = "") +
	scale_fill_manual(values = c("dodgerblue3", "#5B9AD4", "#ED7D31"), labels = c("1st Trimester", "2nd Trimester", "Term"), name = "") +
	theme(legend.justification = c(1, 0), legend.position = c(1, 0)) +
	ggtitle("Descriptive stats for piRNAs")

# 1B. Features by trimester
features.sex <- 
	t(aggregate(t(log10(Exp + 1)),
								by = list(desMat$Sex), 
													 FUN = function(x) {
													 	c(count = sum(x >=1), 
      	mean = mean(x, na.rm = T), 
      	median = median(x, na.rm = T), 
      	max = max(x, na.rm = T),
      	sd = sd(x, na.rm = T))
													 	}))
head(features.sex)
features.sex <- as.data.frame(features.sex[-1,])
str(features.sex)
features.sex$V1 <- as.numeric(features.sex$V1)
features.sex$V2 <- as.numeric(features.sex$V2)
str(features.sex)

category <- do.call(rbind, strsplit(rownames(features.sex), "\\."))
features.sex <- cbind(features.sex, category[,2])
head(features.sex)	
str(features.sex)
features.sex.melt <- melt(features.sex, id.vars = "category[, 2]")
colnames(features.sex.melt) <- c("stat", "variable", "value")

ggplot(features.sex.melt, aes(x = as.numeric(value), group = variable, fill = variable)) +
	geom_histogram(bins = 30, position = "dodge") +
	theme_minimal() +
	facet_wrap( ~ stat, nrow = 2) +
	labs(x = "") +
	scale_fill_manual(values = c("red", "#5B9AD4"), labels = c("Female", "Male"), name = "") +
	theme(legend.justification = c(1, 0), legend.position = c(1, 0)) +
	ggtitle("Descriptive stats for piRNAs")
```


```{r}
sampleMelt.log <- melt(log10(Exp+1))
head(sampleMelt.log)
sampleMerge.log <- merge(sampleMelt.log, desMat, by.x = "X2", by.y = 0)

ggplot(sampleMerge.log, aes(x = reorder(X2, GA), y = value, color = Sex, fill = Sex))+
	geom_violin() +
	theme(axis.text.x = element_text(angle = 90, vjust = 0))

# Stats by column
col.statsLog <- 
apply(log10(Exp + 1), 2, function(x) {
      c(count = sum(x>0)/915, 
      	mean = mean(x, na.rm = T), 
      	median = median(x, na.rm = T), 
      	max = max(x, na.rm = T),
      	sd = sd(x, na.rm = T),
      	mad = mad(x, na.rm = T))
}
)

dat <-melt(t(col.statsLog))
datMelt <- merge(dat, desMat, by.x = "X1", by.y = 0)
head(datMelt)

ggplot(datMelt, aes(x = reorder(ID_Rob, trimester), y = value)) +
	theme_minimal() +
	geom_bar(stat = "identity", aes(fill = as.factor(trimester))) +
	facet_wrap( ~ X2, nrow = 2) +
	theme(axis.text.x = element_text(angle = 90, vjust = 0)) + 
	#scale_x_discrete(labels = reorder(desMat$ID_Rob, desMat$trimester), name = "") +
	scale_y_continuous("log10(RPKM)") +
	scale_fill_manual(values = c("dodgerblue3", "#5B9AD4", "#ED7D31"), labels = c("1st Trimester", "2nd Trimester", "Term"), name = "") +
	theme(legend.justification = c(0, 1), legend.position = c(0, 1))

ggplot(datMelt, aes(x = reorder(ID_Rob, as.numeric(Sex)), y = value)) +
	theme_minimal() +
	geom_bar(stat = "identity", aes(fill = as.factor(Sex))) +
	facet_wrap( ~ X2, nrow = 2) +
	theme(axis.text.x = element_text(angle = 90, vjust = 0)) + 
	#scale_x_discrete(labels = reorder(desMat$ID_Rob, desMat$trimester), name = "") +
	scale_y_continuous("log10(RPKM)") +
	scale_fill_manual(values = c("red", "dodgerblue3"), labels = c("Female", "Male"), name = "") +
	theme(legend.justification = c(0, 1), legend.position = c(0, 1))

ggplot(datMelt, aes(x = reorder(ID_Rob, as.numeric(LibraryID)), y = value)) +
	theme_minimal() +
	geom_bar(stat = "identity", aes(fill = as.factor(LibraryID))) +
	facet_wrap( ~ X2, nrow = 2) +
	theme(axis.text.x = element_text(angle = 90, vjust = 0)) + 
	#scale_x_discrete(labels = desMat$ID_Rob, name = "") +
	scale_y_continuous("log10(RPKM)") +
	#scale_fill_manual(name = "") +
	theme(legend.justification = c(0, 1), legend.position = c(0, 1)) +
	ggtitle("Stats by technical factor 1")

ggplot(datMelt, aes(x = reorder(ID_Rob, as.numeric(Flowcell)), y = value)) +
	theme_minimal() +
	geom_bar(stat = "identity", aes(fill = as.factor(Flowcell))) +
	facet_wrap( ~ X2, nrow = 2) +
	theme(axis.text.x = element_text(angle = 90, vjust = 0)) + 
#	scale_x_discrete(labels = desMat$ID_Rob, name = "") +
	scale_y_continuous("log10(RPKM)") +
	#scale_fill_manual(name = "") +
	theme(legend.justification = c(0, 1), legend.position = c(0, 1)) +
	ggtitle("Stats by technical factor 2")

ggplot(datMelt, aes(x = reorder(ID_Rob, ProcessingTime), y = value)) +
	theme_minimal() +
	geom_bar(stat = "identity", aes(fill = as.factor(trimester))) +
	facet_wrap( ~ X2, nrow = 2) +
	theme(axis.text.x = element_text(angle = 90, vjust = 0)) + 
	#scale_x_discrete(labels = datMelt$ProcessingTime, name = "") +
	scale_y_continuous("log10(RPKM)") +
	#scale_fill_manual(name = "") +
	theme(legend.justification = c(0, 1), legend.position = c(0, 1)) +
	ggtitle("Stats by technical factor 2")
```


```{r}
source("Y:\\RobinsonLab\\Magda\\Project---microRNA\\heat.scree.R")

desMat$trimester.2 <- as.factor(desMat$trimester)

meta_categorical<-desMat[,c(1,2,15,12)]  # input column numbers in meta that contain categorical variables
meta_continuous<-desMat[,c(10,11,13)] # input column numbers in meta that contain continuous variables
colnames(meta_categorical)<-c( "LibraryID","FlowCell","Sex", "Trimester (factor)")
colnames(meta_continuous)<-c("GA","Trimester (numeric)", "Processing Time")

princomp <- princomp(log(Exp+1))
Loadings.2<-as.data.frame(unclass(princomp$loadings))
vars.2 <- princomp$sdev^2
Importance.2<-vars.2/sum(vars.2)

heat_scree_plot(Loadings.2, Importance.2)
```

```{r, eval = F}
modMat <- cbind(desMat, princomp$loadings[rownames(desMat), 1:14])
modMat$trimester <- as.factor(modMat$trimester)

# Cell means paramaterization
(mod.2<- model.matrix(~ 0 + trimester + Sex + Comp.1 + Comp.2, modMat))
dmFit.2 <- lmFit(log10(Exp + 1), mod.2)
(contrast.matrix.2 <- makeContrasts(sex = SexMALE, 
																		T1_T2 = trimester1 - trimester2,
																		T1_T3 = trimester1 - trimester3,
																		T3_T2 = trimester3 - trimester2,
                               			levels=colnames(mod.2)))
fitContrasts.2 <- contrasts.fit(dmFit.2, contrast.matrix.2)
ebFit.2<-eBayes(fitContrasts.2)

		# piRNAs where sex is significant
		sex.2 <- topTable(ebFit.2, coef = "sex", adjust.method="BH", number=Inf)
		head(sex.2)
		hist(sex.2$P.Value, breaks = 100)

		# piRNAs where something changes betweeb trimesters 
		All.T <- topTable(ebFit.2, coef = grep("_T", colnames(coef(ebFit.2))), adjust.method="BH", number=Inf)
		head(All.T)
		hist(All.T$P.Value, breaks = 100)

#stable RNAs
topTable.Exp26 <- All.T[rownames(Exp26),]
hist(topTable.Exp26$P.Value)
sum(topTable.Exp26$P.Value > 0.1)
(61/6029)*100

```


```{r Hits mod1, eval = F}
dim(All.T[All.T$adj.P.Val <= 0.05,])
topGA <- All.T[All.T$adj.P.Val <= 0.05,]
		datGA <- Exp[rownames(All.T[All.T$adj.P.Val <= 0.05,]),]

		dim(datGA.expressed <- datGA[rowSums(datGA > 0) > 2,])

		hitsGA <- 
			t(aggregate(t(log10(datGA.expressed + 1)),
										by = list(desMat$trimester), 
															 FUN = function(x) {
															 	c( 
		      	
		      	mean = mean(x, na.rm = T), 
		      			      	sd = sd(x, na.rm = T))
															 	}))

head(hitsGA)
hitsGA <- as.data.frame(hitsGA[-1,])
str(hitsGA)
hitsGA$V1 <- as.numeric(hitsGA$V1)
hitsGA$V2 <- as.numeric(hitsGA$V2)
str(hitsGA)

fold <- hitsGA
fold$T1.T2 <- fold$V1 / fold$V2
fold$T1.T3 <- fold$V1 / fold$V3
fold$T3.T2 <- fold$V3 / fold$V2

summary(unlist(fold[,4:6]))
range(fold$T1.T3)
range(fold$T2.T3)


head(fold)

category <- do.call(rbind, strsplit(rownames(hitsGA), "\\."))
colnames(category) <- c("piRNA", "stat")
hitsGA <- cbind(hitsGA, category[,1:2])
head(hitsGA)	
str(hitsGA)
hitsGA.melt <- melt(hitsGA, id.vars = c("piRNA", "stat"))
head(hitsGA.melt)

ggplot(hitsGA.melt, aes(x = piRNA, y = value, group = variable, color = variable)) +
	 geom_point(position = position_jitter(width = .5), size = 5, alpha = 0.8) +
	theme_minimal() +
	facet_wrap( ~ stat, nrow = 2) +
	labs(x = "") +
	scale_color_manual(values = c("dodgerblue3", "lightgreen", "#ED7D31"), labels = c("1st Trimester", "2nd Trimester", "Term"), name = "") +
	#theme(legend.justification = c(0, 1), legend.position = c(0, 1)) +
	ggtitle("Descriptive stats for piRNAs") + 
	theme(axis.text.x = element_text(angle = 90, vjust = 0))

# The fold change parameter was calculated by dividing the mean expression value of tumours by the mean expression value in non-malignant tissues



hist(fold$T1.T2, breaks = 30, xlim = c(0,3))
hist(fold$T3.T2, breaks = 30, xlim = c(0,3))
hist(fold$T1.T3, breaks = 30, xlim = c(0,3))
```


```{r}
#All data - small model
(mod_sm<- model.matrix(~ as.factor(trimester) + Sex, desMat))
dmFit_sm <- lmFit(log10(Exp + 1), mod_sm)
colnames(coef(dmFit_sm))
ebFit_sm<-eBayes(dmFit_sm)

# 1. Residual standard error
mod_sm.rse <- dmFit_sm$sigma
hist(mod_sm.rse)

# 2. R-Squared
mod_sm.sst <- rowSums(log10(Exp + 1)^2)
mod_sm.ssr <- mod_sm.sst-dmFit_sm$df.residual * dmFit_sm$sigma^2
mod_sm.rs <- (mod_sm.ssr/mod_sm.sst)
hist(mod_sm.rs)


# piRNAs where sex is significant 
sex <- topTable(ebFit_sm, coef = "SexMALE", adjust.method="BH", number=Inf)
head(sex)
hist(sex$P.Value, breaks = 100)

# piRNAs where T2 is significant (STAT540 - Lecture 10, pg 55)
T2 <- topTable(ebFit_sm, coef = "as.factor(trimester)2", adjust.method="BH", number=Inf)
head(T2)
hist(T2$P.Value, breaks = 100)

# piRNAs where T3 significant (en masse) (STAT540 - Lecture 10, pg 55)
T3 <- topTable(ebFit_sm, coef = "as.factor(trimester)3", adjust.method="BH", number=Inf)
head(T3)
hist(T3$P.Value, breaks = 100)

dim(T3[T3$adj.P.Val <= 0.05,])


Exp[rownames(T3[T3$adj.P.Val <= 0.05,]),]

```



```{r Limma RAW, eval = F}
# Number of NA by trimester
counts.by.T <- t(aggregate.data.frame(t(is.na(Exp.NA)), by = list(desMat$trimester), FUN = sum))
counts.by.T <- counts.by.T[-1, ]
head(counts.by.T)

# % of NA by trimester
nor.counts  <- cbind(counts.by.T[,1] / 6, counts.by.T[,2] / 10, counts.by.T[,3] / 10)
head(nor.counts)

# 1/6 = .1666667, 2/10 = .2
piRNA_GA <- Exp[nor.counts[,1] <= .17 | nor.counts[,2] <= .20 | nor.counts[,3] <= .20,] 
dim(piRNA_GA) #187

(mod <- model.matrix(~ as.factor(trimester) + Sex + LibraryID, desMat))
dmFit <- lmFit(log10(piRNA_GA + 1), mod)
colnames(coef(dmFit))
ebFit<-eBayes(dmFit)

# piRNAs where T2 is significant (STAT540 - Lecture 10, pg 55)
T2 <- topTable(ebFit, coef = "as.factor(trimester)2", adjust.method="BH", number=Inf)
head(T2)
hist(T2$P.Value, breaks = 100)
abline(h=1.87,col="red")

# piRNAs where T3 significant (en masse) (STAT540 - Lecture 10, pg 55)
T3 <- topTable(ebFit, coef = "as.factor(trimester)3", adjust.method="BH", number=Inf)
head(T3)
hist(T3$P.Value, breaks = 100)
abline(h=1.87,col="red")
```

```{r Limma RAW - Sex, eval=FALSE}
with(desMat, table(Sex))

# Number of NA by Sex
counts.by.sex <- t(aggregate.data.frame(t(is.na(Exp.NA)), by = list(as.numeric(desMat$Sex)), FUN = sum))
counts.by.sex <- counts.by.sex[-1, ]
head(counts.by.sex)

# % of NA by trimester
nor.counts.sex  <- cbind(counts.by.sex[,1] / 16, counts.by.sex[,2] / 10)
head(nor.counts.sex)

# 12/16 = .75, 2/10 = .2
piRNA_sex <- Exp[nor.counts.sex[,1] <= .75 | nor.counts.sex[,2] <= .20,] 
dim(piRNA_sex) #307

(mod_sex <- model.matrix(~ Sex + as.factor(trimester) + LibraryID, desMat))
dmFit_sex <- lmFit(log10(piRNA_sex + 1), mod_sex)
colnames(coef(dmFit_sex))
ebFit_sex<-eBayes(dmFit_sex)

# piRNAs where sex is significant 
sex <- topTable(ebFit_sex, coef = "SexMALE", adjust.method="BH", number=Inf)
head(sex)
hist(sex$P.Value, breaks = 100)
abline(h=3.07,col="red")
```