rm(list=ls())setwd("/Users/nini/Desktop/2021lab/Hao/color_change_graph") #To update
library(ggplot2)library(Rmisc)library(gdata)library(RColorBrewer)library(paletteer)
data = read.csv("./genes_matrix_csv/expression_matrix.csv", h=F)data = data[,-1]col = read.csv("./genes_matrix_csv/columns_metadata.csv", h=T)row = read.csv("./genes_matrix_csv/rows_metadata.csv", h=T)
inv = c("1p22.1","1q31.3","2p22.3","2q22.1","3q26.1","6p21.33","6q23.1","7p14.3","7p11.2","7q11.22","7q36.1","8p23.1","11p12","11q13.2","12q13.11","12q21.2","14q23.3","17q21.31","21q21.3","Xq13.2")
age = c("Early prenatal (0-12 pcw)","Mid prenatal (13-24 pcw)","Late prenatal (25 pcw-birth)","Infancy (birth-1 yr)","Early childhood (1-6 yrs)","Late childhood (7-12 yrs)","Adolescence (13-18 yrs)","Youth (19-35 yrs)","Mid adulthood (> 35 yrs)")
result = matrix(NA,length(age)*length(inv),5)colnames(result) = c("Age","RPKM","Lower","Upper","Inv")
for (j in 1:length(inv)) { gene_list = unique(read.table(paste0("./Genes/",inv[j],".txt"),h=T)[,5]) gene = intersect(gene_list,row$gene_symbol) num = cbind(col, t(data[match(gene,row$gene_symbol),])) colnames(num)[9:length(num)] = gene num[num$age=="8 pcw"|num$age=="9 pcw"|num$age=="12 pcw",]$age ="Early prenatal (0-12 pcw)" num[num$age=="13 pcw"|num$age=="16 pcw"|num$age=="17 pcw"|num$age=="19 pcw"|num$age=="21 pcw"|num$age=="24 pcw",]$age ="Mid prenatal (13-24 pcw)" num[num$age=="25 pcw"|num$age=="26 pcw"|num$age=="35 pcw"|num$age=="37 pcw",]$age ="Late prenatal (25 pcw-birth)" num[num$age=="4 mos"|num$age=="10 mos",]$age ="Infancy (birth-1 yr)" num[num$age=="1 yrs"|num$age=="2 yrs"|num$age=="3 yrs"|num$age=="4 yrs",]$age ="Early childhood (1-6 yrs)" num[num$age=="8 yrs"|num$age=="11 yrs",]$age ="Late childhood (7-12 yrs)" num[num$age=="13 yrs"|num$age=="15 yrs"|num$age=="18 yrs",]$age ="Adolescence (13-18 yrs)" num[num$age=="19 yrs"|num$age=="21 yrs"|num$age=="23 yrs"|num$age=="30 yrs",]$age ="Youth (19-35 yrs)" num[num$age=="36 yrs"|num$age=="37 yrs"|num$age=="40 yrs",]$age ="Mid adulthood (> 35 yrs)" for (i in 1:length(age)) { result[(9*(j-1)+i),1] = age[i] result[(9*(j-1)+i),5] = inv[j] k=num[num$age==age[i],][,gene] k <- unmatrix(k,byrow=T) result[(9*(j-1)+i),2] = mean(k) result[(9*(j-1)+i),3] = CI(k, ci=0.95)[3] result[(9*(j-1)+i),4] = CI(k, ci=0.95)[1] } }
result = data.frame(result)level_order = age
####Combined####result[result$Inv=="Xq13.2",]$Inv = "23q13.2"level_order2 = invlevel_order2[20] = "23q13.2"
one <- "#170f11" #blacktwo <- "#7A3F00" #brownthree <- "#0471A6" #bluefour <- "#60C7FB" #lightbluefive <- "#C89933" #light brownsix <- "#335C3B" #dark green seven <- "#D7C9AA" #light browneight <- "#A5D5BF" #light greennine<- "#553EA3" #purple ten <- "#9A89D2" #light purpleeleven <- "#9CDCFC"twelve<- "#C0B5E3"thirteen<- "#FFA5AB"fourteen<- "#FDC05E"fifteen<- "#DB162F"sixteen<- "#93396A"seventeen<- "#8F8F8F"eighteen<- "#E89005"nineteen<- "#FADF63"twenty<- "#ABD888"
#color = c(brewer.pal(12,"Paired"),brewer.pal(8,"Accent"))#6 12 10 15 7col = c(rep(six,9),rep(twelve,9),rep(ten,9),rep(fifteen,9),rep(seven,9), #19 16 1 13 4 rep(nineteen,9),rep(sixteen,9),rep(one,9),rep(thirteen,9),rep(four,9), #8 5 11 20 3 rep(eight,9),rep(five,9),rep(eleven,9),rep(twenty,9),rep(three,9), #14 17 18 9 2 rep(fourteen,9),rep(seventeen,9),rep(eighteen,9),rep(nine,9),rep(two,9))
p <-ggplot(data=data.frame(result), aes(x=factor(Age, level=level_order), y=as.numeric(RPKM), group=Inv)) + geom_line(color= col) + labs(x="Age", y = "Mean brain expression (RPKM)")+ geom_ribbon(aes(ymin=as.numeric(result[,3]), ymax=as.numeric(result[,4]), fill=factor(result[,5],level = level_order2)), linetype=2, alpha=0.1) + #alpha 0.2 theme_classic() + theme(axis.text.x = element_text(angle = 60, hjust=1, size=16, face=2), axis.text.y = element_text(size=25, face=2), axis.title.x = element_text(size=20,face=2), axis.title.y = element_text(size=20,face=2), legend.margin = margin(0,0,0,-50), legend.title = element_text(size=16, face=2), legend.text = element_text(size=12, face=1)) + # #16 1 8 5 scale_fill_manual(values=c(sixteen, one, eight, five, #11 20 3 17 eleven, twenty,three, seventeen, #14 18 9 2 fourteen, eighteen, nine, two, #6 12 10 15 six, twelve,ten, fifteen, #7 19 13 4 seven, nineteen,thirteen, four),labels = inv, guide = guide_legend(title = "Inversion", override.aes = list(alpha = 1))) #alpha for the labels
ggsave("Age_vs_MeanBrainExpression.png",width = 11, height = 14)#####################scale_fill_discrete(name = "Inversion", labels = inv)# scale_fill_manual(values=col,name = "Inversion", labels = inv)+#scale_colour_manual(values=col)+# scale_fill_manual(values=col,name = "Inversion", labels = inv)+
#n$data