#数据下载
rm(list = ls())
options(stringsAsFactors = F)
library(GEOquery)
gse = "GSE6281"
eSet <- getGEO(gse,
destdir = '.',#下载到的位置， "."表示当前目录
getGPL =F )#不要帮我下载好注释的文件
#(1)提取表达矩阵exp
exp <- exprs(eSet[[1]])
exp[1:4,1:4]
range(exp)
boxplot(exp)
#exp = log2(exp+1)先用range(exp)看看它的范围再决定要不要log
#再用boxplot(exp)看看箱式图中间的线整不整齐，如果不整齐的话
#加上一句代码exp=normalizeBetweenArrays(exp)矫正不整齐的样本,或者直接去除这几个样本
#可以用sva::ComBat()和lima::removeBatchEffect()矫正批次效应
#(2)提取临床信息
pd <- pData(eSet[[1]])
write.csv(pd,file = "linchuangxinxi.csv")
#(3)调整pd的行名顺序与exp列名完全一致
p = identical(rownames(pd),colnames(exp));p
if(!p) exp = exp[,match(rownames(pd),colnames(exp))]
#(4)提取芯片平台编号
gpl <- eSet[[1]]@annotation#@相当于列表里的列取子集
save(gse,pd,exp,gpl,file = "step1output.Rdata")
#1.group_list(实验分组)和ids(芯片注释)，每次都需要改
rm(list = ls())
load(file = "step1output.Rdata")
library(stringr)
#group_list----
#第一类，现成的某一列或在某列中包含。
pd$title
#第二类，自己生成
group_list=c(rep("control",times=3),rep("treat",times=3))
group_list
#第三类，ifelse
library(stringr)
group_list=ifelse(str_detect(pd$title,"control"),"control","treat")
#设置参考水平，对照在前，处理在后
group_list = factor(group_list,
levels = c("control","treat"))#levels后面第一个是参考水平
View(group_list)
group_list1 <- if(str_detect(pd$characteristics_ch1.3,"Petrolatum")){
"control"} else if(str_detect(pd$characteristics_ch1.3,"Nickel")){
"nickel"} else{"others"}
#2.ids-----------------
#方法1 BioconductorR包
gpl
#http://www.bio-info-trainee.com/1399.html
if(!require(hgu133plus2.db))BiocManager::install("hgu133plus2.db")
library(hgu133plus2.db)
ls("package:hgu133plus2.db")#表示看看这个R包里有哪些函数
ids <- toTable(hgu133plus2SYMBOL)
head(ids)
# 方法2 读取gpl页面的soft文件，按列取子集
# 方法3 官网下载
# 方法4 自主注释
save(exp,group_list,ids,file = "step2output.Rdata")
#这个脚本最重要
#基因去重复（一个基因对应三种探针）有三种办法：1、取平均值2、随机保留一个3、留下一个最大值
rm(list = ls())
load(file = "step1output.Rdata")
load(file = "step2output.Rdata")
#输入数据：exp和group_list
#Principal Component Analysis
#http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/112-pca-principal-component-analysis-essentials
dat=as.data.frame(t(exp))
library(FactoMineR)#画主成分分析图需要加载这两个包
library(factoextra)
# pca的统一操作走起
dat.pca <- PCA(dat, graph = FALSE)
pca_plot <- 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 = "Groups"
)
pca_plot
ggsave(plot = pca_plot,filename = paste0(gse,"PCA.png"))
save(pca_plot,file = "pca_plot.Rdata")
#热图
cg=names(tail(sort(apply(exp,1,sd)),1000))
n=exp[cg,]
#绘制热图
annotation_col=data.frame(group=group_list)
rownames(annotation_col)=colnames(n)
library(pheatmap)
pheatmap(n,
show_colnames =F,
show_rownames = F,
annotation_col=annotation_col,
scale = "row")
dev.off()
#co=cor(exp)；pheatmap(co)这两句代码可做相关性热图
rm(list = ls())
load(file = "step2output.Rdata")
library(dplyr)
exp1 <- as.data.frame(exp)
exp2 <- mutate(exp1,probe_id=rownames(exp))
exp3 <- inner_join(exp2,ids,by="probe_id")
exp4 = subset(exp3,symbol == "CCL13")
write.csv(exp4,file = "CCL13.csv")
#差异分析，用limma包来做
#需要表达矩阵和group_list，不需要改
library(limma)
design=model.matrix(~group_list)
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))
#或者 deg$probe_id=rownames(deg)
head(deg)
#2.加symbol列，火山图要用
deg <- inner_join(deg,ids,by="probe_id")#和merge功能类似但Innerjoin不改变顺序谁写在前面以谁为准
head(deg)
#按照symbol列去重复
deg <- deg[!duplicated(deg$symbol),]
#3.加change列,标记上下调基因
logFC_t=1
P.Value_t = 0.05
k1 = (deg$P.Value < P.Value_t)&(deg$logFC < -logFC_t)#可以用table(k1)看看具体有几个
k2 = (deg$P.Value < P.Value_t)&(deg$logFC > logFC_t)
change = ifelse(k1,"down",ifelse(k2,"up","stable"))#或者用case_when()
table(change)
deg <- mutate(deg,change)
#4.加ENTREZID列，用于富集分析（symbol转entrezid，然后inner_join）
library(ggplot2)
library(clusterProfiler)
library(org.Hs.eg.db)
s2e <- bitr(deg$symbol,
fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = org.Hs.eg.db)#人类
#其他物种http://bioconductor.org/packages/release/BiocViews.html#___OrgDb
deg <- inner_join(deg,s2e,by=c("symbol"="SYMBOL"))
save(group_list,deg,logFC_t,P.Value_t,file = "step4output.Rdata")
