2018-12-22日总结
2018-12-23 本文已影响34人
小梦游仙境
step1:output.Rdata
rm(list = ls()) ## 魔幻操作，一键清空~
options(stringsAsFactors = F)# 注意查看下载文件的大小，检查数据
f='GSE76275_eSet.Rdata'
library(GEOquery)# 这个包需要注意两个配置，一般来说自动化的配置是足够的。
#Setting options('download.file.method.GEOquery'='auto')
#Setting options('GEOquery.inmemory.gpl'=FALSE)
if(!file.exists(f)){
gset <- getGEO('GSE76275', destdir=".",
AnnotGPL = F, ## 注释文件
getGPL = F) ## 平台文件
save(gset,file=f) ## 保存到本地
}
load('GSE76275_eSet.Rdata') ## 载入数据
class(gset)
length(gset)
> class(gset[[1]])# 因为这个GEO数据集只有一个GPL平台，所以下载到的是一个含有一个元素的list
[1]"ExpressionSet"
> ?ExpressionSet #会在help弹出用Biobase这个包，a现在是一个很复杂的对象，取复杂的对象用a“a@”来取
a=gset[[1]]
dat=exprs(a) ## a现在是一个对象，每一个对象都有一个固定的函数，先中a这个对象对应的函数就是exprs。此为得到所有基因在所有样本的表达矩阵
dim(dat)# 看维度
# [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array
dat[1:4,1:4]##看前4行，前4列。。得到表达矩阵，那么就可以比较一个基因在TNBC和non-TNBC中是否有差异，进而也可以比较所有基因在TNBC和non-TNBC中是否有差异。
pd=pData(a) ## 挑选一些感兴趣的临床表型，pData可以拿到分组信息（目前基因名看不懂，需要把探针名转化为基因名）
trait=pd[,51:53]
head(trait)
trait$T=substring(trait[,2],2,2)
trait$N=substring(trait[,2],4,4)
trait$M=substring(trait[,2],6,6)
colnames(trait)=c('age','tmn','bmi','T','M','N')
head(trait)
save(trait,file='trait.Rdata')
group_list = ifelse(pd$characteristics_ch1.1=='triple-negative status: not TN',
'noTNBC','TNBC')
table(group_list)
save(dat,group_list,file = 'step1-output.Rdata')
## 下面的代码是更为详细的注释
######### 导入数据后查看数据情况，行、列名，数据维度、目标数据的大致情况，比如我们想看不同分组的各个基因的表达量，那么你就要查看分组在哪一行/列，怎么分组的，基因是什么情况，只要在掌握已有数据架构的情况下，才能随心所欲的调整表达矩阵（以下代码可以自行增减，达到查看数据情况的目的即可）
b = gset[[1]] ## 降级提取b
exprSet=exprs(b) ## 获取表达矩阵（如下图2），发现已是log处理后的数据。通常表达矩阵的原始数字从0但好几百万都不等，需要进行归一化处理。14488875 elements
pdata = pData(b) ## 使用函数?pData获取样本临床信息（如性别、年龄、肿瘤分期等等） 265 obs. of 69
colnames(pdata) ## 查看列名，即查看所包含的临床信息类型，找到triple negative status在第67列
length(colnames(pdata)) ## 查看pdata列数 69
pdata[,67] ## 查看第67列triple negative status的数据情况,发现按照"TN"，"not TN"排列
group_list = as.character(pdata[, 67]) ## 将67列改成字符
dim(exprSet) ## 查看矩阵的维度 54675行 265列
table(group_list) ## 查看67列信息
image
![]
step2-check
画PCA
google搜索：pca ggpubr
rm(list = ls()) ## 魔幻操作，一键清空~
options(stringsAsFactors = F)
load(file = 'step1-output.Rdata')# 每次都要检测数据
dat[1:4,1:4]
## 下面是画PCA的必须操作，需要看说明书。
dat=t(dat) #基因名和样本名转换过来
dat=as.data.frame(dat)
dat=cbind(dat,group_list)
library("FactoMineR")
library("factoextra")
# The variable group_list (index = 54676) is removed
# before PCA analysis
dat.pca <- PCA(dat[,-54676], graph = FALSE)
fviz_pca_ind(dat.pca,
geom.ind = "point", # show points only (nbut not "text")
col.ind = dat$group_list, # color by groups
palette = c("#00AFBB", "#E7B800"),
addEllipses = TRUE, # Concentration ellipses
legend.title = "Groups"
)
ggsave('all_samples_PCA.png')
画热图
rm(list = ls()) ## 魔幻操作，一键清空~
load(file = 'step1-output.Rdata')
dat[1:4,1:4]
dat[1,]
system.time(apply(dat,1,sd))
system.time(for(i in 1:nrow(dat)){
sd(dat[i,])
})
cg=names(tail(sort(apply(dat,1,sd)),1000))
library(pheatmap)
pheatmap(dat[cg,],show_colnames =F,show_rownames = F)
n=t(scale(t(dat[cg,])))
n[n>2]=2
n[n< -2]= -2
n[1:4,1:4]
pheatmap(n,show_colnames =F,show_rownames = F)#去掉下面的东西要根据包来调
ac=data.frame(g=group_list)
##接下来要标分组信息看说明书
rownames(ac)=colnames(n)##
## 可以看到TNBC具有一定的异质性，拿它来区分乳腺癌亚型指导临床治疗还是略显粗糙。
## 首先TNBC本身就可以继续分组，其次那些不是TNBC的病人跟T属于NBC病人的表达量是无法区分的。
pheatmap(n,show_colnames =F,show_rownames = F,annotation_col=ac,filename = 'heatmap_top1000_sd.png')
image
老大中间隐藏的代码过程
##选择变化最剧烈的基因 用sd标准差
注：现在基因名已经通过t（dat）而转换到第一列
sd（dat[1,]）#第一个基因在所有样本中的表达量，第一个基因在第一行
sd（dat[1,]）#第二个基因在所有样本中的表达量，第二个基因在第二行
apply（dat,1,sd）#对每一行 进行操作，“1”就是行
system.time(apply(dat,1,sd))
system.time(for(i in 1:nrow(dat)){sd(dat[i,])})
fivenum(apply(dat,1,sad))五分位数
tail(sort(apply(dat,1,sd)),1000)挑大的，挑1000个
cg=names(tail(sort(apply(dat,1,sd)),1000))
step3-DEG差异分析
rm(list = ls()) ## 魔幻操作，一键清空~
options(stringsAsFactors = F)
load(file = 'step1-output.Rdata')
dat[1:4,1:4]
table(group_list)
boxplot(dat[1,]~group_list)#1可修改为2...
t.test(dat[1,]~group_list)
library(ggpubr)
df=data.frame(gene=dat[1,],stage=group_list)
p < -ggboxplot(df,x="stage",y="gene",
color="stage",palette="jco",add="jitter")
也可以在外面改循环
rm(list = ls()) ## 魔幻操作，一键清空~
options(stringsAsFactors = F)
load(file = 'step1-output.Rdata')
dat[1:4,1:4]
table(group_list)
boxplot(dat[1,]~group_list)
t.test(dat[1,]~group_list)
bp=function(g){
library(ggpubr)
df=data.frame(gene=g,stage=group_list)
p <- ggboxplot(df, x = "stage", y = "gene",
color = "stage", palette = "jco",
add = "jitter")
# Add p-value
p + stat_compare_means()
}
bp(dat[1,])
bp(dat[2,])
bp(dat[3,])
bp(dat[4,])
image
image
limma包**
library(limma)
design=model.matrix(~factor( group_list ))
fit=lmFit(dat,design)
fit=eBayes(fit)
options(digits = 4)
#topTable(fit,coef=2,adjust='BH')
topTable(fit,coef=2,adjust='BH')
bp(dat['241662_x_at',])
bp(dat['207639_at',])
deg=topTable(fit,coef=2,adjust='BH',number = Inf) ###Inf:无穷
save(deg,file='deg.Rdata')
rm(list=ls())
load(file="deg.Rdata")
head(deg)
deg$probe_id=rownames(deg) ##增加deg中的rownames这一列到最后一列
library(hgu133plus2.db)
ids=toTable(hgu133plus2SYMBOL)#toTable:将基因名和探针对应起来
deg=merge(deg,ids,by='probe_id')
火山图
## for volcano
if(T){
nrDEG=deg
head(nrDEG)
attach(nrDEG)
plot(logFC,-log10(P.Value))
library(ggpubr)
df=nrDEG
df$v= -log10(P.Value)
ggscatter(df, x = "logFC", y = "v",size=0.5)
df$g=ifelse(df$P.Value>0.01,'stable',
ifelse( df$logFC >1.5,'up',
ifelse( df$logFC < -1.5,'down','stable') )
)
table(df$g)
df$name=rownames(df)
ggscatter(df, x = "logFC", y = "v",size=0.5,color = 'g')
ggscatter(df, x = "logFC", y = "v", color = "g",size = 0.5,
label = "symbol", repel = T,
#label.select = rownames(df)[df$g != 'stable'] ,
label.select = c('PROM1','AGR3','AGR2'),#从
palette = c("#00AFBB", "#E7B800", "#FC4E07") )
if(T){
nrDEG=deg
head(nrDEG)
attach(nrDEG)
plot(logFC,-log10(P.Value))
library(ggpubr)
df=nrDEG
df$v= -log10(P.Value)
ggscatter(df, x = "logFC", y = "v",size=0.5)
df$g=ifelse(df$P.Value>0.01,'stable',
ifelse( df$logFC >1.5,'up',
ifelse( df$logFC < -1.5,'down','stable') )
)
table(df$g)
df$name=rownames(df)
ggscatter(df, x = "logFC", y = "v",size=0.5,color = 'g')
ggscatter(df, x = "logFC", y = "v", color = "g",size = 0.5,
label = "symbol", repel = T,
#label.select = rownames(df)[df$g != 'stable'] ,
label.select = c('PROM1','AGR3','AGR2'),
palette = c("#00AFBB", "#E7B800", "#FC4E07") )
ggscatter(df, x = "AveExpr", y = "logFC",size = 0.2)
接下来,不知道是什么了，下面接着上面的代码
ggscatter(df, x = "AveExpr", y = "logFC",size = 0.2)
df$p_c = ifelse(df$P.Value<0.001,'p<0.001',
ifelse(df$P.Value<0.01,'0.001<p<0.01','p>0.01'))
table(df$p_c )
ggscatter(df,x = "AveExpr", y = "logFC", color = "p_c",size=0.2,
palette = c("green", "red", "black") )
}
image
![]
热图-第一步
image
if(T){
load(file = 'step1-output.Rdata')
# 每次都要检测数据
dat[1:4,1:4]
table(group_list)
x=deg$logFC
names(x)=deg$probe_id
cg=c(names(head(sort(x),100)),##取前100个探针
names(tail(sort(x),100)))##取后100个探针
library(pheatmap)
pheatmap(dat[cg,],show_colnames =F,show_rownames = F)
n=t(scale(t(dat[cg,])))
热图-第二步
n[n>2]=2
n[n< -2]= -2
n[1:4,1:4]
pheatmap(n,show_colnames =F,show_rownames = F)
热图-第三步
ac=data.frame(g=group_list)
rownames(ac)=colnames(n)
pheatmap(n,show_colnames =F,
show_rownames = F,
cluster_cols = F, ###不排序
annotation_col=ac)
}
ac=data.frame(g=group_list)
rownames(ac)=colnames(n)
pheatmap(n,show_colnames =F,
show_rownames = F,
##删掉了 cluster_cols = F,
annotation_col=ac)
}
step4:anno注释
rm(list=ls())
load(file="deg.Rdata")
head(deg)
deg$probe_id=rownames(deg) ##增加deg中的rownames这一列到最后一列
library(hgu133plus2.db)
ids=toTable(hgu133plus2SYMBOL)#toTable:将基因名和探针对应起来
deg=merge(deg,ids,by='probe_id')
rm(list = ls()) ## 魔幻操作，一键清空~
load(file = 'deg.Rdata')
head(deg)
## 不同的阈值，筛选到的差异基因数量就不一样，后面的超几何分布检验结果就大相径庭。
logFC_t=1
deg$g=ifelse(deg$P.Value>0.01,'stable',
ifelse( deg$logFC > logFC_t,'UP',
ifelse( deg$logFC < -logFC_t,'DOWN','stable') )
)
table(deg$g)
library(ggplot2)
library(clusterProfiler)
library(org.Hs.eg.db)
df <- bitr(unique(deg$symbol), fromType = "SYMBOL",
toType = c( "ENTREZID"),
OrgDb = org.Hs.eg.db)
head(df)
DEG=deg
head(DEG)
DEG=merge(DEG,df,by.y='SYMBOL',by.x='symbol')
head(DEG)
save(DEG,file = 'anno_DEG.Rdata')
gene_up= DEG[DEG$g == 'UP','ENTREZID']
gene_down=DEG[DEG$g == 'DOWN','ENTREZID']
gene_diff=c(gene_up,gene_down)
gene_all=as.character(DEG[ ,'ENTREZID'] )
data(geneList, package="DOSE")
head(geneList)
boxplot(geneList)
boxplot(DEG$logFC)
geneList=DEG$logFC
names(geneList)=DEG$ENTREZID
geneList=sort(geneList,decreasing = T)
## KEGG pathway analysis
### 做KEGG数据集超几何分布检验分析，重点在结果的可视化及生物学意义的理解。
if(T){
### over-representation test
kk.up <- enrichKEGG(gene = gene_up,
organism = 'hsa',
universe = gene_all,
pvalueCutoff = 0.9,
qvalueCutoff =0.9)
head(kk.up)[,1:6]
kk.down <- enrichKEGG(gene = gene_down,
organism = 'hsa',
universe = gene_all,
pvalueCutoff = 0.9,
qvalueCutoff =0.9)
head(kk.down)[,1:6]
kk.diff <- enrichKEGG(gene = gene_diff,
organism = 'hsa',
pvalueCutoff = 0.05)
head(kk.diff)[,1:6]
kegg_diff_dt <- as.data.frame(kk.diff)
kegg_down_dt <- as.data.frame(kk.down)
kegg_up_dt <- as.data.frame(kk.up)
down_kegg<-kegg_down_dt[kegg_down_dt$pvalue<0.05,];down_kegg$group=-1
up_kegg<-kegg_up_dt[kegg_up_dt$pvalue<0.05,];up_kegg$group=1
source('functions.R')
g_kegg=kegg_plot(up_kegg,down_kegg)
print(g_kegg)
ggsave(g_kegg,filename = 'kegg_up_down.png')
step5:WGCNA
genecard搜索
