芯片分析流程-CSDN博客

本文链接：https://blog.csdn.net/m0_68600081/article/details/133955959

microarray

芯片分析

rm(list=ls())
setwd()

install.packages('AnnoProbe')
library(AnnoProbe)
#更新镜像库
devtools::install_git("https://gitee.com/jmzeng/GEOmirror")
#使用中国镜像下载GEO数据
GSE8146 <- AnnoProbe::geoChina(gse='GSE8146', mirror = 'tencent', destdir = '.')

这是Jimmy的包

library(GEOquery)
data=GSE8146

counts=(exprs(data[[1]]))[,1:10]#获取表达矩阵
dim(counts)
pData=pData(data[[1]])#获取样本信息，里面有很多，需要仔细观察下
sample=pData$geo_accession[1:10]
group=as.factor(rep(c('young','old'),each=5))
design_df=data.frame(sample=sample,
                                group=group)#创建分组信息

##注释##
#BiocManager::install('mgu74av2.db')
library(mgu74av2.db)#加载芯片注释包
columns(mgu74av2.db) #查看芯片注释包提供的注释信息
#showQCData("mgu74av2", datacache)
ids_GPL81=toTable(mgu74av2SYMBOL)

counts=as.data.frame(counts)
counts$probe_id=as.character(rownames(counts)) #将行名转化为probe_id
exprSet <- merge(ids_GPL81,counts , by = "probe_id") #合并探针注释和表达矩阵
table(is.na(exprSet))
head(exprSet)
exprSet <- exprSet[order(rowMeans(exprSet[, 3:12]), 
                                           decreasing = TRUE), ]
exprSet <- exprSet[!duplicated(exprSet$symbol), ]
rownames(exprSet) <- exprSet$symbol
exprSet <- exprSet[, -c(1, 2)]

#数据标准化、归一化处理
library(DESeq2)
p_before_norm=boxplot(exprSet)#具体美化的数据再自行百度了
exprSet=as.matrix(exprSet)
head(exprSet)
exprSet<-ceiling(exprSet) #这一步是将数据取整数,这里很奇怪的为什么count有小数点，这需要你去了解上游，现在一般用Hista2不用出现小数点的情况
head(exprSet)  
exprSet_rlog=rlog(exprSet)#用的是DEseq2包里面的函数来标准化
p_after_norm=boxplot(exprSet_rlog)
dev.off()#多按几次直到报错


#过滤
keep_feature <- rowSums(exprSet_rlog) >=30 #＜0的都过滤掉，这里的30自己去想一下
table(keep_feature)
exprSet_rlog=exprSet_rlog[keep_feature,]
p_after_norm2=boxplot(exprSet_rlog)
dev.off()#多按几次直到报错

#处理批次效应#批次效应自己去学下，这里没存在批次效应
library(limma)
#removeBatchEffect()
#batch <- c(rep("A",10),rep("B",20))
#exp2 <- removeBatchEffect(exp, batch)
#par(mfrow=c(1,2))  # 展示的图片为一行两列
#boxplot(as.data.frame(exp),main="Original")
#boxplot(as.data.frame(exp2),main="Batch corrected")

###查看一下数据##主要三大图，PCA、热图、聚类图
#第一步查看一下常见的基因，比如内参
exprSet_qc<-exprSet_rlog
exprSet_qc['GAPDH',]#这里没发现GAPDH，没事
boxplot(exprSet_qc['GAPDH',])

exprSet_qc['Actb',]
boxplot(exprSet_qc['Actb',])


##第二步用boxplot##
library(reshape2)

exprSet_qc_L=melt(exprSet_qc)
colnames(exprSet_qc_L) <- c('probe','sample','value')
head(exprSet_qc_L)
exprSet_qc_L$group=rep(group,each=nrow(exprSet_qc))
head(exprSet_qc_L)
### ggplot2 
library(ggplot2)
p=ggplot(exprSet_qc_L,
         aes(x=sample,y=value,fill=group))+geom_boxplot()
print(p)


##第三步用PCA## 

# BiocManager::install('ggfortify')
library(ggfortify)
exprSet_PCA=as.data.frame(t(exprSet_qc))
exprSet_PCA$group=design_df$group
autoplot(prcomp( exprSet_PCA[,1:(ncol(exprSet_PCA)-1)] ), data=exprSet_PCA,colour = 'group')

library("FactoMineR")#画主成分分析图需要加载这两个包
library("factoextra") 
dat.pca <- PCA(exprSet_PCA[,-ncol(exprSet_PCA)], graph = FALSE)#现在dat最后一列是group，需要重新赋值给一个dat.pca,这个矩阵是不含有分组信息的
fviz_pca_ind(dat.pca,
             geom.ind = "point", # show points only (nbut not "text")
             col.ind =group , # color by groups
             palette = c("#00AFBB", "#E7B800"),
             addEllipses = TRUE, # Concentration ellipses
             legend.title = "Groups"
)

#3D_PCA#
library(scatterplot3d)
library(dplyr)
data.pca <- prcomp(t(exprSet_qc),scale. = TRUE)
pca.result<-as.data.frame(data.pca$x)
pca.result$group<-group
scatterplot3d(pca.result[,1:3])
pca.result = pca.result %>% mutate(colour = case_when(
  pca.result$group == "young" ~ "#CC0000",
  pca.result$group == "old" ~ "#2f5688",
)) #定义分组并设置颜色

scatterplot3d(pca.result[,1:3],color=pca.result$colour,
              pch = 16,angle=30,
              cex.symbols = 2,
              box=T,type="p",
              main = "3D PCA Plot",
              lty.hide=2,lty.grid = 2)

legend("bottom",c("young","old"),fill=c('#CC0000',"#2f5688"),
       horiz = T)#添加图例,位置自己再设置一下
#直接在上面标注#
sd3 <- scatterplot3d(pca.result[,1:3],color=pca.result$colour,
                     pch = 16,angle=30,
                     cex.symbols = 2,
                     box=T,type="p",
                     main = "3D PCA Plot",
                     lty.hide=2,lty.grid = 2)
text(sd3$xyz.convert(pca.result[,c("PC1","PC2","PC3")]+0.05),
     labels = rownames(pca.result),
     cex=1,col="black")


#第四步用hclust##
#层次聚类图#
#首先针对样本做个系统聚类树
install.packages('WGCNA')
library(WGCNA)
datExpr_tree<-hclust(dist(t(exprSet_qc)), method = "ward.D")#这里使用的是欧氏距离，ward.D聚类
par(mar = c(0,5,2,0))
plot(datExpr_tree, main = "Sample clustering", sub="", xlab="", cex.lab = 2, 
     cex.axis = 1, cex.main = 1,cex.lab=1)

#针对前面构造的样品矩阵添加对应颜色
sample_colors <- numbers2colors(as.numeric(factor(design_df$group)), 
                                colors = c("blue","red"),signed = FALSE)

par(mar = c(1,4,3,1),cex=0.8)
plotDendroAndColors(datExpr_tree, sample_colors,
                    groupLabels = colnames(sample),
                    cex.dendroLabels = 0.8,
                    marAll = c(1, 4, 3, 1),
                    cex.rowText = 0.01,
                    main = "Sample dendrogram and trait ")

dev.off())


library(RColorBrewer)
#第五步用热图#
#列名与列名之间的关系
#colnames(exprSet_GSE12480)=design_df_GSE12480$group
sampleDists <- dist(t(exprSet_qc))#dist默认计算矩阵行与行的距离， 因此需要转置 了解一下这里使用什么方法来算距离，我记得都默认欧式距离，这个我也有笔记
sampleDistMatrix <- as.matrix(sampleDists)  
colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)  #选取热图的颜色
p0 <- pheatmap::pheatmap(sampleDistMatrix,
                         fontsize=7,
                         clustering_distance_rows=sampleDists,
                         clustering_distance_cols=sampleDists,
                         angle_col=45,
                         col=colors)
#ggsave(p0,filename = 'check_dist.pdf',width = 7.5,height =6)
dev.off()


#基因与列名之间的关系#
library(melbench)
library(cluster)
library(pheatmap)
pheatmap(exprSet_qc,scale = 'row',cluster_rows = T,clustering_method = 'average',
         show_rownames = F,cluster_cols=datExpr_tree)

library(limma)
####limma
data_DEG=exprSet_rlog

design <- model.matrix(~0+factor(design_df$group))
colnames(design)=levels(factor(design_df$group))
rownames(design)=colnames(data_DEG)
design

## 下面的 contrast.matrix 矩阵非常重要，制定了谁比谁这个规则
exp="old"
ctr="young"


contrast.matrix <- makeContrasts(contrasts=paste0(exp,'-',ctr),  #"exp/ctrl"
                                 levels = design)
contrast.matrix 
##step1
fit <- lmFit(data_DEG,design)
##step2
fit2 <- contrasts.fit(fit, contrast.matrix) ##这一步很重要，大家可以自行看看效果
fit2 <- eBayes(fit2)  ## default no trend !!!
##eBayes() with trend=TRUE
##step3
tempOutput = topTable(fit2, coef=1, n=Inf)
nrDEG = na.omit(tempOutput) 
write.csv(nrDEG,"limma.results.csv",quote = F)
head(nrDEG)

## volcano plot
DEG=nrDEG
#logFC_cutoff <- with(DEG,mean(abs( logFC)) + 2*sd(abs( logFC)) )
logFC_cutoff=1
pvalue_cutoff=0.05

DEG$change = as.factor(ifelse(DEG$P.Value <pvalue_cutoff  & abs(DEG$logFC) > logFC_cutoff,
                              ifelse(DEG$logFC > logFC_cutoff ,'UP','DOWN'),'NOT')
)
this_tile <- paste0('Cutoff for logFC is 0.5',#round(logFC_cutoff,3),
                    '\nThe number of up gene is ',nrow(DEG[DEG$change =='UP',]) ,
                    '\nThe number of down gene is ',nrow(DEG[DEG$change =='DOWN',])
)
this_tile
head(DEG)
ggplot(
  #设置数据
  DEG, 
  aes(x = logFC, 
      y = -log10(P.Value), 
      colour=change)) +
  geom_point(alpha=0.4, size=3.5) +
  scale_color_manual(values=c("#546de5", "#d2dae2","#ff4757"))+
  
  # 辅助线
  geom_vline(xintercept=c(-1,1),lty=4,col="black",lwd=0.8) +
  geom_hline(yintercept = -log10(pvalue_cutoff),lty=4,col="black",lwd=0.8) +
  
  # 坐标轴
  labs(x="log2(fold change)",
       y="-log10 (p-value)")+
  theme_bw()+
  
  # 图例
  theme(plot.title = element_text(hjust = 0.5), 
        legend.position="right", 
        legend.title = element_blank()
  )
  
dev.off()