芯片数据质控
1.前后对比箱线图
## --------------------------------------------数据检验和聚类------------------------------
exprSet_L <- melt(exprSet)
colnames(exprSet_L) <- c("sample", "value")
exprSet_L$group <- rep(pdata$group, each = nrow(exprSet))
exprSet_L[1:10, ]
## 箱线图
p21 <- ggplot(exprSet_L, aes(x = sample, y = value, fill = group)) +
geom_boxplot() +
theme(
axis.title.x = element_text(face = "italic"),
axis.text.x = element_text(angle = 45, vjust = 0.5),
legend.position = "top"
)
p21 # 数据不太规整
# #ggsave(p21,filename = ".\\plots\\bar_before.pdf",width = 18,height = 9,units = "cm")
# ##标准化后再图
# exprSet_1 <- normalizeBetweenArrays(exprSet) %>% data.frame()
# exprSet_L <- melt(exprSet_1)
# colnames(exprSet_L) = c('sample','value')
# exprSet_L$group = rep(group_info$group,each = nrow(exprSet))
# p22 <- ggplot(exprSet_L,aes(x = sample, y = value,fill = group)) +geom_boxplot()+
# theme(axis.title.x = element_text(face = 'italic'),
# axis.text.x = element_text(angle = 45 , vjust = 0.5),
# legend.position = 0)
# p22 #表准化后规整多了，但知道效果咋样
#
# plot_qc <- plot_grid(p21,p22,labels = "AUTO",nrow = 2)
#
# ggsave(plot_qc,filename = ".\\plots\\plot_qc.pdf",width = 18,height = 18,units = "cm")
# 把标准化后的数据赋值给原始数据，暂时不做
# exprSet <- exprSet_1
2.提琴图
## 小提琴图
# p1 <- ggplot(exprSet_L,aes(x = sample, y = value,fill = group)) +geom_violin()
# p1
3.密度图
## 密度图
# p3 <- ggplot(exprSet_L,aes(value,col=group)) +geom_density() + facet_wrap(~sample,nrow = 2)
# print(p3)
4.聚类图1
## hclust聚类图
colnames(exprSet) <- paste(pdata$title1, sep = "") ## 样本名称改变
# hc <- hclust(dist(t(exprSet))) #hclust需要用plot画图，此例不用
dist.r <- dist(t(exprSet))
res <- hcut(dist.r, k = 2, stand = TRUE)
# Visualize
cluster_unbatch <- fviz_dend(res,
# 加边框
rect = TRUE,
# 边框颜色
rect_border = "cluster",
# 边框线条类型
rect_lty = 2,
# 边框线条粗细
lwd = 1.2,
# 边框填充
rect_fill = T,
# 字体大小
cex = 0.8,
# 字体颜色
color_labels_by_k = T,
# 平行放置
horiz = T,
k_colors = "lancet", #指定色板
labels_track_height = 7,
main = element_blank()
)+theme(text = element_text(size = 8))
cluster_unbatch
# 数据不可用,需要校正批次效应
# 把样本名改回来
colnames(exprSet) <- pdata$geo_accession
5.PCA图1
# 画主成分分析图需要加载这两个包
library(FactoMineR)
library(factoextra)
# pca
dat.pca <- PCA(as.data.frame(t(exprSet)), graph = FALSE)
pca_plot1 <- fviz_pca_ind(dat.pca,
geom.ind = "point", # show points only (nbut not "text")
col.ind = group_list, # color by groups
# palette = c("#00AFBB", "#E7B800"), #修改颜色
addEllipses = TRUE, # Concentration ellipses
# legend.title = "GROUP"
) + theme(
legend.position = c(0.80, 1),
text = element_text(size = 8),
legend.key.size = unit(3, "mm"),
legend.title = element_blank()
)
pca_plot1
6.处理批次效应
library(sva)
model <- model.matrix(~ as.factor(pdata$group)) # "group"为分组所在列，pdata是表型文件(元数据）
### 重点来了！！！
### 重点来了！！！
### 重点来了！！！
# 方法一(ComBat)
combat_edata <- ComBat(dat = exprSet, batch = pdata$batch, mod = model) # exprSet是表达矩阵，batch是批次信息列，在表型文件里，下同
# 使用sva包中的combat函数!
# 依次是什么数据、批次来自于哪、想保留什么信息
7.再次查看聚类图
# 再次聚类查看
colnames(combat_edata) <- paste(pdata$title1, sep = "") ## 样本名称再次改变
hc <- hclust(dist(t(combat_edata)))
dist.r <- dist(t(combat_edata))
res <- hcut(dist.r, k = 2, stand = TRUE)
# Visualize
cluster_batched <- fviz_dend(res,
# 加边框
rect = TRUE,
# 边框颜色
rect_border = "cluster",
# 边框线条类型
rect_lty = 2,
# 边框线条粗细
lwd = 1.2,
# 边框填充
rect_fill = T,
# 字体大小
cex = 0.8,
# 字体颜色
color_labels_by_k = T,
# 平行放置
horiz = T,
k_colors = "lancet", #指定色板
labels_track_height = 7,
main = element_blank()
)+theme(text = element_text(size = 8))
cluster_batched
plot_batch <- plot_grid(cluster_unbatch,cluster_batched,labels = "AUTO",label_size = 10)
ggsave(plot_batch,filename = ".\\plots\\cluster_batched.pdf",width = 180,height = 89,units = "mm")
8.再次作图PCA
# pca
dat.pca <- PCA(as.data.frame(t(combat_edata)), graph = FALSE)
pca_plot2 <- fviz_pca_ind(dat.pca,
geom.ind = "point", # show points only (nbut not "text")
col.ind = group_list, # color by groups
# palette = c("#00AFBB", "#E7B800"), #修改颜色
addEllipses = TRUE, # Concentration ellipses
# legend.title = "GROUP"
) + theme(
legend.position = c(0.80, 1),
text = element_text(size = 8),
legend.key.size = unit(3, "mm"),
legend.title = element_blank()
)
pca_plot2
plot_batch <- plot_grid(pca_plot1,pca_plot2,labels = "AUTO",label_size = 10)
ggsave(plot_batch,filename = ".\\plots\\pca_batched.pdf",width = 180,height = 89,units = "mm")
9.后续处理
# 把校正批次效应后的文件赋给表达矩阵
exprSet <- combat_edata %>% data.frame()
# 再次改回样本名
colnames(exprSet) <- pdata$geo_accession
# 这样写就行，不用加其他参数，供CIBERSORT使用,生成TXT后需要到文件里把第一行向后tab一下，否则报错。
write.table(exprSet,paste0(".\\data\\",gse_id,"_exprset_to_cibersort.txt"),row.names = T,quote = F,sep = "\t")
质控处理完毕，准备差异表达分析。
