options("repos"="https://mirrors.ustc.edu.cn/CRAN/")
if(!require("BiocManager")) install.packages("BiocManager",update = F,ask = F)
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")

cran_packages <- c('tidyr',
                   'tibble',
                   'dplyr',
                   'stringr',
                   'ggplot2',
                   'ggpubr',
                   'factoextra',
                   'FactoMineR',
                   'devtools',
                   'cowplot',
                   'patchwork',
                   'basetheme',
                   'ggthemes',
                   'paletteer',
                   'AnnoProbe',
                   'tinyarray') 
Biocductor_packages <- c('GEOquery',
                         'hgu133plus2.db',
                         'ggnewscale',
                         "limma",
                         "impute",
                         "GSEABase",
                         "GSVA",
                         "clusterProfiler",
                         "org.Hs.eg.db",
                         "preprocessCore",
                         "enrichplot",
                         "ggplotify")

for (pkg in cran_packages){
  if (! require(pkg,character.only=T) ) {
    install.packages(pkg,ask = F,update = F)
    require(pkg,character.only=T) 
  }
}


for (pkg in Biocductor_packages){
  if (! require(pkg,character.only=T) ) {
    BiocManager::install(pkg,ask = F,update = F)
    require(pkg,character.only=T) 
  }
}


for (pkg in c(Biocductor_packages,cran_packages)){
  require(pkg,character.only=T) 
}
library(GEOquery)
gse_number = "GSE43292"
gset = getGEO('GSE43292', destdir='.',getGPL = F)
gset=gset[[1]]
pdata=pData(gset)
library("tidyverse")
group_list=ifelse(str_detect(pdata$title,"plaque"),
                  "Atherosclerosis",
                  "Control")
group_list=factor(group_list,
                  levels = c("Control","Atherosclerosis"))
exprSet=exprs(gset)
boxplot(exprSet,outline=FALSE, notch=T,col=group_list, las=2)
library(limma) 
exprSet=normalizeBetweenArrays(exprSet)
x=boxplot(exprSet,outline=FALSE, notch=T,col=group_list, las=2)

library(tinyarray)
if(!require(hugene10sttranscriptcluster.db))BiocManager::install("hugene10sttranscriptcluster.db")
library(hugene10sttranscriptcluster.db)
library(hugene10sttranscriptcluster.db)
ls("package:hugene10sttranscriptcluster.db")
ids <- toTable(hugene10sttranscriptclusterSYMBOL)
head(ids)

library(limma)
design=model.matrix(~group_list)
fit=lmFit(exprSet,design)
fit=eBayes(fit)
deg=topTable(fit,coef=2,number = Inf)
library(dplyr)
deg <- mutate(deg,probe_id=rownames(deg))
head(deg)
ids = ids[!duplicated(ids$symbol),]
deg <- inner_join(deg,ids,by="probe_id")
head(deg)
nrow(deg)
logFC_t=0.5
P.Value_t = 0.05
k1 = (deg$adj.P.Val < P.Value_t)&(deg$logFC < -logFC_t)
k2 = (deg$adj.P.Val< P.Value_t)&(deg$logFC > logFC_t)
deg <- mutate(deg,change = ifelse(k1,"down",ifelse(k2,"up","stable")))
head(deg)
table(deg$change)
library(clusterProfiler)
library(org.Hs.eg.db)
s2e <- bitr(deg$symbol, 
            fromType = "SYMBOL",
            toType = "ENTREZID",
            OrgDb = org.Hs.eg.db)
dim(deg)
deg <- inner_join(deg,s2e,by=c("symbol"="SYMBOL"))
dim(deg)
write.table(deg, file = "20221019-DEG-adj--p.Value.csv", sep = ",", col.names = NA)
length(unique(deg$symbol))
library(dplyr)
library(ggplot2)
dat  = deg[!duplicated(deg$symbol),]
deg<-read.csv("20221019-DEG-adj--p.Value.csv")
length(unique(deg$symbol))

library(dplyr)
library(ggplot2)
dat  = deg[!duplicated(deg$symbol),]
p <- ggplot(data = dat, 
            aes(x = logFC, 
                y = -log10(adj.P.Val))) +
  geom_point(alpha=10, size=1.2, 
             aes(color=change)) +
  ylab("-log10(adj.P.Val)")+
  scale_color_manual(values=c("blue", "grey","red"))+
  geom_vline(xintercept=c(-logFC_t,logFC_t),lty=4,col="black",lwd=0.8) +
  geom_hline(yintercept = -log10(P.Value_t),lty=4,col="black",lwd=0.8) +
  theme_bw()
p
head(dat)
exp = exprSet[dat$probe_id,]
rownames(exp) = dat$symbol
if(T){
  cg = dat$symbol[dat$change !="stable"]
  length(cg)
}else{
  x=dat$logFC[dat$change !="stable"]
  names(x)=dat$symbol[dat$change !="stable"]
  cg=names(c(head(sort(x),10),tail(sort(x),10)))
  length(cg)
}
n=exp[cg,]
dim(n)
library(pheatmap)
annotation_col=data.frame(group=group_list)
annotation_col = data.frame(
  Group = factor(rep(c("Control", "Atherosclerosis"), 32)))
rownames(annotation_col)=colnames(n)
write.table(n, file = "nnew.csv", sep = ",", col.names = NA)
nnew <- read.table(file = "C:\\Users\\22788\\Desktop\\GSE43292 NEW\\1\\shishi.txt", sep = "\t", header=T)
rownames(nnew) = nnew[,1]
nn = nnew[,-1]
head(nn)

heatmap_plot <- pheatmap(nn,show_colnames =T,show_rownames= F,cluster_cols=F,cluster_rows=T,
                         scale = "row",clustering_method = "complete",fontsize_row=1,
                         annotation_col=annotation_col,
                         breaks = seq(-2.7,2.5,length.out =110))
heatmap_plot

##
library(GEOquery)
gse_number = "GSE28829"
gset = getGEO('GSE28829', destdir='.',getGPL = F)
gset=gset[[1]]
pdata=pData(gset)
library("tidyverse")
group_list=ifelse(str_detect(pdata$title,"advanced"),
                  "early plaque",
                  "advanced plaque")
group_list=factor(group_list,
                  levels = c("early plaque","advanced plaque"))
exprSet=exprs(gset)
boxplot(exprSet,outline=FALSE, notch=T,col=group_list, las=2)
library(limma)
exprSet=normalizeBetweenArrays(exprSet)
x=boxplot(exprSet,outline=FALSE, notch=T,col=group_list, las=2)
library(tinyarray)
if(!require(hgu133plus2.db))BiocManager::install("hgu133plus2.db")
library(hgu133plus2.db)
library(hgu133plus2.db)
ls("package:hgu133plus2.db")
ids <- toTable(hgu133plus2SYMBOL)
head(ids)
library(limma)
design=model.matrix(~group_list)
fit=lmFit(exprSet,design)
fit=eBayes(fit)
deg=topTable(fit,coef=2,number = Inf)
library(dplyr)
deg <- mutate(deg,probe_id=rownames(deg))
head(deg)
ids = ids[!duplicated(ids$symbol),]
deg <- inner_join(deg,ids,by="probe_id")
head(deg)
nrow(deg)
logFC_t=0.5
P.Value_t = 0.05
k1 = (deg$adj.P.Val < P.Value_t)&(deg$logFC < -logFC_t)
k2 = (deg$adj.P.Val< P.Value_t)&(deg$logFC > logFC_t)
deg <- mutate(deg,change = ifelse(k1,"down",ifelse(k2,"up","stable")))
head(deg)
table(deg$change)
library(clusterProfiler)
library(org.Hs.eg.db)
s2e <- bitr(deg$symbol, 
            fromType = "SYMBOL",
            toType = "ENTREZID",
            OrgDb = org.Hs.eg.db)
dim(deg)
deg <- inner_join(deg,s2e,by=c("symbol"="SYMBOL"))
dim(deg)
length(unique(deg$symbol))
library(dplyr)
library(ggplot2)
dat  = deg[!duplicated(deg$symbol),]
deg<-read.csv("1019DEG-adj--p.Value.csv")
length(unique(deg$symbol))

library(dplyr)
library(ggplot2)
dat  = deg[!duplicated(deg$symbol),]
p <- ggplot(data = dat, 
            aes(x = logFC, 
                y = -log10(adj.P.Val))) +
  geom_point(alpha=10, size=1.2, 
             aes(color=change)) +
  ylab("-log10(adj.P.Val)")+
  scale_color_manual(values=c("blue", "grey","red"))+
  geom_vline(xintercept=c(-logFC_t,logFC_t),lty=4,col="black",lwd=0.8) +
  geom_hline(yintercept = -log10(P.Value_t),lty=4,col="black",lwd=0.8) +
  theme_bw()
p
head(dat)
exp = exprSet[dat$probe_id,]
rownames(exp) = dat$symbol
if(T){
  cg = dat$symbol[dat$change !="stable"]
  length(cg)
}else{
  x=dat$logFC[dat$change !="stable"]
  names(x)=dat$symbol[dat$change !="stable"]
  cg=names(c(head(sort(x),10),tail(sort(x),10)))
  length(cg)
}
n=exp[cg,]
dim(n)
library(pheatmap)
annotation_col=data.frame(group=group_list)
annotation_col = data.frame(
  Group = factor(rep(c("early plaque","advanced plaque"),times = c(13,16))))
rownames(annotation_col)=colnames(n)
write.table(n, file = "nnew.csv", sep = ",", col.names = NA)
nnew <- read.table(file = "C:\\Users\\22788\\Desktop\\1-GSE 28829 NEW\\1\\nnnew.txt", sep = "\t", header=T)
rownames(nnew) = nnew[,1]
nn = nnew[,-1]
head(nn)

heatmap_plot <- pheatmap(nn,show_colnames =T,show_rownames= F,cluster_cols=F,cluster_rows=T,
                         scale = "row",clustering_method = "complete",fontsize_row=1,
                         annotation_col=annotation_col,
                         breaks = seq(-2.7,2.5,length.out =110))
heatmap_plot

library(GEOquery)
gse_number = "GSE54666"
gset = getGEO('GSE54666', destdir='.',getGPL = F)
gset=gset[[1]]
pdata=pData(gset)
library("tidyverse")
View(pdata)
group_list=ifelse(str_detect(pdata$title,"oxLDL"),
                  "ox-LDL macrophage",
                  "Control")
group_list=factor(group_list,
                  levels = c("Control","ox-LDL macrophage"))
exprSet=exprs(gset)
View(exprSet)
boxplot(exprSet,outline=FALSE, notch=T,col=group_list, las=2)
library(limma)
exprSet=normalizeBetweenArrays(exprSet)
x=boxplot(exprSet,outline=FALSE, notch=T,col=group_list, las=2)
View(pdata)
View(exprSet)
exprSet=exprs(gset)
boxplot(exprSet,outline=FALSE, notch=T,col=group_list, las=2)
library(limma)
exprSet=normalizeBetweenArrays(exprSet)
x=boxplot(exprSet,outline=FALSE, notch=T,col=group_list, las=2)
library(tinyarray)
if(!require(illuminaHumanv4.db))BiocManager::install("illuminaHumanv4.db")
library(illuminaHumanv4.db)
ls("package:illuminaHumanv4.db")
ids <- toTable(illuminaHumanv4SYMBOL)
head(ids)
library(limma)
design=model.matrix(~group_list)
fit=lmFit(exprSet,design)
library(limma)
design=model.matrix(~group_list)
fit=lmFit(exprSet,design)
View(design)
View(pdata)
group_list=ifelse(str_detect(pdata$title,"oxLDL"),
                  "ox-LDL macrophage",
                  "Control")
group_list=factor(group_list,
                  levels = c("Control","ox-LDL macrophage"))
exprSet=exprs(gset)
boxplot(exprSet,outline=FALSE, notch=T,col=group_list, las=2)
library(limma)
exprSet=normalizeBetweenArrays(exprSet)
x=boxplot(exprSet,outline=FALSE, notch=T,col=group_list, las=2)
library(tinyarray)
if(!require(illuminaHumanv4.db))BiocManager::install("illuminaHumanv4.db")
library(illuminaHumanv4.db)
ls("package:illuminaHumanv4.db")
ids <- toTable(illuminaHumanv4SYMBOL)
head(ids)
library(limma)
design=model.matrix(~group_list)
fit=lmFit(exprSet,design)
library(limma)
design=model.matrix(~group_list)
fit=lmFit(exprSet,design)
fit=eBayes(fit)
deg=topTable(fit,coef=2,number = Inf)
library(dplyr)
deg <- mutate(deg,probe_id=rownames(deg))
head(deg)
ids = ids[!duplicated(ids$symbol),]
deg <- inner_join(deg,ids,by="probe_id")
head(deg)
nrow(deg)
logFC_t=0.5
P.Value_t = 0.05
k1 = (deg$adj.P.Val < P.Value_t)&(deg$logFC < -logFC_t)
k2 = (deg$adj.P.Val< P.Value_t)&(deg$logFC > logFC_t)
deg <- mutate(deg,change = ifelse(k1,"down",ifelse(k2,"up","stable")))
head(deg)
table(deg$change)
library(clusterProfiler)
library(org.Hs.eg.db)
s2e <- bitr(deg$symbol,
            fromType = "SYMBOL",
            toType = "ENTREZID",
            OrgDb = org.Hs.eg.db)
dim(deg)
deg <- inner_join(deg,s2e,by=c("symbol"="SYMBOL"))
dim(deg)
data1 <- read.delim("D:/jia/pca/data1.txt", row.names=1)
tdata<-t(data1)
tdist<-vegdist(tdata，na.rm = T)
group<-c("Ox-LDL macrophage","Control","Ox-LDL macrophage","Control","Ox-LDL macrophage","Control","Ox-LDL macrophage","Control","Ox-LDL macrophage","Ox-LDL macrophage","Control","Control")
anosin<-anosim(tdist,group)
summary(anosin)
deg<-read.csv("1019DEG-adj--p.Value.csv")
length(unique(deg$symbol))

library(dplyr)
library(ggplot2)
dat  = deg[!duplicated(deg$symbol),]
p <- ggplot(data = dat, 
            aes(x = logFC, 
                y = -log10(adj.P.Val))) +
  geom_point(alpha=10, size=2, 
             aes(color=change)) +
  ylab("-log10(adj.P.Val)")+
  scale_color_manual(values=c("blue", "grey","red"))+
  geom_vline(xintercept=c(-logFC_t,logFC_t),lty=4,col="black",lwd=0.8) +
  geom_hline(yintercept = -log10(P.Value_t),lty=4,col="black",lwd=0.8) +
  theme_bw()
p
head(dat)
exp = exprSet[dat$probe_id,]
rownames(exp) = dat$symbol
if(T){
  cg = dat$symbol[dat$change !="stable"]
  length(cg)
}else{
  x=dat$logFC[dat$change !="stable"] 
  names(x)=dat$symbol[dat$change !="stable"] 
  cg=names(c(head(sort(x),10),tail(sort(x),10)))
  length(cg)
}
n=exp[cg,]
dim(n)
library(pheatmap)
annotation_col=data.frame(group=group_list)

annotation_col = data.frame(
  Group = factor(rep(c("Control", "ox-LDL macrophage"), 6)))
rownames(annotation_col)=colnames(n) 
write.table(n, file = "nnew.csv", sep = ",", col.names = NA)
nnew <- read.table(file = "C:\\Users\\22788\\Desktop\\GSE546666  oxldl  micro\\oxldl\\2Nnew.txt", sep = "\t", header=T)
rownames(nnew) = nnew[,1]
nn = nnew[,-1]
head(nn)
heatmap_plot <- pheatmap(nn,show_colnames =T,show_rownames= F,cluster_cols=F,cluster_rows=T,
                         scale = "row",clustering_method = "complete",fontsize_row=1,
                         annotation_col=annotation_col,
                         breaks = seq(-3,3,length.out =110))
heatmap_plot