admin管理员组文章数量:1437369
GSE75421芯片数据分析
想着去年学习的时候，也是学到芯片数据分析处理，然后说要去学习一下TCGA数据分析，然后就没有然后了
那作为承接去年的笔记，今年还是先跑一跑芯片数据练习，然后这次真的要啃下来TCGA部分内容了！要是再不啃下来，欢迎大家评论区捶我！
今年还是有进步的，至少把小洁老师提供的课后资料全部跑通学习了一遍！
芯片数据集——GSE75421简介
用到的是课上学员找的例子，然后因为是4分组，所以可以用来分别做两分组和多分组的练习，今天先整理一下2分组的分析数据，算是和去年的芯片数据分析流程初探，简称芯片数据再探！
首先通过GEO数据编号，确定是array芯片数据，物种是Mus musculus，所以在后续进行富集分析的时候要注意选择正确的物种
总共是有12个sample，选择normal和tumor组进行第一次两分组的练习
step1-数据下载和整理
因为是2016年的数据，所以直接使用geoChina下载即可
代码语言：javascript代码运行次数：0运行复制#下载数据
library(AnnoProbe)
eSet = geoChina("GSE75421")
#提取数据
eSet = eSet[[1]]
下载得到的数据
整理后
提取表达矩阵和临床信息
提取表达矩阵
代码语言：javascript代码运行次数：0运行复制#(1)提取表达矩阵exp
exp <- exprs(eSet)
dim(exp) #[1] 24479 12
range(exp)#看数据范围决定是否需要log，是否有负值，异常值，如有负值，结合箱线图进一步判断
#[1] 1.1694 13.8943
boxplot(exp,las = 2) #看是否有异常样本
提取分组信息
代码语言：javascript代码运行次数：0运行复制#(2)提取临床信息
pd <- pData(eSet)
if(T){
k1 = str_detect(pd$title,"normal");table(k1)
k2 = str_detect(pd$title,"tumor");table(k2)
pd = pd[k1|k2,]
}
#(3)让exp列名与pd的行名顺序完全一致
p = identical(rownames(pd),colnames(exp));p
if(!p) {
s = intersect(rownames(pd),colnames(exp))
exp = exp[,s]
pd = pd[s,]
}
获取芯片平台编号
代码语言：javascript代码运行次数：0运行复制gpl_number <- eSet@annotation;gpl_number
#[1] "GPL18802"
step2-整理分组信息和探针注释
分组信息整理
代码语言：javascript代码运行次数：0运行复制# 1.Group----
library(stringr)
k = str_detect(pd$title,"normal");table(k) #不在title就在pd的其他列
Group = ifelse(k,"normal","tumor")
# 需要把Group转换成因子，并设置参考水平，指定levels
Group = factor(Group,levels = c("normal","tumor"))
Group
#⭐检查自己得到的分组是否正确，可以肉眼确认一下即可
data.frame(pd$title,Group)
探针注释信息整理——下载并读取GPL网页的表格文件，按列取子集
代码语言：javascript代码运行次数：0运行复制#根据芯片平台获取并下载注释文件
library(tinyarray)
get_gpl_txt(gpl_number,download = T) #获取表格文件的下载链接
#读取文件并提取需要的行进行整理
a = read.delim("GPL18802.txt",skip = 8,comment.char = "!")
colnames(a)
ids = a[,1:2]
#如果有空行怎么去掉
k = ids$GeneSymbol!=""
ids = ids[k,]
colnames(ids)=c("probe_id","symbol")
注释文件
ids
step3-pca及heatmap图
PCA图_如果少于4个点就不会画出圈，是正常的。因为圈是置信区间，样本太少无法计算，不是必须的。
代码语言：javascript代码运行次数：0运行复制1.PCA 图----
dat=as.data.frame(t(exp))
library(FactoMineR)
library(factoextra)
dat.pca <- PCA(dat, graph = FALSE)
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"
)
2.top 1000 sd 热图
代码语言：javascript代码运行次数：0运行复制g = names(tail(sort(apply(exp,1,sd)),1000)) #day7-apply的思考题
n = exp[g,]
library(pheatmap)
annotation_col = data.frame(row.names = colnames(n),
Group = Group)
pheatmap(n,
show_colnames =F,
show_rownames = F,
annotation_col=annotation_col,
scale = "row", #按行标准化，只保留行内差别，不保留行间差别，会把数据范围缩放到大概-5~5之间
breaks = seq(-3,3,length.out = 100) #设置色带分布范围为-3~3之间，超出此范围的数字显示极限颜色
)
如果不需要聚类，可以调整pheatmap的参数
step4-差异分析
代码语言：javascript代码运行次数：0运行复制#差异分析
library(limma)
design = model.matrix(~Group)
fit = lmFit(exp,design)
fit = eBayes(fit)
deg = topTable(fit,coef = 2,number = Inf)
#为deg数据框添加几列
#1.加probe_id列，把行名变成一列
library(dplyr)
deg = mutate(deg,probe_id = rownames(deg))
#2.加上探针注释
ids = distinct(ids,symbol,.keep_all = T)
deg = inner_join(deg,ids,by="probe_id")
nrow(deg) #如果行数为0就是你找的探针注释是错的。
#3.加change列,标记上下调基因
logFC_t = 1
p_t = 0.05
k1 = (deg$P.Value < p_t)&(deg$logFC < -logFC_t)
k2 = (deg$P.Value < p_t)&(deg$logFC > logFC_t)
deg = mutate(deg,change = ifelse(k1,"down",ifelse(k2,"up","stable")))
table(deg$change)
#火山图
library(ggplot2)
ggplot(data = deg, aes(x = logFC, y = -log10(P.Value))) +
geom_point(alpha=0.4, size=3.5, aes(color=change)) +
scale_color_manual(values=c("blue", "grey","red"))+
geom_vline(xintercept=c(-logFC_t,logFC_t),lty=4,col="black",linewidth=0.8) +
geom_hline(yintercept = -log10(p_t),lty=4,col="black",linewidth=0.8) +
theme_bw()
差异基因热图
代码语言：javascript代码运行次数：0运行复制exp = exp[deg$probe_id,]
rownames(exp) = deg$symbol
diff_gene = deg$symbol[deg$change !="stable"]
n = exp[diff_gene,]
library(pheatmap)
annotation_col = data.frame(group = Group)
rownames(annotation_col) = colnames(n)
pheatmap(n,show_colnames =F,
show_rownames = F,
scale = "row",
#cluster_cols = F,
annotation_col=annotation_col,
breaks = seq(-3,3,length.out = 100)
)
step5-富集分析
加ENTREZID列，用于富集分析
物种注释查看其他物种.html#___OrgDb
这个数据集是小鼠的，所以用到的是org.Mm.eg.db
library(clusterProfiler)
library(org.Mm.eg.db)
s2e = bitr(deg$symbol,
fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = org.Mm.eg.db)
#⭐检查行数，减少太多说明不正常
nrow(deg)
deg = inner_join(deg,s2e,by=c("symbol"="SYMBOL"))
#多了一些行少了一些行都正常，SYMBOL与ENTREZID不是一对一的。
nrow(deg)
clusterProfiler富集分析
代码语言：javascript代码运行次数：0运行复制#加载需要的R包
library(clusterProfiler)
library(ggthemes)
library(org.Mm.eg.db)
library(dplyr)
library(ggplot2)
library(stringr)
library(enrichplot)
整理需要的数据
代码语言：javascript代码运行次数：0运行复制#(1)输入数据
gene_diff = deg$ENTREZID[deg$change != "stable"]
#(2)富集-小鼠 Mus musculus
ekk <- enrichKEGG(gene = gene_diff,organism = 'mmu')
#其他物种.html
ekk <- setReadable(ekk,OrgDb = org.Mm.eg.db,keyType = "ENTREZID")
#如果ekk是空的，这句就会报错，因为没富集到任何通路。
ego <- enrichGO(gene = gene_diff,OrgDb= org.Mm.eg.db,
ont = "ALL",readable = TRUE)
#setReadable和readable = TRUE都是把富集结果表格里的基因名称转为symbol
class(ekk)
dotplot或者barplot可视化
代码语言：javascript代码运行次数：0运行复制barplot(ego, split = "ONTOLOGY") +
facet_grid(ONTOLOGY ~ ., space = "free_y",scales = "free_y")
barplot(ekk)
dotplot(ego, split = "ONTOLOGY") +
facet_grid(ONTOLOGY ~ ., space = "free_y",scales = "free_y")+
dotplot(ekk)
晚些试试看跑一下多分组分析的流程，hhh这次可是记住了小洁老师说的，不运行的代码就注释掉，运行完之后要全选运行检查代码！
本文参与 腾讯云自媒体同步曝光计划，分享自微信公众号。原始发表：2025-04-26，如有侵权请联系 cloudcommunity@tencent 删除数据芯片异常数据分析表格本文标签：
版权声明：本文标题：GSE75421芯片数据分析 内容由网友自发贡献，该文观点仅代表作者本人， 转载请联系作者并注明出处：http://www.betaflare.com/biancheng/1747497664a2700151.html， 本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌抄袭侵权/违法违规的内容，一经查实，本站将立刻删除。
发表评论