rm(list=ls())
phe <- read.csv("phe(2).csv")
DAS28_ESR <- phe[,1:2]
fenceng <- read.table("exp.txt")
data <- readRDS("expr_data.rds")
table(rownames(data) %in% colnames(fenceng))
##DEGS
library(limma)
group_data <- phe[,c(1,3)]
data <- data.frame(t(data))
data$X <- rownames(data)
data <- merge(group_data,data,by="X")
data_sorted <- data[order(data[[2]], decreasing = TRUE), ]
rownames(data_sorted) <- data_sorted[,1]
group <- factor(data_sorted[,2],levels = c("RA","healthy control"))
data1 <- data.frame(t(data_sorted[,-c(1,2)]))
dimnames=list(rownames(data1),colnames(data1))
exprSet=matrix(as.numeric(as.matrix(data1)),nrow=nrow(data1),dimnames=dimnames)
exprSet_symbol1 <- aggregate(x = exprSet,
by = list(rownames(exprSet)),
FUN = mean)
rownames(exprSet_symbol1)=exprSet_symbol1[,1]
exprSet_symbol1=exprSet_symbol1[,2:ncol(exprSet_symbol1)]
dimnames=list(rownames(exprSet_symbol1),colnames(exprSet_symbol1))
exprSet=matrix(as.numeric(as.matrix(exprSet_symbol1)),nrow=nrow(exprSet_symbol1),dimnames=dimnames)
exprSet=exprSet[rowMeans(exprSet)>0,]
res.pca <- prcomp(t(exprSet), scale = TRUE)
library(factoextra)
fviz_pca_ind(res.pca,col.ind = group)
design <- model.matrix(~group)
colnames(design) <- levels(group)
design
fit <- lmFit(exprSet,design)
fit2 <- eBayes(fit)
allDiff=topTable(fit2,adjust='fdr',coef=2,number=Inf)
saveRDS(allDiff,"allDiff.rds")
#allDiff <- readRDS("allDiff.rds")
library(dplyr)
diffgene <- allDiff %>%
filter(adj.P.Val < 0.05) %>%
filter(abs(logFC) >0.5)
LDYDEGS <- read.csv("李东宇DEGS.csv",row.names = 1)
table(rownames(diffgene) %in% rownames(LDYDEGS))
JIAOJI_DEGS <- diffgene[rownames(diffgene) %in% rownames(LDYDEGS),]
saveRDS(diffgene,file = "diffgene.rds")
saveRDS(JIAOJI_DEGS,"JIAOJI_DEGS.rds")
saveRDS(exprSet,"exprSet.rds")
给我讲解以上的R代码