Loading in files

metadata <- read.delim(file="/Users/wede/Documents/University/Major_research_project/rnadyn/samples_metadata.txt", sep="|")
exons <- read.delim(file="/Users/wede/Documents/University/Major_research_project/rnadyn/rembrandts_results/exons.txt")
introns <- read.delim(file="/Users/wede/Documents/University/Major_research_project/rnadyn/rembrandts_results/introns.txt")
stability <-read.delim(file="/Users/wede/Documents/University/Major_research_project/rnadyn/rembrandts_results/stability.txt")

Convert to dataframes

df_exons <- as.matrix(exons[,-1])
rownames(df_exons) <- exons[,1]

df_introns <- as.matrix(introns[,-1])
rownames(df_introns) <- introns[,1]
df_stability <- as.matrix(stability[,-1])
rownames(df_stability) <- stability[,1]

Metadata data types

metadata[,c(16:17,19,21:22,33:38,40)] <- sapply(metadata[,c(16:17,19,21:22,33:38,40)], as.numeric)
## Warning in lapply(X = X, FUN = FUN, ...): NAs introduced by coercion

## Warning in lapply(X = X, FUN = FUN, ...): NAs introduced by coercion

## Warning in lapply(X = X, FUN = FUN, ...): NAs introduced by coercion

## Warning in lapply(X = X, FUN = FUN, ...): NAs introduced by coercion

## Warning in lapply(X = X, FUN = FUN, ...): NAs introduced by coercion

## Warning in lapply(X = X, FUN = FUN, ...): NAs introduced by coercion

Load packages for PCA plots

library(ggfortify)
## Loading required package: ggplot2
library(ggplot2)
library(grid)

Exon PCA

df_exon_pca <- prcomp(df_exons)
df_exon_pca <- as.data.frame(df_exon_pca$rotation)

cols = colnames(metadata)
exon_plots_pc12 = list()
exon_plots_pc34 = list()


for (i in 1:41) {
  p = ggplot(df_exon_pca, aes(x=PC1,y=PC2,label=rownames(df_exon_pca))) + geom_point()+ geom_text(aes(color=metadata[,i]), size=2,nudge_y=-0.01) + labs(x="PC1 (36.4%)", y="PC2 (25.1%)") + ggtitle(cols[i])
  exon_plots_pc12[[i]] = p
}

pdf("exon_plots_pc12.pdf")
for (i in 5:41) {
    print(exon_plots_pc12[[i]])
}

for (i in 1:41) {
  p = ggplot(df_exon_pca, aes(x=PC3,y=PC4,label=rownames(df_exon_pca))) + geom_point()+ geom_text(aes(color=metadata[,i]), size=2,nudge_y=-0.01) + labs(x="PC3 (9.6%)", y="PC4 (8.0%)") + ggtitle(cols[i])
  exon_plots_pc34[[i]] = p
}

pdf("exon_plots_pc34.pdf")
for (i in 5:41) {
    print(exon_plots_pc34[[i]])
}

dev.off()
## quartz_off_screen 
##                 2

Intron PCA

df_intron_pca <- prcomp(df_introns)
df_intron_pca <- as.data.frame(df_intron_pca$rotation)

intron_plots_pc12 = list()
intron_plots_pc34 = list()

for (i in 1:41) {
  p = ggplot(df_intron_pca, aes(x=PC1,y=PC2,label=rownames(df_intron_pca))) + geom_point() + geom_text(aes(color=metadata[,i]), size=2,nudge_y=-0.01) + labs(x="PC1 (34.3%)", y="PC2 (27.2%)") + ggtitle(cols[i])
  intron_plots_pc12[[i]] = p
}

pdf("intron_plots_pc12.pdf")
for (i in 5:41) {
    print(intron_plots_pc12[[i]])
}

for (i in 1:41) {
  p = ggplot(df_intron_pca, aes(x=PC3,y=PC4,label=rownames(df_intron_pca))) + geom_point() + geom_text(aes(color=metadata[,i]), size=2,nudge_y=-0.01) + labs(x="PC3 (9.8%)", y="PC2 (5.1%)") + ggtitle(cols[i])
  intron_plots_pc34[[i]] = p
}

pdf("intron_plots_pc34.pdf")
for (i in 5:41) {
    print(intron_plots_pc34[[i]])
}
dev.off()
## quartz_off_screen 
##                 2

Stability PCA

df_stability_pca <- prcomp(df_stability)
df_stability_pca <- as.data.frame(df_stability_pca$rotation)

stability_plots_pc12 = list()
for (i in 1:41) {
  p = ggplot(df_stability_pca, aes(x=PC1,y=PC2,label=rownames(df_stability_pca))) + geom_point() + geom_text(aes(color=metadata[,i]), size=2,nudge_y=-0.01) + labs(x="PC1 (38.0%)", y="PC2 (13.5%)", fill = cols[i]) + ggtitle(cols[i])
  stability_plots_pc12[[i]] = p
}

pdf("stability_plots_pc12.pdf")
for (i in 5:41) {
    print(stability_plots_pc12[[i]])
}

stability_plots_pc34 = list()
for (i in 1:41) {
  p = ggplot(df_stability_pca, aes(x=PC3,y=PC4,label=rownames(df_stability_pca))) + geom_point() + geom_text(aes(color=metadata[,i]), size=2,nudge_y=-0.01) + labs(x="PC3 (6.8%)", y="PC4 (4.3%)", fill = cols[i]) + ggtitle(cols[i])
  stability_plots_pc34[[i]] = p
}

pdf("stability_plots_pc34.pdf")
for (i in 5:41) {
    print(stability_plots_pc34[[i]])
}

stability_plots_pc56 = list()
for (i in 1:41) {
  p = ggplot(df_stability_pca, aes(x=PC5,y=PC6,label=rownames(df_stability_pca))) + geom_point() + geom_text(aes(color=metadata[,i]), size=2,nudge_y=-0.01) + labs(x="PC5 (3.8%)", y="PC6 (3.0%)", fill = cols[i]) + ggtitle(cols[i])
  stability_plots_pc56[[i]] = p
}

pdf("stability_plots_pc56.pdf")
for (i in 5:41) {
    print(stability_plots_pc56[[i]])
}
dev.off()
## quartz_off_screen 
##                 2

Selected plots exons

# PC1 - Platform/project
ggplot(df_exon_pca, aes(x=PC1,y=PC2,label=rownames(df_exon_pca))) + geom_point() + geom_text(aes(color=metadata[,5]), size=2,nudge_y=-0.01) + labs(x="PC1 (36.4%)", y="PC2 (25.1%)", fill = cols[5]) + ggtitle(cols[5])

ggplot(df_exon_pca, aes(x=PC1,y=PC2,label=rownames(df_exon_pca))) + geom_point() + geom_text(aes(color=metadata[,30]), size=2,nudge_y=-0.01) + labs(x="PC1 (36.4%)", y="PC2 (25.1%)", fill = cols[30]) + ggtitle(cols[30])

#PC3/PC4: Sex/platform
ggplot(df_exon_pca, aes(x=PC3,y=PC4,label=rownames(df_exon_pca))) + geom_point() + geom_text(aes(color=metadata[,5]), size=2,nudge_y=-0.01) + labs(x="PC3 (9.6%)", y="PC4 (8.0%)", fill = cols[5]) + ggtitle(cols[5])

ggplot(df_exon_pca, aes(x=PC3,y=PC4,label=rownames(df_exon_pca))) + geom_point() + geom_text(aes(color=metadata[,7]), size=2,nudge_y=-0.01) + labs(x= "PC3 (9.6%)", y="PC4 (8.0%)", fill = cols[7]) + ggtitle(cols[7])

ggplot(df_exon_pca, aes(x=PC3,y=PC4,label=rownames(df_exon_pca))) + geom_point() + geom_text(aes(color=metadata[,30]), size=2,nudge_y=-0.01) + labs(x= "PC3 (9.6%)", y="PC4 (8.0%)", fill = cols[7]) + ggtitle(cols[30])

# Selected plots introns

# PC1 - Platform/project/prep
ggplot(df_intron_pca, aes(x=PC1,y=PC2,label=rownames(df_intron_pca))) + geom_point() + geom_text(aes(color=metadata[,5]), size=2,nudge_y=-0.01) + labs(x="PC1 (34.3%)", y="PC2 (27.2%)") + ggtitle(cols[5])

ggplot(df_intron_pca, aes(x=PC1,y=PC2,label=rownames(df_intron_pca))) + geom_point() + geom_text(aes(color=metadata[,30]), size=2,nudge_y=-0.01) + labs(x="PC1 (34.3%)", y="PC2 (27.2%)") + ggtitle(cols[30])

ggplot(df_intron_pca, aes(x=PC1,y=PC2,label=rownames(df_intron_pca))) + geom_point() + geom_text(aes(color=metadata[,29]), size=2,nudge_y=-0.01) + labs(x="PC1 (34.3%)", y="PC2 (27.2%)") + ggtitle(cols[29])

#PC3/PC4: sex
ggplot(df_intron_pca, aes(x=PC3,y=PC4,label=rownames(df_intron_pca))) + geom_point() + geom_text(aes(color=metadata[,5]), size=2,nudge_y=-0.01) + labs(x="PC3 (9.8%)", y="PC4 (5.1%)", fill = cols[5]) + ggtitle(cols[5])

ggplot(df_intron_pca, aes(x=PC3,y=PC4,label=rownames(df_intron_pca))) + geom_point() + geom_text(aes(color=metadata[,7]), size=2,nudge_y=-0.01) + labs(x="PC3 (9.8%)", y="PC4 (5.1%)", fill = cols[7]) + ggtitle(cols[7])

# Selected plots stability

# PC2 - Platform/project/prep
ggplot(df_stability_pca, aes(x=PC1,y=PC2,label=rownames(df_stability_pca))) + geom_point() + geom_text(aes(color=metadata[,5]), size=2,nudge_y=-0.01) + labs(x="PC1 (38.0%)", y="PC2 (13.5%)") + ggtitle(cols[5])

#PC3: Platform/sex
ggplot(df_stability_pca, aes(x=PC3,y=PC4,label=rownames(df_stability_pca))) + geom_point() + geom_text(aes(color=metadata[,5]), size=2,nudge_y=-0.01) + labs(x="PC3 (6.8%)", y="PC4 (4.3%)", fill = cols[5]) + ggtitle(cols[5])

ggplot(df_stability_pca, aes(x=PC3,y=PC4,label=rownames(df_stability_pca))) + geom_point() + geom_text(aes(color=metadata[,30]), size=2,nudge_y=-0.01) + labs(x="PC3 (6.8%)", y="PC4 (4.3%)", fill = cols[30]) + ggtitle(cols[30])

ggplot(df_stability_pca, aes(x=PC3,y=PC4,label=rownames(df_stability_pca))) + geom_point() + geom_text(aes(color=metadata[,7]), size=2,nudge_y=-0.01) + labs(x="PC3 (6.8%)", y="PC4 (4.3%)", fill = cols[7]) + ggtitle(cols[7])

# Stability heatmap

library(RColorBrewer)
library(pheatmap)
coul <- colorRampPalette(rev(brewer.pal(8, "RdYlBu")))(26)
pheatmap(df_stability, scale="row",col=coul)

#library("heatmaply")
#heatmaply(df_stability)
exoncounts <- read.delim(file="/Users/wede/Documents/University/Major_research_project/rnadyn/exoncounts.txt")
introncounts <- read.delim(file="/Users/wede/Documents/University/Major_research_project/rnadyn/introncounts.txt")