其他
最简单的芯片挖掘也会出错(菜鸟团周一数据挖掘专栏第？期)
你好啊，我是Christine，这里是生信菜鸟团周一数据挖掘专栏，欢迎和我们一起学习！
我们周日的文献阅读专栏：菜鸟团一周文献推荐（No.35）已经到了35期，实际上我也记不清楚现在周一数据挖掘专栏是第几期，希望读者可以帮我汇总一下！
本周我们带大家看一个“诡异”的数据集，主要是为了传达Jimmy老师的一个观念：公共数据不全是对的，大家要破除迷信，敢于质疑。
数据：GSE128502
实验设计：鼻咽癌细胞系 + 敲低PNPD的细胞系，每组各3个生物学重复，共6个样本。hug133plus3基因表达芯片。
文章的分析流程：差异表达基因 + 通路网络分析
下载数据
### 下载GSE128502表达矩阵
f='./GSE128502_eSet.Rdata'
library(GEOquery)
if(!file.exists(f)){
gset <- getGEO('GSE128502', destdir=".",
AnnotGPL = F, ## 注释文件
getGPL = F) ## 平台文件
gset #查看数据
save(gset,file=f) ## 保存到本地
}
load('./GSE128502_eSet.Rdata') ## 载入数据
### 取表达矩阵
a=gset[[1]]
expr=exprs(a)
dim(expr)
expr[1:4,1:4]
### 取表型信息
pd=pData(a)
## 芯片注释
BiocManager::install("hgu133plus2.db")
library("hgu133plus2.db")
id2sb <- toTable(hgu133plus2SYMBOL)
dim(id2sb)
## 去掉表达矩阵中没有基因注释的探针
expr <- expr[rownames(expr) %in% id2sb$probe_id,]
dim(expr)
# 同一个基因只保留平均表达量最大的探针
identical(rownames(expr),id2sb$probe_id) #确保表达矩阵的行名与探针注释一致
id2sb$mean <- apply(expr, 1, mean) #计算平均表达量
id2sb <- id2sb[order(id2sb$symbol,id2sb$mean,decreasing = T),] #按symbol和平均表达量排序
id2sb <- id2sb[!duplicated(id2sb$symbol),] #去除重复
expr <- expr[id2sb$probe_id,]
rownames(expr) <- id2sb$symbol
save(expr,pd,file = "./DEG_input.Rdata")
差异分析
# 差异分析
library(limma)
pd$title
group_list <- c(rep("CK",3),rep("Treat",3))
#分组矩阵
design <- model.matrix(~0+factor(group_list))
colnames(design) <- levels(factor(group_list))
rownames(design) <- rownames(pd)
#比较矩阵
cont.matrix <- makeContrasts(contrasts="Treat-CK",levels = design)
#差异表达分析
fit <- lmFit(expr, design)
fit2 <-contrasts.fit(fit,cont.matrix)
fit2 <- eBayes(fit2)
#提取结果
deg <- topTable(fit2, coef="Treat-CK", n=Inf)
head(deg)
#下调最显著的10个基因
head(rownames(deg[order(deg$logFC),]),10)
#上调最显著的10个基因
head(rownames(deg[order(deg$logFC,decreasing = T),]),10)
虽然PDPN的表达量的确下调了，但是校正后的P值居然都等于1？说明结果很不可靠。
而此时上下调最显著的基因和作者的结果类似，分析流程应该是对的，那很可能是数据本身有问题了。
检查数据集
样本分组有没有问题？
对样本的表达进行聚类和PCA分析，发现确实有点问题，样本并没有按照对照组和处理组分开。
e <- expr
colnames(e) <- c("CK1","CK2","CK3","Treat1","Treat2","Treat3")
# 层次聚类
d <- dist(t(e))
clusters <- hclust(d)
plot(clusters)
#主成分分析
library(FactoMineR)
library(factoextra)
df_pca <- PCA(t(e),graph = F)
fviz_pca_ind(df_pca,
geom="point",
col.ind = group_list,
legend.title = "Groups")是不是样本名标错了？
处理组PDPN确实表达量降低了，所以应该没有标错。
会不会是某个样本芯片有问题？
检查各样本的全部基因表达量，非常一致，没问题。
library(RColorBrewer)
cols <- brewer.pal(6, "Set1")
boxplot(e,col=cols)
既然数据没有问题，那应该是实验的问题了，我觉得可能是：
细胞间自身差异过大，以至于超过处理干扰所带来的影响。
siRNA敲低不足，芯片数据看到PDPN的表达只降低了一半。
PDPN基因敲低可能确实对其它基因表达没什么影响。
最后，顺便画个图看看作者筛选到的差异基因，也确实不大对……
## 筛选出差异基因
up10_paper <- c("IFI27","IFI44L","IFI6","OAS1","TRIM22",
"IFITM1","OAS2","STAT1","CMPK2","OAS3")
down10_paper <- c("PDPN", "NALCN", "TMTC4", "ZFP1", "LINC00973",
"JPH3", "MMP12", "RGL3") # BOP1和SMIM2-AS1不在我注释的基因中，暂且忽略
df_plot <- as.data.frame(e[c(up10_paper,down10_paper),])
df_plot$group <- c(rep("up",10),rep("down",8))
## 画图
library(hrbrthemes)
library(GGally)
library(viridis)
ggparcoord(df_plot,
columns = 1:6, groupColumn = 7,
showPoints = TRUE,
scale = "globalminmax",
title = "Parallel Coordinate Plot for the top DEGs",
alphaLines = 0.3) +
scale_color_manual(values = c("blue","red")) +
theme_ipsum()+
theme(plot.title = element_text(size=10))
总结：
挖掘公共数据一定要记得检查数据！
其实在做差异分析之前，就应该检查分组是否符合预期，比如挑选方差最大的1000个基因进行PCA或者聚类分析。