From 3d125d409ca8acb0bbf3771b3431954e22f51bb5 Mon Sep 17 00:00:00 2001 From: "Hao.Chen" <11779668+ewail@users.noreply.github.com> Date: Mon, 27 Apr 2020 21:11:23 +0800 Subject: [PATCH 1/5] update from liuwei --- CVDSBA_UMAP.Rmd | 180 +++++++++++++++++++ CVDSBA_heatmap.Rmd | 323 +++++++++++++++++++++++++++++++++ CVDSBA_mfuzz.Rmd | 363 +++++++++++++++++++++++++++++++++++++ common.R | 408 ++++++++++++++++++++++++++++++++++++++++++ common_change_label.R | 231 ++++++++++++++++++++++++ 5 files changed, 1505 insertions(+) create mode 100644 CVDSBA_UMAP.Rmd create mode 100644 CVDSBA_heatmap.Rmd create mode 100644 CVDSBA_mfuzz.Rmd create mode 100644 common.R create mode 100644 common_change_label.R diff --git a/CVDSBA_UMAP.Rmd b/CVDSBA_UMAP.Rmd new file mode 100644 index 0000000..f5860ef --- /dev/null +++ b/CVDSBA_UMAP.Rmd @@ -0,0 +1,180 @@ + +```{r metabo} + + +rm(list=ls()) +library(umap) +source("common.R") +source("datamining_library_ge20200306.R") +library(readr) +library(plyr) +library(readxl) +library(stringr) +library(magrittr) +library(Mfuzz) + +df<-read.csv("meta_matrix_log2delectNA08anddrug36_NAmin20200329(1).csv" ) + +row.names(df)<-c(as.matrix(df$X)) + +df2<-data.frame(t(df[-1])) + +df3<-data.frame(row.names(df2),df2) + +df3$row.names.df2.<-gsub("\\d+","",df3$row.names.df2.,perl = T) + +df1<-df3 + +names(df1)<-c("label",c(as.matrix(names(df1[,-c(1)])))) + +set.seed(2020) + +color<-c("firebrick1","forestgreen","dodgerblue","navyblue") + +pdf("UMAP_2020_allprotein_meta_847feature.pdf") + +drawUMAP(df1[,],ptColors = color,rowNormalization =F, colNormalization = T) + + +``` + + + + +```{r proteome} +rm(list=ls()) +library(umap) +source("common.R") +source("datamining_library_ge20200306.R") +library(readr) +library(plyr) +library(readxl) +library(stringr) +library(magrittr) +library(Mfuzz) + +info <- read_xlsx("sampleinfo2.xlsx") + +info$TMT <- gsub("^b","F",info$TMT) + + ################d +df2<-read.table("proteomic_matrix_ratio20200328.txt",header = T,sep = "\t",row.names = 1) + +df3<-data.frame(names(df2),t(df2)) + +names(df3)<-c("TMT",as.matrix(names(df3[,-1]))) + + +df4<-merge(info,df3,by.x = "TMT") + + +df1<-df4 +names(df1)<-c('patient',"label",c(as.matrix(names(df1[,-c(1:2)])))) + +df2<-data.frame(df1[-c(which(df1$label==c("NA"))),]) + +row.names(df2)<-c(as.matrix(df2$patient)) + +set.seed(2020) + +color<-c("firebrick1","forestgreen","dodgerblue","navyblue") + +pdf("UMAP_2020_allprotein.pdf") + +drawUMAP(df2[,-1],ptColors = color,rowNormalization =F, colNormalization = T) + +``` + +```{r proteome_severe_nonsevere_label} +rm(list=ls()) +library(umap) +source("common_change_label.R") +source("datamining_library_ge20200306.R") +library(readr) +library(plyr) +library(readxl) +library(stringr) +protein<-read.csv("proteomic_training_severe_nonsevere.csv",row.names = 1) +color=c("#DA942C","#8D2926") +set.seed(2020) +pdf("proteomic_training_severe_nonsevereABCD.pdf") +drawUMAP(protein,ptColors = color,rowNormalization =F, colNormalization = T) + +rm(list=ls()) +library(umap) +source("common_change_label.R") +source("datamining_library_ge20200306.R") +library(readr) +library(plyr) +library(readxl) +library(stringr) +protein<-read.csv("proteomic_training_severe_nonsevere.csv",row.names = 1) +protein<-protein[1:18,] +color=c("#DA942C") +set.seed(2020) +pdf("proteomic_training_nonsevere.pdf") +drawUMAP(protein,ptColors = color,rowNormalization =F, colNormalization = T) +rm(list=ls()) +# library(umap) +source("common_change_label.R") +source("datamining_library_ge20200306.R") +library(readr) +library(plyr) +library(readxl) +library(stringr) +protein<-read.csv("proteomic_training_severe_nonsevere.csv",row.names = 1) +protein<-protein[19:31,] +color=c("#8D2926","black") +set.seed(2020) +pdf("proteomic_training_severe.pdf") +drawUMAP(protein,ptColors = color,rowNormalization =F, colNormalization = T) + +``` + +```{r metabo_severe_nonsevere_label} +rm(list=ls()) +library(umap) +source("common_change_label.R") +source("datamining_library_ge20200306.R") +library(readr) +library(plyr) +library(readxl) +library(stringr) +protein<-read.csv("meta_training_severe_nonsevere.csv",row.names = 1) +color=c("#DA942C","#8D2926") +set.seed(2020) +pdf("meta_training_severe_nonsevere_ABCD.pdf") +drawUMAP(protein,ptColors = color,rowNormalization =F, colNormalization = T) +rm(list=ls()) +library(umap) +source("common_change_label.R") +source("datamining_library_ge20200306.R") +library(readr) +library(plyr) +library(readxl) +library(stringr) +protein<-read.csv("meta_training_severe_nonsevere.csv",row.names = 1) +protein<-protein[1:18,] +color=c("#DA942C") +set.seed(2020) +pdf("meta_training_nonsevere.pdf") +drawUMAP(protein,ptColors = color,rowNormalization =F, colNormalization = T) +rm(list=ls()) +# library(umap) +source("common_change_label.R") +source("datamining_library_ge20200306.R") +library(readr) +library(plyr) +library(readxl) +library(stringr) +protein<-read.csv("meta_training_severe_nonsevere.csv",row.names = 1) +protein<-protein[19:31,] +color=c("#8D2926","black") +set.seed(2020) +pdf("meta_training_severe.pdf") +drawUMAP(protein,ptColors = color,rowNormalization =F, colNormalization = T) + +``` + + + diff --git a/CVDSBA_heatmap.Rmd b/CVDSBA_heatmap.Rmd new file mode 100644 index 0000000..24e35d9 --- /dev/null +++ b/CVDSBA_heatmap.Rmd @@ -0,0 +1,323 @@ + + + +```{r heatmap_pathway_proteomics} + +rm(list = ls()) +library(readr) +library(plyr) +library(readxl) +library(stringr) +library(magrittr) +source("datamining_library_ge20200306.R") +library(RColorBrewer) +library(pheatmap) + + +matrix<-read.table("proteomic_matrix_delect5andNA15_ratio20200329.txt",header = T,sep = "\t",row.names = 1) + +protein<-read_excel("Fig3_heatmap_prots.xlsx",sheet=1) + + +Fig3prot<-c(as.matrix(protein$UniprotID)) + + +info<-read_excel("sTable1_20200329-1900.xlsx",sheet = 2) + +info2<-info[,c(3,4,14,15)] + +info2.1<-info[,c(3,5,14,15)] + + + +names(info2)<-c("type","patient","sex","age") + +names(info2.1)<-c("type","patient","sex","age") + + +matrix2<-data.frame(names(matrix),t(matrix)) + + + + +names(matrix2)<-c('patient',c(as.matrix(names(matrix2[,-c(1)])))) + + + +matrix3<-merge(info2,matrix2,by.x = "patient") + +matrix4<-merge(info2.1,matrix2,by.x = "patient") + + +matrix5<-rbind(matrix4,matrix3) + + +matrix5$type<-gsub("\\d+","",perl = T,matrix5$type) + + +matrix5$sex<-gsub("1","M",perl = T,matrix5$sex) + +matrix5$sex<-gsub("0","F",perl = T,matrix5$sex) + +# matrix5.1<-matrix5[-c(which(matrix5$type==c("jbdz"))),] + + + + +write.csv(matrix5,"heatmapfig2.csv",row.names = F) + + + +matrix5.1<-matrix5[] + + + +label<-matrix5.1[,c(1:4)] + +matrix6<-matrix5.1[,-c(1,3:4)] + +row.names(matrix6)<-c(as.matrix(label[,1])) + + +matrix7<-data.frame(matrix6[,c(1,which(names(matrix6)%in%Fig3prot))]) + + + +matrix7$type<-factor(matrix7$type,levels = c("jkdz","jbdz","PT","ZX")) + + + +matrix7<-matrix7[order(matrix7[,1]),] + + + + + + + +annotation_col<- data.frame(type = factor(label$type,levels = c("jkdz","jbdz","PT","ZX")), + + sex=label$sex, + age=label$age, + + + row.names = label$patient) + +type_color <- c("#85B22E","#5F80B4","#E29827","#922927") +names(type_color) <- c("jkdz","jbdz","PT","ZX") + +sex_color <- c("red","#016D06") + +names(sex_color) <- c("F","M") + +ann_colors <- list(type=type_color,sex=sex_color) +# +# colors = c( brewer.pal(11,"RdYlGn")[9:2]) + +colors = colorRampPalette(c("blue", "white","red" ))(1000) +matrix8<-data.frame(t(matrix7[,-1])) + + + +matrix9<-data.frame(row.names(matrix8),matrix8) + + +matrix9$row.names.matrix8.<-factor(matrix9$row.names.matrix8.,levels=c(as.matrix(Fig3prot))) + +matrix9<-matrix9[order(matrix9[,1]),] + + +gene<-protein[,c(5,6)] + + +for (i in 1:nrow(gene)) { + + + gene[i,3]<-paste0(gene[i,1],"_",gene[i,2]) + +} + + +row.names(matrix9)<-c(as.matrix(gene[,3])) + + + + +matrix10<-matrix9[,-1] + + + + + + + +matrix11<-data.frame(scale(matrix10,center = T)) + + + + + +pheatmap(matrix11,scale="row",color = colors, annotation_col = annotation_col, + annotation_colors = ann_colors, fontsize_col = 10, cluster_rows = T, cluster_cols = F,show_rownames =T, show_colnames = F,fontsize = 10,cellwidth=10,cellheight=10,filename = paste0("Heatmap for cluster ","Fig3_UN_fourtype",".pdf"),main = paste0("Heatmap for clusterFig3","")) + + + + +``` + +```{heatmapp_metabo} + + +rm(list = ls()) +library(readr) +library(plyr) +library(readxl) +library(stringr) +library(magrittr) +source("datamining_library_ge20200306.R") +library(RColorBrewer) +library(pheatmap) + + +all<-read.csv("metabo/meta_matrix_log2delectNA08anddrug36_NAmin20200329(1).csv") + +row.names(all)<-c(as.matrix(all[,1])) +protein<-read_excel("metabo/heatmap_0415.xlsx") +protein70<-c(as.matrix(protein[,1])) + +row.names(protein)<-c(as.matrix(protein[,1])) + +protein70order<- c(as.matrix(names(data.frame(t(protein))))) + +matrix<-data.frame(all[which(all$X%in%protein70),]) + +info<-read_excel("sTable1_20200329-1900.xlsx",sheet = 2) + + +info2<-info[,c(8,14,15)] + + + + +names(info2)<-c("type","sex","age") + + +matrix2<-data.frame(names(matrix[,-1]),t(matrix[,-1])) + + + +names(matrix2)<-c("type",c(as.matrix(matrix[,1]))) + + + + + + +matrix3<-merge(info2,matrix2,by.x = "type") + +matrix5<-data.frame(c(as.matrix(matrix3[,1])),matrix3) + + + + + +matrix5$type<-gsub("\\d+","",perl = T,matrix5$type) + + +matrix5$sex<-gsub("1","M",perl = T,matrix5$sex) + +matrix5$sex<-gsub("0","F",perl = T,matrix5$sex) + +# matrix5.1<-matrix5[-c(which(matrix5$type==c("jbdz"))),] +matrix5.1<-matrix5[,] + +label<-matrix5.1[,c(1:4)] + + + + +matrix6<-matrix5.1[,-c(1,3:4)] + +row.names(matrix6)<-c(as.matrix(label[,1])) + + +# matrix7<-data.frame(matrix6[,c(1,which(names(matrix6)%in%protein70))]) + + + +matrix6$type<-factor(matrix6$type,levels = c("jkdz","jbdz","PT","ZX")) + + + +matrix6<-matrix6[order(matrix6[,1]),] + + + + + + + +annotation_col<- data.frame(type = factor(label$type,levels = c("jkdz","jbdz","PT","ZX")), + + sex=label$sex, + age=label$age, + + + row.names = label$c.as.matrix.matrix3...1...) + +type_color <- c("#85B22E","#5F80B4","#E29827","#922927") +names(type_color) <- c("jkdz","jbdz","PT","ZX") + +sex_color <- c("red","#016D06") + +names(sex_color) <- c("F","M") + +ann_colors <- list(type=type_color,sex=sex_color) +# +# colors = c( brewer.pal(11,"RdYlGn")[9:2]) + +colors = colorRampPalette(c("blue", "white","red" ))(1000) + +matrix7<-data.frame(t(matrix6[,-1])) + + + +matrix8<-data.frame(row.names(matrix7),matrix7) + + + +matrix8$row.names.matrix7.<-factor(matrix8$row.names.matrix7.,levels=protein70order) + + + +matrix8<-matrix8[order(matrix8[,1]),] + + +row.names(matrix8)<-protein70 + +matrix9<-matrix8[,-1] + + +# matrix9<-data.frame(row.names(matrix8),matrix8) + + +# matrix9$row.names.matrix8.<-factor(matrix9$row.names.matrix8.,levels=c(as.matrix(Fig3prot))) + +# matrix9<-matrix9[order(matrix9[,1]),] +# +# +# +# matrix10<-matrix9[,-1] +# +write.csv(matrix9,"metabo/heatmap_meta.csv",row.names = T) +matrix9[is.na(matrix9)]<-0 + + +pheatmap(matrix9,scale = "row",color = colors, annotation_col = annotation_col, + annotation_colors = ann_colors, fontsize_col = 10, cluster_rows = F, cluster_cols = F,show_rownames =T, show_colnames = F,fontsize = 10,cellwidth=10,cellheight=10,filename = paste0("metabo/Heatmap ","Fig3_meta_uncluster_80_four_type_final20200415",".pdf"),main = paste0("Heatmap for clusterFig3_meta_PATHWAY_final","")) + +``` + + + + diff --git a/CVDSBA_mfuzz.Rmd b/CVDSBA_mfuzz.Rmd new file mode 100644 index 0000000..3795e62 --- /dev/null +++ b/CVDSBA_mfuzz.Rmd @@ -0,0 +1,363 @@ +--- +title: "covid19_mfuzz" +author: "liuwe" +date: "2020年3月26日" +output: html_document +--- +```{r mfuzz-proteomics} +rm(list = ls()) +library(readr) +library(plyr) +library(readxl) +library(stringr) +library(magrittr) +library(Mfuzz) +source("datamining_library_ge20200306.R") + +info <- read_xlsx("sampleinfo2.xlsx") +info$TMT <- gsub("^b","F",info$TMT) + + ################d + +df2<-read.table("proteomic_matrix_ratio20200328.txt",header = T,sep = "\t",row.names = 1) + +df3<-data.frame(names(df2),t(df2)) + +names(df3)<-c("TMT",as.matrix(names(df3[,-1]))) + + +df4<-merge(info,df3,by.x = "TMT") + + + +df5<-df4[-c(which(df4$Type=="NA")),-1] + + + + +df6<-aggregate(df5[,colnames(df5)[2:ncol(df5)]],by=list(df5$Type),mean,na.rm= TRUE) + +df6<-df6[c(1,2,3,4),] + +row.names(df6)<-c(as.matrix((df6[,1]))) + + +df7<-data.frame(t(df6[,-1])) + +df7<-data.frame(df7[,c(1,2,3,4)]) + +################## + +set.seed(2021) + +a<-ge.mfuzz.cselection(df7,range = seq(3,10,1)) + +b<-ge.mfuzz.getresult(a,8,"4_protein_mfuzz_ge2021type_8") + + + + +``` + +```{r anova-proteomics} +rm(list = ls()) +library(readr) +library(plyr) +library(readxl) +library(stringr) +library(magrittr) +source("datamining_library_ge20200306.R") + +info <- read_xlsx("sampleinfo2.xlsx") +info$TMT <- gsub("^b","F",info$TMT) + +df2<-read.table("proteomic_matrix_ratio20200328.txt",header = T,sep = "\t",row.names = 1) +df3<-data.frame(names(df2),t(df2)) + +names(df3)<-c("TMT",as.matrix(names(df3[,-1]))) + + +df4<-merge(info,df3,by.x = "TMT") + + + +df5<-df4[-c(which(df4$Type=="NA")),-1] + + +############################################ + +Matrix<-df5 + +for (K in 1:8) { + + + +mfuzz_prot<-read.csv(paste0("4_protein_mfuzz_ge2021type_8/mfuzz_",K,".csv")) + +prot<-as.matrix(mfuzz_prot$X) + +Matrix2<-data.frame(Matrix[,1],Matrix[,c(which(names(Matrix)%in% prot))]) + + + + +col<-ncol(Matrix2) +names(Matrix2)<-c("label",as.matrix(names(Matrix2[,-1]))) +anova<-data.frame(names(Matrix2[,-1])) + + +for (i in 2:col) { + aov<-(summary(aov(Matrix2[,i] ~ Matrix2[,1],Matrix2))[[1]])$`Pr(>F)`[1] + aov<-as.numeric(aov) + # aov<-as.numeric(aov,digits = 4, scientific = F) + anova[c(i-1),2]<-aov +} + +names(anova)<-c('protein',"pvalue") + +anova2<-data.frame(anova[c(which(anova$pvalue<=0.05)),]) + + +write.csv(anova2,paste0("4_protein_mfuzz_ge2021type_8/",K,"_mfuzz_anova_0.05.csv"),row.names = F) + + + +} + +``` + + + +```{r mfuzz-metabo} + +rm(list = ls()) +library(readr) +library(plyr) +library(readxl) +library(stringr) +library(magrittr) +source("datamining_library_ge20200306.R") +library(RColorBrewer) +library(pheatmap) +library(Mfuzz) + +matrix<-read.csv("meta/meta_matrix_log2delectNA08anddrug36_NAmin20200329(1).csv") + +row.names(matrix)<-c(as.matrix(matrix[,1])) + + + + +info<-read_excel("meta/sTable1_20200329-1900.xlsx",sheet = 2) + + +info2<-info[,c(8,14,15)] + + + + +names(info2)<-c("type","sex","age") + + +matrix2<-data.frame(names(matrix[,-1]),t(matrix[,-1])) + + + +names(matrix2)<-c("type",c(as.matrix(matrix[,1]))) + + +matrix3<-merge(info2,matrix2,by.x = "type") + +matrix5<-data.frame(c(as.matrix(matrix3[,1])),matrix3) + + + + + +matrix5$type<-gsub("\\d+","",perl = T,matrix5$type) + + +matrix5$sex<-gsub("1","M",perl = T,matrix5$sex) + +matrix5$sex<-gsub("0","F",perl = T,matrix5$sex) + + + + + +df5<-data.frame(matrix5[,-c(1,3,4)]) + + + +df6<-aggregate(df5[,colnames(df5)[2:ncol(df5)]],by=list(df5$type),mean,na.rm= TRUE) + + + +row.names(df6)<-c(as.matrix(df6$Group.1)) + + +df7.1<-data.frame(t(df6[,-1])) + +df7<-data.frame(df7.1[,c(1,3,4)]) + + +row.names(df7)<-row.names(matrix) + +################## + +set.seed(2021) + +a<-ge.mfuzz.cselection(df7,range = seq(3,10,1)) + +b<-ge.mfuzz.getresult(a,8,"meta/mfuzz_ge2021type_8_meta_new") + + + + +``` + +```{r anova-metabo} +rm(list = ls()) +library(readr) +library(plyr) +library(readxl) +library(stringr) +library(magrittr) +library(Mfuzz) +source("datamining_library_ge20200306.R") +rm(list = ls()) +library(readr) +library(plyr) +library(readxl) +library(stringr) +library(magrittr) +source("datamining_library_ge20200306.R") +library(RColorBrewer) +library(pheatmap) + + +matrix<-read.csv("meta/meta_matrix_log2delectNA08anddrug36_NAmin20200329(1).csv") + +row.names(matrix)<-c(as.matrix(matrix[,1])) + + + + +info<-read_excel("meta/sTable1_20200329-1900.xlsx",sheet = 2) + + +info2<-info[,c(8,14,15)] + + + + +names(info2)<-c("type","sex","age") + + +matrix2<-data.frame(names(matrix[,-1]),t(matrix[,-1])) + + + +names(matrix2)<-c("type",c(as.matrix(matrix[,1]))) + + + + + + +matrix3<-merge(info2,matrix2,by.x = "type") + +matrix5<-data.frame(c(as.matrix(matrix3[,1])),matrix3) + + + + + +matrix5$type<-gsub("\\d+","",perl = T,matrix5$type) + + +matrix5$sex<-gsub("1","M",perl = T,matrix5$sex) + +matrix5$sex<-gsub("0","F",perl = T,matrix5$sex) + + + + + +df5<-data.frame(matrix5[,-c(1,3,4)]) + + + + + + + +############################################ + +df6<-data.frame(df5[-c(which(df5$type=="jbdz")),]) + +Matrix<-df6 + + + + +for (K in 1:8) { + +mfuzz_prot<-read.csv(paste0("meta/mfuzz_ge2021type_8_meta_new/mfuzz_",K,".csv")) + + +row.names(mfuzz_prot)<-c(as.matrix(mfuzz_prot[,1])) + + +mfuzz_prot2<-data.frame(t(mfuzz_prot)) + +prot<-as.matrix(names(mfuzz_prot2)) + + + + + + + +Matrix2<-data.frame(Matrix[,1],Matrix[,c(which(names(Matrix)%in% prot))]) + + +col<-ncol(Matrix2) +names(Matrix2)<-c("label",as.matrix(names(Matrix2[,-1]))) +anova<-data.frame(names(Matrix2[,-1])) + + +for (i in 2:col) { + aov<-(summary(aov(Matrix2[,i] ~ Matrix2[,1],Matrix2))[[1]])$`Pr(>F)`[1] + aov<-as.numeric(aov) + # aov<-as.numeric(aov,digits = 4, scientific = F) + anova[c(i-1),2]<-aov +} + +names(anova)<-c('protein',"pvalue") + +anova2<-data.frame(anova[c(which(anova$pvalue<=0.05)),]) + + +write.csv(anova2,paste0("meta/mfuzz_ge2021type_8_meta_new/",K,"_mfuzz_anova_0.05.csv"),row.names = F) + + +} + + + +``` + + + + + + + + + + + + + + diff --git a/common.R b/common.R new file mode 100644 index 0000000..e1c858e --- /dev/null +++ b/common.R @@ -0,0 +1,408 @@ +library("FactoMineR") +library("factoextra") +library("grid") +library("data.table") +library("ggplot2") +library(pheatmap); +library("RColorBrewer"); +library("tsne"); + + +###tSNE plot DE: dataframe with rows as sample and columns as features,a column with name 'label' is required, represents the label of samples +# lblColors: a vector containing colors for each label, like lblColors=c(M = "forestgreen", N = "gray0", P="firebrick", C= "red",A="blue") +drawTSNE <- function(DF,ptColors,rowNormalization=F,colNormalization=F,perplexity=10,strTitle='tSNE'){ + M <- DF[,colnames(DF)!='label'] + if(rowNormalization){M <- data.frame(t(apply(M,1,function(v){(v-mean(v,na.rm=T))/sd(v,na.rm=T)})))} + if(colNormalization){M <- apply(M,2,function(v){(v-mean(v,na.rm=T))/sd(v,na.rm=T)})} + M[is.na(M)] <- 0 + indx <- match('label',colnames(M)) + clnames <- colnames(DF)[colnames(DF)!='label'] + tsn = tsne(M,perplexity =perplexity) + cnames <- colnames(M) + tsn <- data.frame(tsn,DF$label) + colnames(tsn) <- c("X","Y","label") + rownames(tsn) <- rownames(M) + #tsn <- tsn[-c(which(tsn$X==min(tsn$X)),which(tsn$Y==min(tsn$Y))),] + #tsn <- tsn[-c(which(tsn$X==max(tsn$X)),which(tsn$Y==max(tsn$Y))),] + #lblColors <- c(N='#537e35',M='#e17832',A='#f5b819',C='#5992c6',P='#282f89',W='mediumorchid3') + + #lblColors <- c(A='#537e35',M='#e17832',N='#f5b819',B='#5992c6',C='#282f89',W='mediumorchid3') + p <- ggplot(tsn, aes(x=X, y=Y, colour=label)) + geom_point(size=4) + p <- p + theme( panel.grid.major = element_blank(), + panel.grid.minor = element_blank(), + panel.border = element_blank(), + plot.title = element_text(size=15), + axis.line.x = element_line(color="black", size = 0.5), + axis.line.y = element_line(color="black", size = 0.5), + panel.background = element_blank()) + p <- p + labs(title=strTitle) + p <- p + scale_colour_manual(values=ptColors) + p +} + +drawPCA<- function(DF,ptColors,rowNormalization=F,colNormalization=F,strTitle=NULL){ +## M is a matrix or dataframe, rows are samples and columns are features, rownames are sample names + M <- DF[,colnames(DF)!='label'] + if(rowNormalization){ + M <- data.frame(t(apply(M,1,function(v){(v-mean(v,na.rm=T))/sd(v,na.rm=T)}))) + #print('Normalization by row is done!') + } + if(colNormalization){ + M <- apply(M,2,function(v){(v-mean(v,na.rm=T))/sd(v,na.rm=T)}) + } + clnames <- colnames(DF)[colnames(DF)!='label'] + M[is.na(M)] <- 0 + m1 <- prcomp(M,colNormalization); + Y <- scale(M, m1$center, m1$scale) %*% m1$rotation + Y <- Y[,c(1,2)] + + Y <- data.frame(Y,DF$label); + colnames(Y) <- c("PC1","PC2","label") + if(is.null(strTitle)){ + strTitle <- sprintf("PCA:%d features",length(clnames)) + } + eigs <- m1$sdev^2 + percentages <- eigs[1:2] / sum(eigs) + # lblColors <- c(N='#537e35',M='#e17832',A='#f5b819',C='#5992c6',P='#282f89',W='mediumorchid3') + #lblColors <- c(training='#537e35',validation='#e17832',A='#f5b819',C='#5992c6',P='#282f89',W='mediumorchid3') + p <- ggplot(Y, aes(x=PC1, y=PC2, colour=label)) + geom_point(size=4) + p <- p + theme( panel.grid.major = element_blank(), + panel.grid.minor = element_blank(), + panel.border = element_blank(), + axis.line.x = element_line(color="black", size = 0.25), + axis.line.y = element_line(color="black", size = 0.25), + plot.title = element_text(size=16), + panel.background = element_blank()) + + strLabx <- sprintf("PC1(%4.2f%%)",percentages[1]*100) + p <- p + labs(x =strLabx,y = sprintf("PC2(%4.2f%%)",percentages[2]*100), + title =strTitle) + p <- p + scale_colour_manual(values=ptColors) + + p +} + +# parse gene symbol from uniprot website +# acc: protein access Id like ' "O00264" "O00468" "O14773" "O14979" "O15230" "O43175" "O43707" +# return: the gene symbol for the acc +geneName <- function(acc){ + kk <- read.table(sprintf('https://www.uniprot.org/uniprot/%s.fasta',acc),sep="\n",stringsAsFactors=F,header=F) + a <- strsplit(kk[1,]," ")[[1]] + b <- a[grepl("GN=",a)] + strsplit(b,"=")[[1]][2] +} +getSubset <- function(Labels,M0,missRate,imputation=F){ + lbls <- Labels + lbls <- if(length(Labels)==1){lbls <- unlist(strsplit(lbls,""))}else{lbls <- Labels} + tmp <- M0[M0$label %in% lbls,] + R0 <- apply(tmp,2,function(v){sum(is.na(v))/length(v)*100}) + tmp <- tmp[,R0<=missRate] + + #if(!imputation){tmp[is.na(tmp)] <- 0} + if(imputation){ #### kNN imputation + K=3 + t1 <- sapply(1:dim(tmp)[1],function(i) { + lbl <- tmp[i, 'label'] + v0 <- tmp[i, colnames(tmp) != 'label'] + feature0 <- names(v0)[is.na(v0)] + if(length(feature0)>0){ + k0 <- tmp[-i,] + k0 <- k0[k0$label==lbl,-1] + k0 <- rbind(v0,k0) + k0 <- apply(k0,2,function(v){(v-mean(v,na.rm=T))/sd(v,na.rm=T)}) + k0[is.na(k0)] <- 0 + + D0 <- as.matrix(dist(k0)) + d0 <- order(D0[1,])[2:(K+1)] + if(length(feature0)==1){ + v0[feature0] <- mean(tmp[d0,feature0],na.rm=T) + }else{ + v0[feature0] <- apply(tmp[d0,feature0],2,function(v){mean(v,na.rm=T)}) + } + } + v0 + }) + + t1 <- apply(t(t1),2,unlist) + rownames(t1) <- rownames(tmp) + t1 <- data.frame(t1) + clnames <- colnames(t1) + t1$label <- tmp$label + t1 <- t1[,c('label',clnames)] + #print(sprintf("get subset for %s with %d rows and %d features",paste0(lbls,collapse=","),dim(t1)[1],dim(t1)[2])) + tmp <- t1 + } + return(tmp) +} + + + RFScore <- function(feature,M1=t1,nFolds=10,it=10){ + #tmp <- data.frame(t(apply(M1[,feature],1,function(v0){v <- as.numeric(v0); (v-mean(v))/sd(v)}))) + tmp <- M1[,c('label',feature)] + clsses <- unique(tmp$label) + #rownames(tmp) <- sapply(1:dim(tmp)[1],function(v){paste0('R',v)}) + sampleNumber <- sapply(clsses,function(v){sum(tmp$label==v)}) + names(sampleNumber) <- clsses + maxLabel <- names(sampleNumber)[max(sampleNumber)==sampleNumber] + minLabel <- names(sampleNumber)[min(sampleNumber)==sampleNumber] + size <- min(sampleNumber) + + avgAcc <- 0; + probs <- rep(0,max(sampleNumber)) + names(probs) <- rownames(tmp)[tmp$label==maxLabel] + its <- 0 + indx2 <- rownames(tmp)[tmp$label==minLabel] + + prediction <- matrix(0,nrow=dim(M1)[1],ncol=2) + rownames(prediction) <- rownames(M1) + + while(sum(probs==0)>0){ + indx1 <- names(sort(probs))[1:size] + probs[indx1] <- probs[indx1]+1 + tmp1 <- tmp[c(indx1,indx2),] + for(fd in 1:nFolds){ + folds <- createFolds(tmp1$label,nFolds) + for(fold in folds){ + valids <- tmp1[fold,] + rownames(valids) <- rownames(tmp1)[fold] + + trains <- tmp1[setdiff(1:dim(tmp1)[1],fold),] + trains$label <- as.factor(trains$label) + + tmpRF <- randomForest(label ~ . ,data=trains,importance=T,ntree=1000,nodesize=5) + predicted <- predict(tmpRF,valids,type='prob') + prediction[rownames(predicted),] <- predicted+prediction[rownames(predicted),] + + } + } + }#while + colnames(prediction) <- colnames(predicted) + prediction <- data.frame(t(apply(prediction,1,function(v){v/sum(v)}))) + prediction$predicted <- as.character(apply(prediction,1,function(v){names(v)[v==max(v)]})) + prediction$observed <- M1$label + rtv <-sprintf("%s %4.3f",paste0(feature,collapse=","),sum(prediction$observed==prediction$predicted)/dim(M1)[1]*100) + #prediction + rtv + } +############################################################################################# +RFPredict <- function(trainM,testM){ + + clsses <- unique(trainM$label) + #rownames(tmp) <- sapply(1:dim(tmp)[1],function(v){paste0('R',v)}) + sampleNumber <- sapply(clsses,function(v){sum(tmp$label==v)}) + names(sampleNumber) <- clsses + maxLabel <- names(sampleNumber)[max(sampleNumber)==sampleNumber] + minLabel <- names(sampleNumber)[min(sampleNumber)==sampleNumber] + size <- min(sampleNumber) + + avgAcc <- 0; + probs <- rep(0,max(sampleNumber)) + names(probs) <- rownames(tmp)[tmp$label==maxLabel] + its <- 0 + indx2 <- rownames(tmp)[tmp$label==minLabel] + + prediction <- matrix(0,nrow=dim(M1)[1],ncol=2) + rownames(prediction) <- rownames(M1) + + while(sum(probs==0)>0){ + indx1 <- names(sort(probs))[1:size] + probs[indx1] <- probs[indx1]+1 + tmp1 <- tmp[c(indx1,indx2),] + for(fd in 1:nFolds){ + folds <- createFolds(tmp1$label,nFolds) + for(fold in folds){ + valids <- tmp1[fold,] + rownames(valids) <- rownames(tmp1)[fold] + + trains <- tmp1[setdiff(1:dim(tmp1)[1],fold),] + trains$label <- as.factor(trains$label) + + tmpRF <- randomForest(label ~ . ,data=trains,importance=T,ntree=1000,nodesize=5) + predicted <- predict(tmpRF,valids,type='prob') + prediction[rownames(predicted),] <- predicted+prediction[rownames(predicted),] + + } + } + }#while + colnames(prediction) <- colnames(predicted) + prediction <- data.frame(t(apply(prediction,1,function(v){v/sum(v)}))) + prediction$predicted <- as.character(apply(prediction,1,function(v){names(v)[v==max(v)]})) + prediction$observed <- M1$label + rtv <-sprintf("%s %4.3f",paste0(feature,collapse=","),sum(prediction$observed==prediction$predicted)/dim(M1)[1]*100) + #prediction + rtv + } + + +############################################################################################################### +drawVolcano <- function(df1,pdfPath,strTitle="volcano plot",outFile=NULL){ + label <- df1$label + lbls = unique(df1$label) + table(df1$label) # A143,C75 + fc <- apply(2^df1[,colnames(df1)!='label'],2,function(x){ + log2( mean(na.omit(x[label==lbls[1]])) /mean(na.omit(x[label==lbls[2]]))) + }) + + df1[is.na(df1)] <- 0 + pValue <- apply(df1[,colnames(df1)!='label'], 2, function(v) { + p1 <- t.test(v[label == lbls[1]], v[label == lbls[2]], paired = F, var.equal = F)$p.value + p1 <- p.adjust(p1,method="BH") + p1 + }) + + pdf(pdfPath) + plot(fc, -log10(pValue), col = '#00000033', pch = 19,xlab = 'log2(FC)', ylab = '-log10(p-value)', main = strTitle) + abline(h = 1.3, v = c(-log2(1.5),log2(1.5)), lty = 2, lwd = 1) + + up <- fc >= log2(1.5) & pValue <= 0.05 + points(fc[up], -log10(pValue[up]), col = 1,bg = brewer.pal(9,"YlOrRd")[6], pch = 21, cex = 2) + + down <- fc <= -log2(1.5) & pValue <= 0.05 + points(fc[down], -log10(pValue[down]), col = 1,bg = brewer.pal(11,"RdBu")[9], pch = 21, cex = 2) + dev.off() + + name = list(up = data.frame(prot = colnames(df1[,colnames(df1)!='label'])[up],fc = fc[up], p_value = pValue[up], + type = rep("Upregulation",sum(up)),stringsAsFactors=F), + down = data.frame(prot = colnames(df1[,colnames(df1)!='label'])[down],fc = fc[down], p_value = pValue[down], + type = rep("Downregulation",sum(down))),stringsAsFactors=F) + name1 = rbind(name[[1]],name[[2]]) + rownames(name1) <- 1:dim(name1)[1] + + if( !is.null(outFile) ){ + write.table(name1,file=outFile,sep="\t",col.names=T,row.names=F,quote=F) + } + #write.xlsx(name1, "Table_1_TPD_volcano_prot_AC10_fc1.5_190131.xlsx",showNA = T,row.names = F) + sum(up) + sum(down) + return(name1) +} + +drawUMAP <- function(M1,ptColors, strTitle="UMAP",rowNormalization=T,colNormalization=F){ + + if(!'label' %in% colnames(M1)){ + print('A column with named label must existed in data frame') + return(NULL) + } + + tmp <- M1[,colnames(M1)!='label'] + if(rowNormalization){ + tmp <- data.frame(t(apply(tmp,1,function(v){(v-mean(v,na.rm=T))/sd(v,na.rm=T)})),stringsAsFactors=F) + rownames(tmp) <- rownames(M1) + } + if(colNormalization){ tmp <- apply(tmp,2,function(v){(v-mean(v,na.rm=T))/sd(v,na.rm=T)}) } + tmp[is.na(tmp)] <- 0 + + obj = umap(d=tmp,method='naive') + clnames <- colnames(tmp) + df1 <- data.frame(obj$layout) + df1$label <- M1$label + colnames(df1) <- c('X','Y','label') + + p <- ggplot(df1, aes(x=X, y=Y, colour=label)) + geom_point(size=4) + p <- p + theme( panel.grid.major = element_blank(), + panel.grid.minor = element_blank(), + panel.border = element_blank(), + plot.title = element_text(size=16), + axis.line.x = element_line(color="black", size = 0.5), + axis.line.y = element_line(color="black", size = 0.5), + + panel.background = element_blank()) + + #strLabx <- sprintf("PC1(%4.2f%%)",percentages[1]*100) + p <- p + labs(title =strTitle) + p <- p + scale_colour_manual(values=ptColors) + p + } + +Volcano <- function(df1,pdfPath,outFile=NULL,thresholdFC=1.5,thresholdPValue=0.05,strTitle="volcano plot"){ + label <- df1$label + lbls = unique(df1$label) + table(df1$label) # A143,C75 + fc <- apply(2^df1[,colnames(df1)!='label'],2,function(x){ + log2( mean(na.omit(x[label==lbls[1]])) /mean(na.omit(x[label==lbls[2]]))) + }) + + df1[is.na(df1)] <- 0 + pValue <- apply(df1[,colnames(df1)!='label'], 2, function(v) { + p1 <- t.test(v[label == lbls[1]], v[label == lbls[2]], paired = F, var.equal = F)$p.value + p1 <- p.adjust(p1,method="BH") + p1 + }) + if(!is.null(pdfPath)){ + pdf(pdfPath) + plot(fc, -log10(pValue), col = '#00000033', pch = 19,xlab = 'log2(FoldChange)', ylab = '-log10(P-value)', main = strTitle) + abline(h = -log10(thresholdPValue), v = c(-log2(thresholdFC),log2(thresholdFC)), lty = 2, lwd = 1) + + up <- fc >= log2(thresholdFC) & pValue <= thresholdPValue + points(fc[up], -log10(pValue[up]), col = 1,bg = brewer.pal(9,"YlOrRd")[6], pch = 21, cex = 2) + + down <- fc <= -log2(thresholdFC) & pValue <= thresholdPValue + points(fc[down], -log10(pValue[down]), col = 1,bg = brewer.pal(11,"RdBu")[9], pch = 21, cex = 2) + dev.off() + } + + name = list(up = data.frame(prot = colnames(df1[,colnames(df1)!='label'])[up],fc = fc[up], p_value = pValue[up], + type = rep("Upregulation",sum(up)),stringsAsFactors=F), + down = data.frame(prot = colnames(df1[,colnames(df1)!='label'])[down],fc = fc[down], p_value = pValue[down], + type = rep("Downregulation",sum(down))),stringsAsFactors=F) + name1 = rbind(name[[1]],name[[2]],stringsAsFactors=F) + rownames(name1) <- 1:dim(name1)[1] + name1 <- name1[order(abs(name1$fc),decreasing=T),] + if( !is.null(outFile) ){ + write.table(name1,file=outFile,sep="\t",col.names=T,row.names=F,quote=F) + } + #write.xlsx(name1, "Table_1_TPD_volcano_prot_AC10_fc1.5_190131.xlsx",showNA = T,row.names = F) + sum(up) + sum(down) + return(name1) + } + +RFScore1 <- function(M1,nFolds=10,itors=10){ + if(! 'label' %in% colnames(M1)){ + print('A column with name "label" is required for labeling the class of each sample') + return(NULL) + } + clsses <- unique(M1$label) + #rownames(tmp) <- sapply(1:dim(tmp)[1],function(v){paste0('R',v)}) + sampleNumber <- sapply(clsses,function(v){sum(M1$label==v)}) + names(sampleNumber) <- clsses + maxLabel <- names(sampleNumber)[max(sampleNumber)==sampleNumber] + minLabel <- names(sampleNumber)[min(sampleNumber)==sampleNumber] + size <- min(sampleNumber) + + probs <- rep(0,max(sampleNumber)) + names(probs) <- rownames(M1)[M1$label==maxLabel] + + indx2 <- rownames(M1)[M1$label==minLabel] + + prediction <- matrix(0,nrow=dim(M1)[1],ncol=2) + rownames(prediction) <- rownames(M1) + + while(sum(probs==0)>0){ + indx1 <- names(sort(probs))[1:size] + probs[indx1] <- probs[indx1]+1 + + tmp1 <- M1[c(indx1,indx2),] + for(i in 1:itors){ + folds <- createFolds(tmp1$label,nFolds) + for(fold in folds){ + valids <- tmp1[fold,] + rownames(valids) <- rownames(tmp1)[fold] + trains <- tmp1[setdiff(1:dim(tmp1)[1],fold),] + trains$label <- as.factor(trains$label) + tmpRF <- randomForest(label ~ . ,data=trains,importance=T,ntree=1000,nodesize=5) + predicted <- predict(tmpRF,valids,type='prob') + prediction[rownames(predicted),] <- predicted+prediction[rownames(predicted),] + } + } + }#while + colnames(prediction) <- colnames(predicted) + prediction <- data.frame(t(apply(prediction,1,function(v){v/sum(v)}))) + prediction$predicted <- as.character(apply(prediction,1,function(v){names(v)[v==max(v)]})) + prediction$observed <- M1$label + feature <- colnames(M1)[colnames(M1) !='label'] + rtv <-sprintf("%s %4.3f",paste0(feature,collapse=","),sum(prediction$observed==prediction$predicted)/dim(M1)[1]*100) + print(rtv) + #prediction + rtv +} diff --git a/common_change_label.R b/common_change_label.R new file mode 100644 index 0000000..b51f5f7 --- /dev/null +++ b/common_change_label.R @@ -0,0 +1,231 @@ +# This file includes PCA UMAP TSNE and Vocano plot +# copyright: wuzhicheng@fudan.edu.cn +# +library("FactoMineR") +library("factoextra") +library("grid") +library("data.table") +library("ggplot2") +library(pheatmap); +library("RColorBrewer"); +library("tsne"); + + +# draw PCA plot + +# +#DF is a dataframe, rows are samples and columns are proteins, rownames are sample names and DF must have a column +#named label to indicate which class the corresponding row below to. +#ptColors?? point colors example as : c('B'='red','M'='black'). the number of colors equals with the number of distinct classes +# strTitle: title text showed on the plot. +# + +drawPCA<- function(DF,ptColors,rowNormalization=T,colNormalization=T,strTitle=NULL){ + M <- DF[,colnames(DF)!='label'] + if(rowNormalization){ + M <- data.frame(t(apply(M,1,function(v){(v-mean(v,na.rm=T))/sd(v,na.rm=T)}))) + #print('Normalization by row is done!') + } + if(colNormalization){ + M <- apply(M,2,function(v){(v-mean(v,na.rm=T))/sd(v,na.rm=T)}) + } + clnames <- colnames(DF)[colnames(DF)!='label'] + M[is.na(M)] <- 0 + m1 <- prcomp(M,colNormalization); + Y <- scale(M, m1$center, m1$scale) %*% m1$rotation + Y <- Y[,c(1,2)] + + Y <- data.frame(Y,DF$label); + colnames(Y) <- c("PC1","PC2","label") + if(is.null(strTitle)){ + strTitle <- sprintf("PCA:%d features",length(clnames)) + } + eigs <- m1$sdev^2 + percentages <- eigs[1:2] / sum(eigs) + # lblColors <- c(N='#537e35',M='#e17832',A='#f5b819',C='#5992c6',P='#282f89',W='mediumorchid3') + #lblColors <- c(training='#537e35',validation='#e17832',A='#f5b819',C='#5992c6',P='#282f89',W='mediumorchid3') + p <- ggplot(Y, aes(x=PC1, y=PC2, colour=label)) + geom_point(size=4) + p <- p + theme( panel.grid.major = element_blank(), + panel.grid.minor = element_blank(), + panel.border = element_blank(), + axis.line.x = element_line(color="black", size = 0.25), + axis.line.y = element_line(color="black", size = 0.25), + plot.title = element_text(size=16), + panel.background = element_blank()) + + strLabx <- sprintf("PC1(%4.2f%%)",percentages[1]*100) + p <- p + labs(x =strLabx,y = sprintf("PC2(%4.2f%%)",percentages[2]*100), + title =strTitle) + p <- p + scale_colour_manual(values=ptColors) + + p +} + +###tSNE plot DE: dataframe with rows as sample and columns as features,a column with name 'label' is required, represents the label of samples +# lblColors: a vector containing colors for each label, like lblColors=c(M = "forestgreen", N = "gray0", P="firebrick", C= "red",A="blue") +drawTSNE <- function(DF,ptColors,rowNormalization=F,colNormalization=F,perplexity=10,strTitle='tSNE'){ + M <- DF[,colnames(DF)!='label'] + if(rowNormalization){M <- data.frame(t(apply(M,1,function(v){(v-mean(v,na.rm=T))/sd(v,na.rm=T)})))} + if(colNormalization){M <- apply(M,2,function(v){(v-mean(v,na.rm=T))/sd(v,na.rm=T)})} + M[is.na(M)] <- 0 + indx <- match('label',colnames(M)) + clnames <- colnames(DF)[colnames(DF)!='label'] + tsn = tsne(M,perplexity =perplexity) + cnames <- colnames(M) + tsn <- data.frame(tsn,DF$label) + colnames(tsn) <- c("X","Y","label") + rownames(tsn) <- rownames(M) + #tsn <- tsn[-c(which(tsn$X==min(tsn$X)),which(tsn$Y==min(tsn$Y))),] + #tsn <- tsn[-c(which(tsn$X==max(tsn$X)),which(tsn$Y==max(tsn$Y))),] + #lblColors <- c(N='#537e35',M='#e17832',A='#f5b819',C='#5992c6',P='#282f89',W='mediumorchid3') + + #lblColors <- c(A='#537e35',M='#e17832',N='#f5b819',B='#5992c6',C='#282f89',W='mediumorchid3') + p <- ggplot(tsn, aes(x=X, y=Y, colour=label)) + geom_point(size=4) + p <- p + theme( panel.grid.major = element_blank(), + panel.grid.minor = element_blank(), + panel.border = element_blank(), + plot.title = element_text(size=15), + axis.line.x = element_line(color="black", size = 0.5), + axis.line.y = element_line(color="black", size = 0.5), + panel.background = element_blank()) + p <- p + labs(title=strTitle) + p <- p + scale_colour_manual(values=ptColors) + p +} + +# parse gene symbol from uniprot website +# acc: protein access Id like ' "O00264" "O00468" "O14773" "O14979" "O15230" "O43175" "O43707" +# return: the gene symbol for the acc +geneName <- function(acc){ + kk <- read.table(sprintf('https://www.uniprot.org/uniprot/%s.fasta',acc),sep="\n",stringsAsFactors=F,header=F) + a <- strsplit(kk[1,]," ")[[1]] + b <- a[grepl("GN=",a)] + strsplit(b,"=")[[1]][2] +} +############################################################################################################### +drawVolcano <- function(df1,pdfPath,strTitle="volcano plot",outFile=NULL){ + label <- df1$label + lbls = unique(df1$label) + table(df1$label) # A143,C75 + fc <- apply(2^df1[,colnames(df1)!='label'],2,function(x){ + log2( mean(na.omit(x[label==lbls[1]])) /mean(na.omit(x[label==lbls[2]]))) + }) + + df1[is.na(df1)] <- 0 + pValue <- apply(df1[,colnames(df1)!='label'], 2, function(v) { + p1 <- t.test(v[label == lbls[1]], v[label == lbls[2]], paired = F, var.equal = F)$p.value + p1 <- p.adjust(p1,method="BH") + p1 + }) + + pdf(pdfPath) + plot(fc, -log10(pValue), col = '#00000033', pch = 19,xlab = 'log2(FC)', ylab = '-log10(p-value)', main = strTitle) + abline(h = 1.3, v = c(-log2(1.5),log2(1.5)), lty = 2, lwd = 1) + + up <- fc >= log2(1.5) & pValue <= 0.05 + points(fc[up], -log10(pValue[up]), col = 1,bg = brewer.pal(9,"YlOrRd")[6], pch = 21, cex = 2) + + down <- fc <= -log2(1.5) & pValue <= 0.05 + points(fc[down], -log10(pValue[down]), col = 1,bg = brewer.pal(11,"RdBu")[9], pch = 21, cex = 2) + dev.off() + + name = list(up = data.frame(prot = colnames(df1[,colnames(df1)!='label'])[up],fc = fc[up], p_value = pValue[up], + type = rep("Upregulation",sum(up)),stringsAsFactors=F), + down = data.frame(prot = colnames(df1[,colnames(df1)!='label'])[down],fc = fc[down], p_value = pValue[down], + type = rep("Downregulation",sum(down))),stringsAsFactors=F) + name1 = rbind(name[[1]],name[[2]]) + rownames(name1) <- 1:dim(name1)[1] + + if( !is.null(outFile) ){ + write.table(name1,file=outFile,sep="\t",col.names=T,row.names=F,quote=F) + } + #write.xlsx(name1, "Table_1_TPD_volcano_prot_AC10_fc1.5_190131.xlsx",showNA = T,row.names = F) + sum(up) + sum(down) + return(name1) +} +# see drawPCA +drawUMAP <- function(M1,ptColors, strTitle="UMAP",rowNormalization=T,colNormalization=F){ + + if(!'label' %in% colnames(M1)){ + print('A column with named label must existed in data frame') + return(NULL) + } + + tmp <- M1[,colnames(M1)!='label'] + if(rowNormalization){ + tmp <- data.frame(t(apply(tmp,1,function(v){(v-mean(v,na.rm=T))/sd(v,na.rm=T)})),stringsAsFactors=F) + rownames(tmp) <- rownames(M1) + } + if(colNormalization){ tmp <- apply(tmp,2,function(v){(v-mean(v,na.rm=T))/sd(v,na.rm=T)}) } + tmp[is.na(tmp)] <- 0 + + obj = umap(d=tmp,method='naive',n_neighbors = 10) + clnames <- colnames(tmp) + df1 <- data.frame(obj$layout) + df1$label <- M1$label + colnames(df1) <- c('X','Y','label') + + p <- ggplot(df1, aes(x=X, y=Y, colour=label)) + geom_point(size=4) + p <- p + theme( panel.grid.major = element_blank(), + panel.grid.minor = element_blank(), + panel.border = element_blank(), + plot.title = element_text(size=16), + axis.line.x = element_line(color="black", size = 0.5), + axis.line.y = element_line(color="black", size = 0.5), + + panel.background = element_blank())+ + geom_text( + label = row.names(M1), + colour = "blue", + size = 2, + nudge_y = 0.1 + ) + + #strLabx <- sprintf("PC1(%4.2f%%)",percentages[1]*100) + p <- p + labs(title =strTitle) + p <- p + scale_colour_manual(values=ptColors) + p + } +Volcano <- function(df1,pdfPath,outFile=NULL,thresholdFC=1.5,thresholdPValue=0.05,strTitle="volcano plot"){ + label <- df1$label + lbls = unique(df1$label) + table(df1$label) # A143,C75 + fc <- apply(2^df1[,colnames(df1)!='label'],2,function(x){ + log2( mean(na.omit(x[label==lbls[1]])) /mean(na.omit(x[label==lbls[2]]))) + }) + + df1[is.na(df1)] <- 0 + pValue <- apply(df1[,colnames(df1)!='label'], 2, function(v) { + p1 <- t.test(v[label == lbls[1]], v[label == lbls[2]], paired = F, var.equal = F)$p.value + p1 <- p.adjust(p1,method="BH") + p1 + }) + if(!is.null(pdfPath)){ + pdf(pdfPath) + plot(fc, -log10(pValue), col = '#00000033', pch = 19,xlab = 'log2(FoldChange)', ylab = '-log10(P-value)', main = strTitle) + abline(h = -log10(thresholdPValue), v = c(-log2(thresholdFC),log2(thresholdFC)), lty = 2, lwd = 1) + + up <- fc >= log2(thresholdFC) & pValue <= thresholdPValue + points(fc[up], -log10(pValue[up]), col = 1,bg = brewer.pal(9,"YlOrRd")[6], pch = 21, cex = 2) + + down <- fc <= -log2(thresholdFC) & pValue <= thresholdPValue + points(fc[down], -log10(pValue[down]), col = 1,bg = brewer.pal(11,"RdBu")[9], pch = 21, cex = 2) + dev.off() + } + + name = list(up = data.frame(prot = colnames(df1[,colnames(df1)!='label'])[up],fc = fc[up], p_value = pValue[up], + type = rep("Upregulation",sum(up)),stringsAsFactors=F), + down = data.frame(prot = colnames(df1[,colnames(df1)!='label'])[down],fc = fc[down], p_value = pValue[down], + type = rep("Downregulation",sum(down))),stringsAsFactors=F) + name1 = rbind(name[[1]],name[[2]],stringsAsFactors=F) + rownames(name1) <- 1:dim(name1)[1] + name1 <- name1[order(abs(name1$fc),decreasing=T),] + if( !is.null(outFile) ){ + write.table(name1,file=outFile,sep="\t",col.names=T,row.names=F,quote=F) + } + #write.xlsx(name1, "Table_1_TPD_volcano_prot_AC10_fc1.5_190131.xlsx",showNA = T,row.names = F) + sum(up) + sum(down) + return(name1) + } From 87b667b97ce3c7babb854dc9015b519aca586262 Mon Sep 17 00:00:00 2001 From: "Hao.Chen" <11779668+ewail@users.noreply.github.com> Date: Fri, 5 Jun 2020 12:24:45 +0800 Subject: [PATCH 2/5] Update --- README.md | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/README.md b/README.md index 581ef15..9ad034a 100644 --- a/README.md +++ b/README.md @@ -15,3 +15,12 @@ Code for processing and analyzing metabonomics data ## datamining_library_ge20200306.R Common self-built function library + +## CVDSBA_UMAP.Rmd +Umap analysis and visualization + +## CVDSBA_heatmap.Rmd +Heatmap analysis and visualization + +## CVDSBA_mfuzz.Rmd +Mfuzz analysis and visualization From d0f2361971865e2365ee0687f8367c228053a6cb Mon Sep 17 00:00:00 2001 From: Allen188 <50605360+Allen188@users.noreply.github.com> Date: Wed, 15 Jun 2022 21:30:13 +0800 Subject: [PATCH 3/5] Add files via upload --- sTable1_20200329-2100.xlsx | Bin 0 -> 29761 bytes 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 sTable1_20200329-2100.xlsx diff --git a/sTable1_20200329-2100.xlsx b/sTable1_20200329-2100.xlsx new file mode 100644 index 0000000000000000000000000000000000000000..842e2a006319524bf1380454ae7cc11a7697144b GIT binary patch literal 29761 zcmeFY^M7RT)-4*_MklG*wr$(CZKGqW!;Wp+?AS&p-5sOTv2Ufn``!CJ_dVzS1LxKc ztJbG5GoLZXm}9P$r<7#DAuvFoKwv;XK!`#52bP*jKtVuop+P{GDiYFHo)FNj z721S0BB#J2kqsGd{kEB)u?WCfz!63`Gtm`%39w4a+OFm`>@lIh7A;~pt*B(mu+rD! zI61yr{n=F+Nv!po*{Fy^>ck@8%KWZfS}IKp+D2_Mr?!P+(F9s;EKFDfi`lcEWlA=; zzuGA3$%f~(R30HVm~zXs{P7H$;roc|ttPFzGsC7k;^o=%X(|Qzo&g0^^zUi`K=unk`~j#!H=Ih9ggu2(J@| zwXM#9KRwMX8#*JVJm;c3d+wg-50+hHtDd;VD=Mg|D|wyvMDJ1G8Z+Iz_2cP^?=^!I zmz;0&-10xIeCmJqFc$X-Wb?4NYk*E3QVA-by!OL{vI$iy6XHL+Z;^&7V%zJ0_~SGD4{t~dGW!b9S1 zXlB=cEtcJf{_01skTb4A!iM(1GAlVoo_`yt)68}?nd?2w2|H9Z{`*@<3{m+G(CGfB zNZjW(9~{6zKv;o}hzvx=%Z}00(Z$xp(b4w9{)*KH95%U7da&05QM@$CXwvQAm7}$y zmBFR)ig(PAq2P21=}6EjRMgJ$k7@*!@~P%FckbYyWX}Uek;o^bPL1|TL*%1G|6EP~ zp8kFGt8zS!mxh+YE9iG{9KwWY5>qeCT;5^>w~#~$jG`0+7vda$&D`VcCQ|)t}PUM{wz-+yOEZ#$=d{0Ju&UWA&0^vI-U%^(5VRdtVJ<97*~R#`^&j z>j@LQsCEV_P5Q`UZe^jq2`z1REAeh*J>j4B`dJNBGC>(lITg=P7`WnRidcG5EVBdW zm>zLyVR_JjcA`Jrcn9WYG$3P2%Ec;?KsIcoFCP-@;eXViRC}1#6RFH5-^3rYZq+C6 zuxfkb(_MlhZ1%;43+u#b_`&S3&Tk=dLK6^`(3P#$xFM>fV;9d+rje3Trt!#U6iNp@ z&0)%RnkTw*S85{_#sA>(HzRQT=|_#IBxv7Iid*Hd!XqgmqN;Le0|$25`Z@NhyWocn z1j>|jZTU~Jc_k;vC%xzAE*-o2A-SMGakNjh7Re%m!509ix#Ts<>yza2sd96|{PDb3 zc5ZL}=d0N=qQJnn`^~XS!`_*})zuUjWXl+xDMIQmXQkg5t-cIzIK%TfuD9$8N`gws z^erBiBt!Wr7P4{`asEV7KyXW{M)v`eWbu05Ip=+yJI%}D$B&rvQVLgi0&hxC^ev+L z)cF-V0EDsla4ckrMz>1H6sr7CfZ0t!V;iRRHH#eixvy_yFQLp`5jS9hI(GKTa38Jz0Z^r<=@y(L4iJ zmSs@ugE)O3GheAg>w|Uz42L^bB!^=tXNR){b}FUPoIB!5LpTYAe@{dXnU?Q^mq!Nc zj7f5zc{lwJdGGaV$mr}Ag`!Xg*-o1UKwlwsnp3t?euJs6mfsp(fDNi0;b~ad7mJvP zHIHx~4mp7;wQHzp*2wvBg4na60@J%uFT1c_tj@8Fy)9DS*lM<<6cvmv+z&dJqa-6kJ;U>z z`A8Jxqy++UY`QCc*w>=3ZXLy_5YS>Z1PLEr7PJF*U22c+{<)6XiEDPjPO&xWp`Rgzzgn8& zbJDDg6|Zqz1$m^Od?e8!zDlHV2#)G-xJ)lZ$}8*C{;l0=7%@3tI?=#bM}7rw(~XMd zh3042o+wkUUPHWn^5-u0bJ0%@Ia?8lh^Jmec6YoOeg9jAf$n=%^8%-%V<5wDfS>qR z^LMo}H+OSo{QJcG;Q~3S@{ao-xYC{n*85Q!&UxS`-GEf4i*{{ajBRtJ43o4rwFVb1 zzum{ngQpdyDWQ_h;0Oj-9>g7bGB%R(GHF<;&#}o+fsVM7tEtWz3SLqgQZ>X^M7!$< z3nLpH>pwm&`qfR3wt~b1gx5n_XKuYe>F$J>aV+ptP}@txA>6({K`>M?-t_hsmSvLz(@rpJZ$WnYMuPK4tk-9#;ut7C3A58Z&v z4c!IdXEHF4Z=YxnB`#rGiBtitmT&elW!}*;Uek5)cW3NO0r9SI$PjVCW}ap~K-BKQKLY-tcx)zsmaZ+ZeWy4XGHT@;< z@%(LaoiV%YuX_%l-8A zqoKk3b6?Mk`}580;cWKw#gv4XSAhGQcbl)Dv-hLlv$t#DwO@n(7T;Qq;%3=%cJAwJ zSsp_kQ62{;QMSmvAw$tOGY{`T*L$blgx~wtQnBY{ik@xRe!kDwuLs6$$bXt2U7DKC z1K<4aAN(Rc0(r@{UZ)oY0|QvRyqTjrnz1@x6Dk5)|WfS$S1-x4% zXKyy`em#39d47LTl256%>Fn$B-w<#VJZG@`Iq9G6et-DodNO;HBl+SzpNXQwZ{x*9 z)3D#P+)6KHsyPmmAasN@bb7<+)``F8;ql@})%AOByWsA6pZsgWH&1lc_o3(Oz{5I2 z4!*~|$ig)gHUAl?WxCKr(cLJ&z9?a{qC9{f#Ah?H?$I;uwT@7<$|rCl zXE*H%^66x;IeTo5xm;V{dt8rT`IttW9e8-&PaZ(S`?|mV@MzoIS$&*Vw6aVQ;J1Ho z@wVDnXDho9CQuEz7rqnT-B=NLS6rNII~dVrm2R}=sUp>D)xUjzojdoOzMMhY9h)x8 z&g(mSGkrdFalKzW?A-n)&enomY;hBl(07=9sQawiNYEX5(%!bo$PxK^xc;ioPrhb2 zIPEyS+qdJt<=-0WnpD#Ia(!{{^J4dQK0Di5^xiUSU#qCQfh4WE@pdkgePp7cDOleA zOJ0@FIwnHR!PZ)4vv0E$xz^_ZZnF=WUJ*>^Yv=ymrRbG)fBU+2n|0^wv(i@N#(9tj z!t)kpTOd*NY*x&LL{C?9gJ9Gvd)xbhM{kcA|7YF7wZMYdW9q>C*efpd1y5H%+*NJY zSld%s)1lWj_=KYeXPu#deIS3+wd1R4Q=gY?8MvGli!Xf4|Oy2+4ez}ym>EwcqSUZeCb<~cgcVp8b(``gTt|8VCIfu&6AqK z91=sjRxsU%9e(_392?uTWqOwHIY3KJ{r`s|^ ztkrv(Y3$V9dVN)8uF*g#QDG5UVA)7kJW#IS1Txqx$h(Re>LY%{q^EoT(r~!359P$= zSGq^Ro;TeN1MI-py;;MasNUD7%LAfk{Fn8-ACA3mi<*uGe(tySFXvsHzapv}4GaP* zL*MQXwByh2JGG_4e-S2g2)2>Ul+HEUPUV5^F%f<9n7X{2zirDrQ{D!Y*CW7lT_S2d zu0m418R)S-=-0hYQFV+{SY(yl9qztxOQ1otq8~Vctr`KCVn;)SFoM}f%}^+JO=C0s zFhh}paYDt(qX(|EpiWz7_tNXPMB@>{YxPY9QVnl^VvqUd4RlMU4k#JpVkI)+fN8wB z3q!@&hO^e?;Z8FgIH#8`9-H;PFTO26hSG~WVH%x@$kBuIlpXLy>eD^Ay7C2ASPBcz z-M1kvNpJbV992Rx`z*Hh>)rE1f)ac|TKdB-+)H8i@s|6?# zJHT*RCE=7|+$U2=(8o}#TuponX-r{xhAhh|_>R@JW;#q@I312>$`YB<)VXB07J_3a zMsLHD6Q7IZ{^PK|QU8Ki(nt7>{imX_f|jws&RoLdo7=eWaOn`PnpQaECnr%s1hgNu zd0&K9SN5YN=Ed2vH_Ds^5OtDYudM5l3QL9cid0Y(LXnpE5s^UkTL2Qa8(29|)Gz@i zvd}y|ztuqWyh)6pSe1i=FtD7=iSliHRT3am0~lG#=K5O3;~HbwIILMF4^qbE6!^BZ z*v%f?v_g&TVOEWW0X7^%JXhpR}Pqb&UFH zg;p~W1>}q{@8%(Te>vJmp`jZBOI}CO(R##&y8IbHxC66tZtloX(a^MK=g@Lw8amR` z?I8lSJ88?5i2a!mDj4_sIU)G)FeN_RpyjoweQXa)P@yA7dT(chMM&cu;o<`4Xg2rH zs%-9qq3;;LTK9I#ruKTV3PRkt(dluizZsE|zII~}X}-vHHRi~bS@3*Ym>I*Q{eyI9 zjtFHi7|94DCSQqXsx&@@3w6@TV}wI1?kzs?sHKu@QNYo@WvZ~-!Azu?kti67Z|gjD zFa52HEi%;@^fo72c4eZEshti58!^DqR**HzU@Fg7;YVYum8EWcxFxvtA+$C;-XOdY0pVVSSE*1qSUh#+=K=OVuYFLMxobOUe( zPkPpkG!JX$O%70~J(Hfa{9(kk=}+dXIf}I9d6Rw5oaO2Fnt0Dt>kbv?m-<& zuoY$sTOKaT5-r{mKQVRxB2SKAq|LukPSV4MuBJe6L}fp~I>Zj9gsOE53Wp4!8t0YJ zjze}UAcn2llhHW0Qr39JR+82$EKh4+uriTveRw#c-E2F;!z2tj70g3n?Nn-ElPnRGRot;`7M`ui$}SBZ{)t7} z^HK9$X+=*G(vjXJTbd_&al8qJJlEFIq+Uvy72eWv$pSWfYKA8<^Kb4&*Tfk5u8}@$ zxheCuNDINRPWHzLUO`;Wd;MB7zR54?7Tp$buIrUMBTVjmdDsld-<)(R8 z_TwN+E_O|`j|KVk!(HouDh-{2V-O0}XL-iXEHh=Me$<25Y}qBGaicHDCNGo)3v*E=AhZVxQ5At7pk=upe%2^ATu$q*u;xX@9aDH0<(MUYJX|@m%_EVRt9ki8 znAW(_K8?W7{*0!G;8}f8+uFO{Ho1_3xYBfKrI-*K#%y)MJ?AVBDtl}Tag@gL;PYUn z<%J5N7{Tk?E=MgYilj~byeieOF`Gj(t?dlaI2Wjof}l-} zu8Qc{aKNBKRMkZ%ypA({GG7WfE@5^OPe)p6P?6Y0CaGo$v0&%?Gv#Ef)WpN|6Gh{B z73THagX(mX^1|LQMz2fdYc^LV<$BbJ7={@FhgC&4xfKFhY$+?KHaPlDopU}RPg}v7 ziO5QkVKk0CLL(B6*O?v=6dqC``31P3Y(~ByFU}B8yTlw+JQtAM5%xi89kcyQcX}84 zTuN;$@rL8pt|EYGy)byyJPsLVlpE6}C}%z%%l9{W=EQ{s_&LAaxBKOgmwu^DFW9GG2+{Fs zaPjswI$&EE-}(hTF>1w<5a=A|1`r=vUk=DAbatf#OC?%^rs4}S+9CQ44#YRJ@YS}s z)Jo75?*7=_7c#pCKhs~AG_z!M#!vtBh-w_mZq6lI;@tS78EhEFI=Sl1knGjwRbegN z;29jJp|Vv-c7+sMkA7c4gg(*o{Qhin1ZNV4!Aa%~B7730xiN={fXLmCiLnavWuAl; zs>Ya}(b2=$lYY}((#Pb`1>+|iht)TWmgp4^Tw+$xJ%O)V>MNY$&!vFge1^@cb-oP@ zGH_!Ahjpj?qPgnQqoVh_N!u=T;krW5{PihNA6!-W4urPEDIG$*DPFeAu!jvxDsbu{!+sUjC~L z6W{0wXd$dPuG{tsL0FZhr-O+jrl$4l)>=fRJ?$d(x*&NbDHh#WV(xYG6IeU+B%m{c zQJJ5<UYTfTAsQYVd9ADj{;zQ8k`#GDLfbAJe7X7kTE$^Q2sF!!qwnN3ik&LpY1G2 zLdCW#acT)7%-7pRQ=6#5e^)xGa}zFBLO0k95IbOx59E3DDE9Lb4IW6=(u;f#Isv{% zv$T;^V+|4Fg-7#I=TKizySgiV^hTx+p?vnc`|OhNsB!VDz?QKjAZ%$^Z@7l9BW6j` z#+YuesGwYuQBdMJy!eX-#b?&SPvxSQAi5~X*V@Jq!2BVYYoNMO%Pmg&6T_CgcBFZA z>>_f^Je1>$$A0*toFH;9ll!+uSN9;FW;?W?^ z&dM!;7v|~qsb`3LWO6!JC7>S|e=;NFuH-8*tBpr5gIdR=lSm$u6`4HaqOmQm;ks8a zVToB>PVLepAKalwk_r1)Z1|*)*ptB6wG6e-AV6ZOT>2Vn664?)5f3Uh{I&6^r3V>f zUY*pXeX#us?ZWbHqpWs)*_W+qu5n(wx$guaC&0|7hk3+K>4DML=#b>a{?dLdFSLT! zh~x%IADU-{z<-BWV$SZM?tKJwhuGot{shBL05h2nFsP~wbPaw$^c2o)25QHh+|>G; z;{E6X*SLt?oKa_fgcOsHt{`}Y{VUSa*a-u&^2*0}AZkXsN6_9Ds67&_0KWv&4{ojj zGgo~OZusdt@jmIk@mh_l3*Uk+!RI9;5c0}E`K2r1gg{=&u1I-18A4gQ4bWwxxV$T1 zlwva5A@2ob4FfSBn`1Ck;*elbNUBzv%rqyx$=5}rVzF>heH)X4m#1C~No5OL-C#||77lsr)FlrXM1iTNcsCP&q%Jxa z%tLrHhcfI5&KyF)>Up^zG$d3ovrwg)_^38cjO_Equ9BRSS8*5J3cgKMdlxFD(Gh9t zIkj5h+Pe=B{p%j6AMUXaQz-jS<}MI_Gv^JlBO`HxASe%AekcU!U13T9QNPoFf!-mz zIl-tFz;x-~iSuYK_XCa`exDtIec;?(*(}-m6L-Pmu=2EXBI#Ow@Uj&w8iFyq>+>9% z8*rMTmm^JO?&KfvYNi~Q081$Q@Vh)+(9{Kl+CmHVdk950ySagBgCdA3a6Vx^=~4pv z37dDY{RU@S<)sbE(J2yi-vom{*@VItMlM{J^Sq*{+d?7pZWaKyNLT2S0$2R#j^IW@ zc4^9@Z6OXc=FVWCMw2x_^n3t9^XuRbTf9J#LRBh;oiVaZtWb&~Dg7E|A1sOJ4yvCp zm@LeERf)s=*}Yj{=9&nfkqFK{4+TWV!UVjK{R0`hi$+i>D-wL3U(LCfYo6Nz+cG zT^HchC_v}1?>B2lEN?~vx+*xX1Gvd}B>_DLcZ^(%vFJ~oP>l0^F52mEQ&4x)oFK-NIZy#*m3 zEmLdx)BFq?8}ewTs(OWwq2WOI`rlNY2Ss1a@(kn8OC;3*)2D9D^Zl~xm8K=0qrRC+XLOp~xIhbt%9CYVi? zuI$FEc=~Kd(c?=e60PO$$?4sI<4W|@4+>eK9sX6SpweHZB9jp?dGY*O2(DK9g=_1D z;Rat>=k2-#nw^3TM&6> zaO6+_ddf6h{s6al~0@i z$&eKgebnh+|Dd1n5Bicoi}ex#t|7i!f~Wzx{x3wRDM8o>R|HOV3Iit~M9K}Flah83 zjiudBS^48Ft}WY$_Rt%=pCsS_!U`rL2tsaXxDtz*5Ya~NUx2i6yKq4RqMs|t8&&;< znkNcKNnRDA597AHAfp({6}wGRuTu~{1ff)(@0gUfi*V_COMV;F`Kbd+)>PH~IIQa{ z%LA<)8RhZzTk8+xK%%9X32|ZFM`7LW2&l9`dQ#%%-3aZ8lMV3rXw8p>jz4bdO$T(@ zrZuWif0m(Ems!?ObliW7{N5pFtL*lBfo(iRwmB`lZwPV#B@0CnDQ*J@K~a?BF;$kE zugfmMNfdQ1$h2 zSYg~cadkO=kAPV36g9RK7QEm(lUE5@enZJ@c;}nXk{t`dxqA{44xf?7Xg>$!r;74$6`!j1jD={q|6^Cv*#-qFU2Sy zqEHV0c*dp0e@<4&^sm*;QlKVPAW>s>c+i^|GB*?bY|gIeQs~cy%49~wvkK|Tf(f4s z%V0j5`vn~6<^AXl=f#$atXVcr;0MsXSAOQi)8OO3SipmRl!&;n5uV zH@5L@t2rU`Eec>x6vNv>leEB`toFO`qbar>dF3U$@Qc4d-A z(UjO!m9bq(iq+_e&p+WWize6cB=f~1`c96X(@qtUz@ulo(*5%s+ph5q!c1fH_U>iskNu&9o!t9h}QbjK` zz@_F9N8CESgw>Pp6872|kIaJ8Dns6>2Jsv)QJ< z&*TDEO2;iVE1CYUm7bj3_$E~h8-MZr(?&VmY@Null^;ydz87Ky!<#XVz&EL1+$&uy z>gfl}>zG*wpQ4BD!F?MT<%UbDL(=Czc-`X#irplOSsmZ0a>9#(Z~ER1%#Q#j8oNi3t-O8 zru%y9H|u~`+`(a$59?p)a11gB^ctH;+PI(F2T3{Xf5ZXRSFXO77AySh1 zw~Df8S0Vnoe`;m%D=DhyG(#sG>nrqC{OHQz2G%HKLgmdZ2hOZ;3_L`;k_3 z9fYWr$gn0xi3@f`0)_s5m0fJ<0(n}zT9G9RoEi>iZDEl*V`QM13b<3PSC9 zGqsvbg-T_sj|}>i14zJ-8b%H0-Uc>B>78#^F5WAhl%M?^Q+ED?43PC=e&WhBAFN*> z?*(1}X2oI|N;4w**!;gkwT5A4o^&yitVY*Kvot_QW@z8rp^^2@G+Sg9{C@-@RB8ST zy5E05ZvwU%R=!kGj22m~Nr$e==+4>6T*#<8L_7<0zlua782QMD=t{K@0Ts*orNzy& z;nMde8*Ge0q8Jj4mqno&VwU~|eB>t1!&6jinf%)hA=b0T`zWIed=g|bPP5-U`vb;s+SJz+Qt4k*UTHGZ(+6y=S#;_O1 z{8xicdt3MY34K2dp@=A8(IQ*GFRVV>f0W6VS+U*BzbfuvYd}{xp zjmEzj=rIfAvbsl-Z5xw<1;s66qE>RTFo|s_l+Sl0%@E=n010ih;)m)J7?M6zZ`wRQ z;gh-F;Sbb=XjJuoG%-TMa%B|MGR}VY10#w*`2L}oVI8%&rc?(_e9uo9Xena4;?o~D zIM?=UtLl9MfDLJ#Sx1+=ry06w;5ruQDoPub|G}H@1$E$dTL*){N6rMVxQ4@uJ1XC) z_iyv4DLhF1XdYK-V&W+maW%(Q_8&SOGfV%_sdd%Ab4sN!F{QPs0y>HPAM3={jr?GW zqa_rC0@!&CW#0+?Vu7p&nZm-Ta>-eWq89Z9+2?KT9C;<_9*+?|cd>tY8 zgV_9&DP5GmB41OTrMtiygmq0tx|H$>@QxQSB^~9Me^Qz#%mR8+7BD4^1zFZ1QX61@ zjoR3&{4cvzv<7O*0mm?=Fw`;yFwMTKLM;%^_}7SPsU-gK*w4!Uc&v^qc@Gk`pU{@>0O_{Y@Fp0-&Hpcqjzo94~};3m`&Eqx3QJAsTY4Z5?UoEY8b( zbg!WYAil1m)HVOc*9{JJQ|BtkSh;y;g{OtVoyS=cHpgp`0$xQ2zquGzz?>dNg0=>s zI8_69bpRywOV-CYf{`4jp(70~7tR>Nn<$C??=gv-C+<)wMviaHl-?DZL8=i@@agX#0rB zOooK5Kp~z9b%g^xF>)?i=A&~p3}NeF@&aKdDEbV9`FY;e40@t+wKn7aca;iUn#y16 z-?kMihhVJA!IwNB0GZ5^jWU7>Ir%aCNq^a=`8zi_1LbglGec2cXV7Kx zW*3qnV*Ia(2xbJ+AqaRKwqVAQsqiDyb1Q<$p%G94EE`*h3he+&07?aB+5|Az56Le} zz+jKeIyS)CQwPXz1;1E}*_-|+otNK|V*Yj_7%lse{CnID&?WM`ySd{b#XywU2IQ|G zL@zje+2;f*+_&xT@&&Toq-z?PfpjbZ;?#u>UWSUHQK!!wh$k#pzyx&FjT3757gXil z9xUSJxd0Smac<~}K)TD5>` zeOZTK{yg~4#PBiq`Ku{ow3OAaj;T|+Dl&hGNd|PVFVaxM9l?hX$DmQPRa$2AmnRVE zy8vDOV9crm*gcAm*0>S{Z90 z{{NLy)#9Z%#RalK%$5SBBVlnZht+qe$|?xa_JJk!7mN^~V&fRxfj!#Qw7C;FnJtw8 zC$q!wjnRL1j#&qK64f3E2A8f5RI-Yh+Z|Q4b}R@{2fAmW(YM%62Rfu>SJUGKHFH}w zZVB7U;%td?XAeN$CmGd$lgIK@1vu-gVRXE5rURJWYuS}jKSEF-Q7mHgLx|~! zlWCp|ft|-Q<~Ui+-1LTyMMtHiM+|%xc!fYE{a-NDe4QZ?QTfb%t%Gr%Cmzzg)Ujh9 z9Ww#EnhHp?ud?#{=zxNq&1%4Z8zju@(43QBpaFXXk)DoO3)XXygi1fEtfdsA7*~%-zQxSt%^%AGS$f*I$c>i%dc1MLiUTvk8Uz z`aN7d;`j$_lG6P6xHZ*u3InWqtCk`uFz4wM;**Po<%?@E>GJ<-&V5HDVrN?FoVZE0 zf^Jz^-gqMy+mWO?ROc3GWTh8-$<*zmdY^>&)r=cR+*n-*<&u8k#7@&Fr4`0Zzu%qq zWe55_bg!RH{hv1qfYbIO*zWuRxyC6<(cr0=lP*!HG`n>6N$NPs#mn609nce44NcuW zb1>-hd+GLHJG*cQ^x}DL_FvNv^mBJO$u;g1bTG&jwA*zc3Rvgk-?b+?Nq`mDd1m~5q_-IywXr3JQV1^v*owNxBPBS9;&!` ze({@#;UZ%C;&7~KEPHHvF`_GQed65F_}Q(^du{v7?zV1*s7N{}b_W#^IT3G*%YykXW(=_1Y$4!4X_qkb_ z+nY1~J^tO-cc!fygU^W(z<48waDDN>@naIr@rr%Zj1yLyvLr!A_l~Lx7gu^G4mLPh zAWftoF})i@#_&l=2A-a+8TplQJ9%z0;V8-CENR81Haj_eTOTXoCAX-^HC`Z>yYpct zH^WCV1AZ_$-LO~b4}VVVS_au-dY~&86`P-;RB9}$IjqAYJl%_D>u$g)$wo42d;|2K zs5KFK!DN?Kb-HL&1XAbVX(zB`re4Kopm zwddD)#z4Q@-)+bO;evs_+>&{NVo?}2Mq-b>S)jThc)S91v|(D248L|z3+ahU@tb}S z_mcW2kPz?1q#G^>1~@;PpxT>@s5vCE|Acy?Oqqljjlkrm%p)WZjwWamP3z8t!vbKq-X)-v5!d(F=AtLDXy)|gD+e3U^X3yR%cxw5LmwJM{@;g)* zyW47yooXb)_Mv|9?8fdDFLmG6_Q9uO8~!9c05O_M8fkT|q+a8kL9;7F4rkkWFx!}^ zT$AM&H(aA}n!ZgmtnANHy_k9XQD_h)p=dH`zK)DNSUh1CPgHW~9Iq}jJmLLmr9lXU zFEk~>6>#s%j`@Pek%jcJ2-fDCMaCId%3)qr63FP zb-e3yv}u@sJI^<>nRvauEej013cf!?HhjDN-sE<&=w-a0=l49Hk{9To<)h_%5j~a{ z`2KV3l1T97XM!UkS92RT57V9~n(3bY$k+}d)({xM-4F0L6SN^1R9P_b{LO++MpAe$ z_uz8))qz%cM`t5k?B5G7t>&iz!-vUS7QJXASWXGd8^&A5L9{juLqkK}EbfXkhZGl3 zB9gToT~rNVTY0lV22rx+oP+6^tJE@NS^E&OivD9r%pTPj{vvdJ7fGm=A@ard0#`hd zT>30z&>AGaM0m-<+v)V6M9S%R&$j>Y#CB?I4@q*MwL%sm{fNlWap(>^KpSC56U9(n z-0PB(U?9E6_A}dm5)mmb%V4P!92IFdt990@L%pHxXS%4roENtv?I4R_$nA<`Mqf9L z*oGNlEt@8bk~{0#M^VA>Ud?nHxXC!5lHjUBW%;9P(V$IX3=!sz#AG3;W@VwTKfF$5(^81SK&kd%rK2ST{ z8GELuJ`0Prw&*_1BYuHqTyiGxwGr{|j9vkC7L1&EKby_)bGdE>d|6J4Ox^r^0c^KS zK3tfeGtaLq5CRF8drs>fLgLdimD$1EI(UbXOYb|Jv&l`yPDWvKAYST2(~XNh$GkAqSrn~!yy9Qx-90vHUDl+hp+CL!8SI;1H%9Gf zfh@VITg_L`ON|qELAT1TYFY5qvltv>ISYJccN}FyBoU0S7GVbKnOgI1TJuULUkY?2 zoPWVF*<)gQ{kr=u-t#L$ctbec;#2CODP~!7W_N5Y+)So~PYP>ACa97LBP$@mO?AD~ z+1tk5x{_bAn4UD&ttpOafhh!gkB)ab{aJgoGc#g|MX|%A1K)-9ByJVeR-B})J(z+* zb#9o~;jzUA<_j$;o9!d|5zf$<7gPu7k}h_kYP&zK^jC+6FUIM(*o_flsXTfq3KLMw zNg2MW!}NyqfoX_r^PdXsX*GjTcC>zE3DG)cIitAoj(f2TkXnat<0yYsAPU}7t#t~h z8ZjiSCR_hB*h3VT&;!nM$ca*{rNu}K9Um1G6Y6p@%^aLvkD2uA(euP;!fLfUdVxDN zo|Gp>$)hmC?CT(VRPEHfA}tY);(FFV!6^e@MkV-Nn&P%4(ei{Vl=XpS=?Klql3U90 zDeq6jsz0JjjWi_=uL|@owYCf$NNCI6GJHRl{Nup~nCSDSSmd?46^5dlnv^%2S?OR9 z6qh(2?A;}}_^%o%1;=w|blipLPGxc*o({IUmo1>%fHSC5AXk~nk?-B^aMviGuYSceoy%wjKPNPU9l(_CsGuyid(l+xaVw7&gBg+EA8cCSD|(l#prC zJd_+cpM;+cG)DU;RG-HEKrfa(?14W-{NLEI$XLu~8!QM&oID5!;y-rm>gH`{?)tIc z`AknYVRIBSfN}0!2<50~Dn?3LPYO<>OBNuc15pqEGcB2l-qrvj+NBv5j&1`Y20|(r zV3iVqC6C{JA!?t0E*jS?|1PK5`Yt|;yJO)-aITUSF&bi@w4ZZ5{Wz1u>KbnLT;u&+ zR&mSgMBG?%g|NKwXLzX9(dQrWFN=~MHRQ63!>|$YM4ufbdtJ)@=x!6e<@3we`rYr9 zoKm1g!jus_vPe@{lu86 z3QIW1>6kD?x>DkuLrZL_ zc>}Qnh%cR|@&<89vts)LbwQ&i+kj!IBm?%%Go?5)x{}pDZ_G2*jH=;zSoG+joFJMOrdJyi5e#=sdyM z23IGKly$wY$Oh~#V6%NPvnO65^||_3UPiur6@Yk*>8s^t4IkAp*2$ivovYC=-!25W z1$7t81n&4a^S+h+jz2I}OZv(rS97k5S*ej*D5HMr$4Ge3+h{MUHmBxJV2;USR65C{ zD!kp`mN{*DTjj8NZjttC(w*oMHx*xA=b&jdM(}jDc<|U@r`=I2?>}jBnWRSqty4mA zbCl73LL!eZDCv~+qqJVM-uz07S|0L?0^es`t*o(AD2}b%Bu_4Y<{!bVr^PkbmzJdW z_@>xgu&6C+IYYi(MhaD<%^SHk#KoE6I8pq7^v27CS_&pxIOFhGDNMVL>l5cQJ!_s! zzTfN_61|*a7(ZB9G=YC!`^TU_7r=+wtec7hrscfp05r_trn>1Q>8rNcFcqeUVcVW6 z^MWe6V45dUdF528_8FOTF`O&bpW@*m@~EeLYi$i#+h!g5$#UCp;n=NNz76<@i!B0= zpAzO4I^srPU{yPfcswA=aahXkY7p)qb8OS;722Jkfzti_d?A(gHWP)~B`7c^X(E6~ zb(YCGr)rL{CK+OsVVciEC{@#y3u0LBS%d%ih)PzaQ?-b_hLgd6eQSAwB5z%*<-fl5 zeXZJJk_P{+_%>-dbN%3vS^L|-Mg!4^=$hI$*{dwL)7Bfkw)3tc2#qM_mAt)ovaJ5T zS1M_v!-7phpBXMV0;1Ij?fTvalKt5O=Il7dJvsF~k!r3D;_ zs*c4m<zeQ=pUY>*uS5VE*Fw8;z_Ws5$RE^wdz*NXqfmuPx=IJ)^s#YSMqyu zXq(RG+E@P0{u)ojV*2Y7z8bthUjd%e6XNo|LqLKsUq#BM}Lw6)SP{*@QnCH$s3>q;qPz-K!f2;CvAV!cHv7Yg+E z6UwIi%TbQ?F3OE@Z94CCN&e&?6F*dU@gss|PxxGD4z{^uwuD74=mBYeI^k9qs0jJ@#A2{BJ}FjOmGtunV3BX2Tx9tNUGjK_Ht|9fwp)BnNGnB|*p{g+OSt zKY_NF!9Xy<&!Pr}ro({DVAZwHas@@X6rlO#tHC}Ym!O%*?Z`^}0`={nehGOiLhy`> zZu@>zlDYR?$zlh?MEo;OvN zxES9AH4{cuIEa&&C`xJo`vf+VP83IIK7!mcX-Y^NbQBGx2pa^2AR#K)If7g-nkf{8 zVPr5sK;WcspY#_(y*mRB%Xj^ikkAx{e_lIB*at6Wpb=7;G{3r1#>OPFwnQNUobFR5 zx?f#21blxepMzX%_LEUpU_VtzQ?(Y@&=x77!JHI3tXVWiS4SR*aLnnJqtTaWxCgDK zydus4U6S9v*~@&is&J^gGzKr`!!xw7LiIt=c|Gh}^L23T&u3vHy z%~42O%16+HB4LXgH8Q)H_Mi~3|H)i6SC>pIOF1Cc4wa9yjO&LAZRk1a_Y0I&v?;9w ze6i(~&12ufDbeTmdz^7X+*))b0_ffZ@wer7Pr7{t4iZsifA5BNdgN04-ZXha4#*3F z46iA5W)+%nmFkN|f^|a0OX2}->M7BRS3eKVzKj-`2O}1~>Y6vHO6j@)&)27vK53)h zshBw?x0YX#J-zq+Ff=hj*D)L{PD&`7C|la#hRXWC*`=q#QW-NI+Z1K8(prj^NYN27FQPN*N+3z2|e99grVrP!4Tv^&#@$Wfdh%DJemt zJ{->W^6RmS*0{#x+0~0H+_~wf8jat0=et)+;WhM^0uBqPoNx+76ix^4U&rj!QnrnA*zn z74;egIejy*dmM3+rN-6c;O^Jt>a4(odHqn@mQZWXi9_BSZ(mj%hF0`1mR&-(IXF}rUiiK9B-ex}{MD7NJVc9PkATjVM_&Svb@eX)_}TxSl$hvVD!bA63@ z@Adn7A1`mJ6>-3_>}j~Q=ts?B|xc0qb2kxl9?YOI9vP8;a7L?+p|jT=Lxy&w0jfm<&xA z^qh10>6OO%OgUff?~}*Vzw1_w=k|@rfc`|-PGY@IksL`N=y-WE{lf=QpRIpedru}4 z85yq^K*`7xk^=z=N0r^n`Z)sod!RLzFTvwbjKLyLU;V9)?0Jig#5TKa`Z{l2=BW1N z%T}wJKL0cK4Z7i#JD6ii=U8^rcvWwHa}6$~D$!0Ire*GQzaZVLyV(igv*HB<^ctIo ztYdjuY3a&YXXSEi_+n)Xo_mu;c2>X+lMi;Oy~^4&Z?CbeIvc&eCUyv?A{wKUKjT!# z1a9Fhn6y)S3A)Tom`bl?&9fkIe*2#{T-=v*vN0e*Kr{(KKrsF_;aM5Gn477(xmY_` z{vGs+wdEZ*HBh>-j{*T+Hb{tlKkI~45}Blaf+$3VLu%=RN%*v8$`s28KsyKcduSqk zw;7FymW)du%FXmwEZ({K)pkztaPm}eyrG@JmoU`dQsuDS%-Qt2?$7z0PLCs@bxDI>ZBj%wH8w-;M?vVtZ<(noA*<*8WxhD~l{BAN_rY;NzW- zQ+JfH=T*n_x!4dUY&a2KxIsHA7({BA85${FFvk|OvSFp)t4cI*mLoP5bUD~nC4)}3 ztvWP_hodY(m_}5~;n~g03&Tz!)#(N(5TvjBAW{UV!F#dfU19|xL`1CLVz)~q%X8O&DEFP+H!$-`SufEo%8CWI|U zagI^EhR85n2tS<26kLe#3ZIxg)pD4#F1`BaBwhOD&kew{@Oj`__?0dTA0r+PQZTcT z*?T=R3W_$Q*bqceI{n)PlO(5gA^51Dhy2-fi%<(^^sG>uS!N&X4yIw(n0&!5In zq^2R|j!mjqBh?8f_T0L=nSZokq!<~IfE$C{A1BuThzbx<7&c_|a}HD(NgOs}z)#t5 z4)mLePVR^wJb;XaICl6uB^LW<+^ATXQ!(Slgd?UzC`gMfsB6T~#3L*8Dujb)7~=WY zghZUFHf#|R(?KbOgRU8LvSDFdqEUM$qP(yK#lb4oPJHt)v9X9mDmf`P08G?iu}%u* zGL-2q96Fl(*ZwT>J&p{Csmf)?;- z{Rdv=!|pY+bH3NvnKQFQ z{PuogZY5p^9-y+Gy0NJrlAa({F)A&`f_kmaO0NXY(~HCe_7d&sQtTC6J#M$l(2?1RYi)T=xzRlzK5zs(mXj^@ zzdW(jT7SxRx;6dts#t-$f9}<`hd!9e5aPA1VOh5})vWxkan`_IyU*fDfxABGaB(S7 zbD>U0Tv&5)J2^{-Jp86ckSjxkJ5bbXRik@Jn#81 zAOirF5B_^Z*uYlL@n!;L|NG)E-L*!i6>Jmk6w#^x&iFe5X;=E92R#j;?{TQvzsGKy z27L2LmH3o~|J-Y=%&kx5>C{Q}KE;a+`U%N6d(PMri9EzGrihr(i}^~mNtbd7o=OAe zz|2)%X!HsWBRom#+i?eet)NM5KZCfBVYW=_Bd@=8Ih9{uIl2VCawz2EkR zg!Sh(CqAt`dZwFjA_|OT{_-(OW?zx|tgc**>p4(N+P%q?U~A;6tb9sYwl9w{K#tEJ^~0V~(**i-5hiZCj}s)yl~ zAFF$^`*Ws26hU#@s0`(6NAXtIDP!+m9<^=@fsMIyo|XxanYmTt3l5Mk5R;Kvi?K27@M^_)&8)W|Rj%`Hsn+^9&@&+>2IM9#>>wB0 zb(ys9?`cv(?`0KcmE}(nqhhSUfAaNsHfSNq{*0^H#-@k19j9w=|1qK&FNdZ(HeC9w z({2L9xXwvh_aU9dc|ZwLz$Xa;k7D@@V}|g=Hyq^ygL@6i38XVU-WPfpoQ0vA7e`$b2@NA}L?prKg|lu}XXNi5`Rin{ zq`~ot$hSm2r_U5NdVjR8Gc$1} zA59#~n5Hu4oeMhEGyo(Mp+u2A0PGm@f**6WdmCXa#L8+`)(iK9sV}Z@XzdNNX1PCB zmEZIV*bQQKh9WvRmDar8`MKT}hbCV_yfLnmT#%sGA+x$d+V0lCoz&*3#td!N^IIj! z`!0{=Co1keuUm|~=kNGoU*Y^bctJ)I)0=c*l#M%UmdM`EfLqdLQc>t^Kh6j70Y-|2 zE*q{JXNxpCKG(^<6imuX*cdQP#H@vGXh zQFYg}q)|{&{4-6&M@Rrzo;}5ts6Y)kv+q2ZEM>A@#eyrdgcZiYR%ha&{37J2`N4I( zofn0I)};BGF5$spl&fkba?9>jMCgHCj0by1pLv>|sIosmEw>~RjODwRLSYT}z!n#0 zQGKQI^C|J=n}t!$T3b=A!Zkx6PCD7@uAgW{&~nAGnBKzPDf*|2`O~JyxsS=Tky?B`4AzV~{3av^4s;+nyr1F(?!UJBq)Bec&7%nVLv zexPKO()}`3CE~y-px(u4ZvAkwtEE`1jXDXp4$aO8u9xH9{r6|aZ^xl~dVY{l*+qzv z#58sGxAz@2ZkZRMJkcxJ@XW`mo8US7QLnbQpV~P`ct5kJyES^xA2+N&5wfY)b%>vJ z0FzgP;EtK4olMg!MieK>=a)Cu!pzL}_m8dGwyKTGdM2YLu3-kxe=pGQfBDQ;3tN~! zg88PP`)h%I<4WgmcPD@O*7@c1vDJVxPI#7%=q~dST(rT3 zmnGh}y=)MO&C_$%(m3YIC~I6GTVaxEe>Txqk&wEb1cWkTTEX4+lMwNNJ>kn8iDns$ z3+Kw59`3h4Ggk7$5hh-70 z=EeI$@ilZa1Wd(em0sT?vES=0`7BJAKTYCUyV%yg`&*$6$Kxb(FcxB#Wmo6HF-+0^ z@BW`6E4@k007#*`rq)qy`im@o`sC5mAM`J@44GWMv#^^XrUeV19sT8;OZLI zsCH}yKe!Cb_;ngFy(8Q)p4NWWr%LlV`u?zzp6z;hY)SLcg2LsR*rd2^RU0ESGl{kB z70`wDp$579#J=@Rpupr|+Wc3h}Oo~liyffTn#Q)j((A+`VliP$eP@_;0n&_T!?|LEve@x z%k_%cOz)g`ba&^}zX4&_@-}U`;^P_A^Fz;$maeiU>+^dj^j zD-_tu83palvT|NI-FL2^BJWf>JG=A{KI-gEE{~tDYBF`x7~!e@u-x!#EwC3fCvtK@B$wo9h%?rmmQSM8lZrH!3chuZ2Lu>DX= zKlWN^Ep)Tc&F6zw@(5L+-C(b_*nP(A)u#mK>QMHyIPGZJqKV3Sw*Z&^vpPp!s5!Ly^A>aO4b*9s@17epIVXJ;-YOxP>>l^0AOXT;=GRTudE+I1+R*J;9I?c<{_HQ{I_<>Cu}q4-k$!=bE&PgO zSs0s?i3MPN6({ZA0r$_}9q&yKul{}SfvAgjGBYKzGs=D0(MP{C*%(=#{-U}B;hb$N zoy=ijeu4vBRE#$*5N2wy0Q#(nGJ?h-fjhy0UgOl4*R1N24%0LF3wq!is3ikT&J z?|xZ!uXWOYPxhH*>QBwXS4W;Q`_y3_)T3#VOTOv6tkK%;T?jalAztT3xh@@(FNrw- z#t}~8r$2&Lt-P&UAs7dAv_#Y&`ZJpFtC9K)$6$3SvKT>`x^_`q9Tt)-V99(k{3a3^ z-`Duz8zYLEF9~HCS%8!&>^89|)U+N~b9@}vb0Kk~w8R@&-qy(SdX%Dzo^;#u@V-^( zre+=@oCg5WDEZWoG1#q1Y@6~WSZXEm<8f~g$8f{ZIoB}#c+X>GSwkd&j@2xZkKp_8 zK~fkm^a|;t%jk!}Xmc2h%?SbRNU?|0fYkt$oC*JYrsy1T{cm)_+jQ#AO&C%v6jJ$` z1UG=O1Ta7}i>8S(i?Jd2r$R6e8QGcK<*hOs+PFCmnC%&cC8r9MYiB7RbbayH3X|NC zux}n2nm2kn%KFYJx)8Mz#rkvQJa7U}dIoQ!)#s*pirHw13YnUE6FI10P7A ze*Q^NiGs05`9KfZx9~e;h(hKkI;D=^GvaV`ja})L4YFwqsf2wh4q%b*BF#L^aDA?p`E)az}jvXpqui*3u$@Rl0AcClxr_mHZrIxx| z4_RW-ll!8LdY!Yx!p!C^JY>@uXYgZ+K)W!Y&Yu$@N&+Ej71`gqm6p`VU6ido<0+Ny zx+CCQmU@lli94?DK(;QRhn`C?fLk!Y${p1hNb*-T#|P)d1zCjeYr>rXaZ@n_{*r)M znO|@rk`C*=VmJ>s^+HYKl&#h3EZ`crTtTVIrLP^t+Hu{vIxt?*V!lcvqrT}~66$dD z5J`}`a$2-3t2Oo~2RK!2cXBwM$SwfG0DwVvLdWGEkg*vf3??i!qcI*Doe}fgepBZz zOs~-atr>2D2*Z-eo>Q{Wbu=!)TyDYKx-oik%1Bz&xGC7!iu{53Y!~@sIeQfoNbFUB zMk&C|n_TRcTPcJn&y)9NohhM`xQLt`b{-41_lAQil&K;p0lo(4DQwt}2s{qHZ^|-n z$$s7;L^BSm0J_lGrdCd5G&v8z=O8d20InLFKR*(l{^)+f!{i2%8IM<CYuEX;w;P>+m^K5RnP&4h!FZ7dOKjQ zJutUUlRlx*0OEyPNgX-;Vm0)Byhaem8l5}G8k{c7v#@c*Gl()TL|K6C0j%069r?{L zz4C)b*eap(C$GZiYXE?Zja-y}!=1Ozd;Isd&L{|wLy(Wh)LRtYpC3pa*&7~p-bnW? zgBXVkwl_U_Gf+6!wtQR%28ZCC+N%J0go19S0Ju>IIS@re?efqj|G~ zXQm8jlN#zK5bVD|V7-po;Pe)1Kz-FAc=EuznNAMTgU_e-?Qn)|wU`QKG$laJR1+|j zQUW+<_6HG-BAJi12vbw&6{;*vhQurNAbKs}WT|BqacX2$NnrAi|n1a+Md& zf13shho<{&105(^WQ|{91_y}bM?$b>jP)XsK;i16z<|w62=%D=r?QYhul`cYv7veY zflrl40R4e@dhJ+E+^?M25UP)hrRA?W#~FCYwF!_l;2JYn=c-Yq+I@^JOls)c5!!AP-LYu= zo>gF06U&?X2gV>(csEPrTJNrgTZ0-ZDh!BGkTzar*&hu@yCajLCi%^gNM07I7eTvk z5fk^@d!lM;hU<9BM_z`$l40#e-#X|`>{*413WdDCKh^g3u!vJ^88aZ|t{SG^9e+X- z%-Jj|al-^HcNJl(x*E+56S*+D=niEQ`vT!vAN{)}mO8K!WT?1c1qvz0{j^8V1{?Ff zsR-~4HY?vEF>_a`riwZWd2^3qDA9xmQyuRI^oYWk5T1Q3Crqm18vuh3e`Z7rOPXy`nl$iWhJ^oi=5w2z^p^3l)^FM{ZmUokQ!1w zmC_1Cr&o+EkMHP(Y_ysc%9DGNJoYsWl`tG`FEQ>C(|+g1IL_QyXZ5Bvp7~>%vM$o?^bQDZ)D|#_+Sa>>134e6rS&%rWpu!#isx zWK0+pO&WyY=63UOg=VLoAO*}=W^v=^a?e8Ct;NIj6zObVJH|70i!x*K&`~ES=tGA5 zwU(!h6uH6`5$Y)AdD7c~1-O+gw6Q}n8AI+1a7{O@5!`Nw*HL0bTZ&AMh12tp6UIAR zpakO+e5NO}1~A-UJu?`ts9Y?}gUTEF#!A7T=&TX6s2Xk-#F}ndGM1}S1i5*#KhX(E zvYyb!*pz4A6b}|@yCc9g&CxJ3WinUyrt`) zkBrZaBHIm?yqW1bU`!EqEE+K{A%7qVhkZ9be<glF~1NGTEJYErA5gUC)>fpcA0(>KX0q;bZyUUg#fp zNJEoosv2OMwHj~$0GPJ<#{H=^%u)w?B2+T8x4+dm-^p{Zy*0JkRGYFQ6C>cCdAcxG zL~;H7-9ScnZ>y!SrpZ8wBwD5b8(Sva_y8@vD9FF$z83WgE~%RR_u) zz=2B$iEb&qvX;Ggn^yrdi=iLSuFwVd_bV!v@jXQMX_$t3Q@4EcDdxyOz%$bhhP_Ds zG&;R*;b_0V;_N6!0x#baaA_*=aZ!JxzHAGQHEeAgRiDMq9xY|(zNSiI z|D4BDCxTat#L_IqZ!>UIO(EYn(v100?{T39DnI+$eqY(>d%Ng)YN*507P_DTiC1Yv z9ODJYTAT9^qZOII<^dG9JB+bci{6bU+ifr>yLqcQ&EzC5s$RlP-DlR@dpePs?Sv;1NYHtH<*&7j0myLZz})Z8|B z+JZd6Ts-JV-|aeVzeOahD7k-aBnYeY%&YL-FcE?)!bDoM+GCwwW}GS8+0~w&)1E%) zK6#&v{g4qo6G4d%lQf9tu~tJYNsnxWq+^$?*kHd+LIZ{oho4jACW4qgZZHG4XyQ82 zY8%6dL%Io{LUWf|!51SAGw%1Ya$c`x5%yk>3TIz>+Jzsm?c6I|ZPn`V=~Fo=S{7|q zS&AAebma{;k_{R2I+V-sW)UDvd<@m_ zEa$z&kaf~J{;tpfroLbJQ1!}|o@p*tGOIWK)y`ly9g8z0Go$>3MnO&Y^aanP*P(SX zD#1Q;XMy#lpV?p@Ppx@tU*9FW2`XW`;Yp`;lweR*4^LLkrjFe^LtNYd?h7g@<|BhAD<%^T5p{ps#14FJhu+YIUX-I7_b|0 zv6GLD#?r;a)_&!~tSwibqruKI&Y*)g(>)_fZt-rlbwq(~g?U;B=~%@yy&?5(geNv8Nj*<{-_)%^>&~H2S;%#^zy?3#0oNq`TgA7EdW!lH zz{(9O`TAxmdH)y%6Bb+brJlztg?t-_u7UDysaZ(AL@(zMf0c# zyKj&B?#0{2oj@c$-cxW9K=vl*W=lvdDuRK@!Tr4tN94~JhmCzGB2*Y);GZr&e7Bc} zRMUOPS)^m5qNbS_{_f+Cj!Q(#CaUFMfn)Q(6kSwzm9BpyBo@l4z?e6aaQ z$^kSOOdYx`s`&XmAz@rqBc{!-69{VFw3-R_AqBX z!_bCtDi%E!5P+o!xQStk7!nuTag@|ZzR28=S(~J#(@wGm)9V$%FIqgM_gu(8ye+i` z(l!JrCp~$MqYm1`pTI#7fj0I)zijA`a;%6NJkImOL2D|)q93L3PSwU|?>Si?U25iJ ztfni>xub+hYKe#T_x$)=Q~jm~_mPBS-GvnfrOkmD~54lm%Rd zLG*c&hV#i8fdnv2p#H59|DYs)kqi?zM%XF{=XY__v$p;}X@iN`-+u`ay=G(_IITNA zb5-gAvth7;XZNf6hja1Gel?D z;l%T`c!dY?rM;vQt8X8EpL<9|f=DJ!6zx-?{^2gj!NI=>oqQT-CUo24Bk`jjbyh6z8T@2h?Crjmcl4{SNMR{G zE^lp)@^LWa^Ne6eazM^{Y-yw=TtFJft;>{!F(qGS$@0zIt8bnljiV!0Fz9EbO+86mP+$-t*Ls5+>Bnuzyxq z2f}{Qe>d#d{T)T)_p)1K&tE35z8BI6LO95?I=ceoN(^FFJB2uCaEOS3=m%%pyjru!g+!1J1sf_#kIV~bu^Jv84^7c;r%!{HEvlSrV>vGs%0 zXUafPrE64SDZ{7(vLpfDJ3+Hv4%hFyhM{!N401jjlvOGVor!*}$~YR+ks+)H^`{;s zR`%RWR}7MeFy(hEhxiJs+Lm;8q>v@A=Ss}KRqSj)uQ{Fl=2llvv#Rw@gA=oK57MtC;Wt9SX2Hrd%(f7!ZtVl zIe^}OM$3QR{}4{^mGnOW{u!g|KZG~$)37P&w}@SL3Gc=&`bCNWJM;Q?AdCM*F}h3q z&&VpjNCAMPd$+{@k2ouLaqh;?_zelRn*5uIzs<|P49j;>?gmBpg<=DHZ2beut%3F~ z%H89>zfh*QZ&B_Z{=Exu_khwb01LibfM3Uz?vmc!`TC1gPW&&@yZd18BHZ2e^$Q_N z`WE4jd%*6Z+%@q3LfKKcMfu;%{dWQH+RT3eZYlnM0DoD~?*iU6LH`2WQMv{EeOCH| zLHaJHZ0PwE` z?OpPJE=hhThl51}`!D(53zWOmcjx|J^uv~aXy3oa`R{H1>zCk_G!kqZA^?B~OPjC& KIRc=Y-2VWM6>mWR literal 0 HcmV?d00001 From df9bb18ed422b031477aa42d60b24aa4c358c8f1 Mon Sep 17 00:00:00 2001 From: June <68281140+xuezhangzhi@users.noreply.github.com> Date: Tue, 30 Apr 2024 11:00:45 +0800 Subject: [PATCH 4/5] Add files via upload --- data/weigang_mata_boxplot.csv | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) create mode 100644 data/weigang_mata_boxplot.csv diff --git a/data/weigang_mata_boxplot.csv b/data/weigang_mata_boxplot.csv new file mode 100644 index 0000000..4f7aef0 --- /dev/null +++ b/data/weigang_mata_boxplot.csv @@ -0,0 +1,17 @@ +meta +1-(1-enyl-palmitoyl)-2-arachidonoyl-GPC (P-16:0/20:4)* +arachidonate (20:4n6) +"bilirubin degradation product, C16H18N2O5 (1)" +"bilirubin degradation product, C16H18N2O5 (2)" +"bilirubin degradation product, C16H18N2O5 (3)" +"bilirubin degradation product, C16H18N2O5 (4)" +"bilirubin degradation product, C17H18N2O4 (1)" +"bilirubin degradation product, C17H18N2O4 (2)" +"bilirubin degradation product, C17H18N2O4 (3)" +"bilirubin degradation product, C17H20N2O5 (1)" +"bilirubin degradation product, C17H20N2O5 (2)" +choline +kynurenine +mannose +serotonin +sphingosine From ad9a83147008307c0d59ccf7463b2d7776a54451 Mon Sep 17 00:00:00 2001 From: June <68281140+xuezhangzhi@users.noreply.github.com> Date: Tue, 6 Aug 2024 15:17:01 +0800 Subject: [PATCH 5/5] Add files via upload --- proteomics_tmt_20200325.Rmd | 42 ++++++------------------------------- 1 file changed, 6 insertions(+), 36 deletions(-) diff --git a/proteomics_tmt_20200325.Rmd b/proteomics_tmt_20200325.Rmd index 9182256..899c98a 100644 --- a/proteomics_tmt_20200325.Rmd +++ b/proteomics_tmt_20200325.Rmd @@ -425,7 +425,7 @@ dev.off() boxplot ```{r} -prot2 <- c("P02743","P02776","P02775","P02765","P35542","O95445","P04004") +prot2 <- c("P01011","Q9UK55","P02741","Q86UD1","P13796","P0DJI8","P0DJI9") label <- info$Type[match(names(df2),info$TMT)] df4 <- df2[prot2,label!="NA"] label4 <- info$Type[match(names(df4),info$TMT)] @@ -435,48 +435,18 @@ for (i in 1:nrow(df.nor)) { nm <- row.names(df.nor)[i] pv <- c() data1 <- data.frame(value=as.numeric(df.nor[i,]),type=label4) - a1 <- 0 - b1 <- 0 - for (a in lav) { + for (a in lav) { for (b in lav) { - if(a>=b){ - next - }else{ + a_type <- which(label4 %in% a) b_type <- which(label4 %in% b) c <- t.test(df.nor[i,a_type],df.nor[i,b_type], paired = F, var.equal = F)$p.value - pv <- paste0(pv,"\n",a,"_",b,":",format(c,digits = 3, scientific = FALSE)) - } - b1=b1+1 + pv <- paste0(pv,paste0(" ",a,"_",b,":",format(c,digits = 3, scientific = FALSE))) + } - a1=a1+1 } - - plot.boxplot <- function(data,x,y,type,filename,title="boxplot"){ - a <- ggplot(data=data, aes(x =x, y =y ,color=type,group=type)) + - geom_jitter(alpha = 0.3,size=3) + - geom_boxplot(alpha = .5,size=1)+ - labs(x="sample",y="value",fill= "type")+ - ggtitle(title)+ - theme_bw() + - theme(panel.border = element_blank())+ - theme(axis.line = element_line(size=1, colour = "black")) + - theme(panel.grid =element_blank())+ - theme(axis.text = element_text(size = 15,colour = "black"),text = element_text(size = 15,colour = "black"))+ - theme(axis.text.x = element_text( hjust = 1,angle = 45))+ - scale_x_discrete(limit=c("jkdz","jbdz","PT","ZX") )+ - scale_color_manual(limits=c("jkdz","jbdz","PT","ZX"), values=c("#85B22E","#5F80B4","#E29827","#922927")) - ggsave(paste0(filename, ".pdf"),plot=a,width=8,height=8) -} - plot.boxplot( data1,data1$type,data1$value,data1$type,paste0("prot7_boxplot/",i, "_boxplot"),title=paste0(nm, pv)) + ge.plot.boxplot( data1,data1$type,data1$value,data1$type,paste0("prot7_boxplot/",nm, "_boxplot"),title=paste0(nm, pv)) } - - -which(grepl("choline",row.names(df2) )) - -``` - ``` -