莎莎写于2020.5.7✍今天上午10点钟准时就醒了，然后和室友磨叽半天去实验室打扫卫生，然后签了字就偷偷溜了😄然后回到中校区买了超级甜的哈密瓜🍈还有火龙果，沃柑，一共是30￥ 关于今天要学习的内容就是：昨天晚上做PPT的时候，我想知道一些母基因的Regulation情况，但是现成的表格信息有很多NA，我就想着要不自己分析。然后我就找到一个示例数据集GSE14520，它对应的两个平台是GPL571和GPL3921

#数据下载
rm(list = ls())
options(stringsAsFactors = F)
# 构建结果储存文件夹
if(!dir.exists("data")) dir.create("data")
if(!dir.exists("figures")) dir.create("figures")
suppressMessages(library(GEOquery))
gse = "GSE14520"
eSet <- getGEO(gse, 
               destdir = '.', 
               getGPL = F)
# class(eSet)  #可以看一下是什么类型的数据：list
# length(eSet)  #看这个list中有多少个个元素：2
# class(eSet[[1]])  #提取列表中的子集

image.png

image.png

可以看到上面有2个elements，有两个元素

#(1)提取表达矩阵exp
exp <- exprs(eSet[[1]])   #GPL3921
exp1 <- exprs(eSet[[2]])   #GPL571
exp[1:4,1:4]
#dim(exp)  #22268行  445列
#exp = log2(exp+1)
#(2)提取临床信息
pd <- pData(eSet[[1]])
# class(pd)  #是一个data.frame
# dim(pd)   #6行 34列
#(3)调整pd的行名顺序与exp列名完全一致：用identical函数
p = identical(rownames(pd),colnames(exp));p  #输出结果是TRUE就不需要调整
if(!p) exp = exp[,match(rownames(pd),colnames(exp))]  #用match函数根据pd的行名来调整exp的列名

image.png

#(4)提取芯片平台编号
gpl <- eSet[[1]]@annotation
gpl1 <- eSet[[2]]@annotation

image.png

可以看到通过提取列表中的子集就可以提取自己想要的平台对应的数据，相当于两个平台当作两个数据集来做就可以了:

exp <- exprs(eSet[[1]])   #GPL3921
exp1 <- exprs(eSet[[2]])   #GPL571

我需要的就是GPL3921,下面开始进行平台注释

平台注释有很多种方法：

#2.ids-----------------
#方法1 BioconductorR包
gpl 
#http://www.bio-info-trainee.com/1399.html
if(!require(hugene10sttranscriptcluster.db))BiocManager::install("hugene10sttranscriptcluster.db")
library(hugene10sttranscriptcluster.db)
ls("package:hugene10sttranscriptcluster.db")
ids <- toTable(hugene10sttranscriptclusterSYMBOL)
head(ids)
# 方法2 读取gpl页面的soft文件，按列取子集
# 方法3 官网下载
# 方法4 自主注释 

我一开始尝试了从GEO网站下载，但是最后压缩文件是不完整的

image.png

正是我需要的

#2.ids-----------------
#如果镜像一直出现问题导致下载包打不开，可以试试下面这个代码,然后再下载
options(BioC_mirror="http://mirrors.cloud.tencent.com/bioconductor")  #这个是Bioconductor的包下载之前设置的镜像
options("repos" = c(CRAN="https://mirrors.aliyun.com/CRAN/")) #这个是CRAN上面的包下载之前要设置的镜像
#方法1 BioconductorR包
gpl 
#http://www.bio-info-trainee.com/1399.html
if(!require(hthgu133a.db))BiocManager::install("hthgu133a.db")
library(hthgu133a.db)
ls("package:hthgu133a.db")
ids <- toTable(hthgu133aSYMBOL)
head(ids)

内心想法：这个 包下的可真不容易啊，镜像也有问题，然后又设置了镜像，最后才下载好!😄

image.png

image.png

SYMBOL是我们需要的注释信息ids

##下面的代码就可以提取我们要的ids
ids <- toTable(hthgu133aSYMBOL)
head(ids)

image.png

然后保存这些变量，继续往下走：开始画PCA和热图

rm(list = ls())  
load(file = "data/step1output.Rdata")
load(file = "data/step2output.Rdata")
#输入数据：exp和group_list
#Principal Component Analysis:PCA详细原理
#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))  #先转置，然后再as.data.frame
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"), #这里可以修改颜色，网页搜索RGB选择子集喜欢的颜色替换即可
                         addEllipses = TRUE, # Concentration ellipses
                         legend.title = "Groups"
)
pca_plot
ggsave(plot = pca_plot,filename = paste0(gse,"PCA.png"))
save(pca_plot,file = "data/pca_plot.Rdata")

PCA图

#热图 
cg=names(tail(sort(apply(exp,1,sd)),1000)) #代码从里向外看，先是apply三个参数设置好，对每一行取标准差sd，然后是sort排序，最后是取最大的1000个基因名
n=exp[cg,]  #用cg作为行名，取子集,然后就得到了标准差最大的1000个基因，接下来用于画图

#绘制热图
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()

热图：可能是没有标准化吧，所以看起来奇奇怪怪的

取差异基因：

rm(list = ls()) 
load(file = "data/step2output.Rdata")
#差异分析，用limma包来做
#需要表达矩阵和group_list，不需要改
suppressMessages(library(limma))
design=model.matrix(~group_list)
fit=lmFit(exp,design)
fit=eBayes(fit)
deg=topTable(fit,coef=2,number = Inf)  #这个deg就是差异分析得到的结果
#class(deg)  #是一个data.frame
#dim(deg)    #22268  6

image.png

image.png

下面是做一些准备：

#为deg数据框添加几列为后面的可视化和富集分析做准备
#1.加probe_id列，把行名变成一列，与ids匹配做准备
#新增一列有两种方法：deg$probe_id=rownames(deg)，第2种方法是下面的mutate
suppressMessages(library(dplyr))
deg <- mutate(deg,probe_id=rownames(deg))
head(deg)  #取了交集之后，肯定会丢失一部分行（基因），就是没有注释的一些基因
#2.加symbol列，火山图要用
deg <- inner_join(deg,ids,by="probe_id")  #probe_id是共有的列名
#上面这一步用inner_join和merge都可以，但是用inner_join的好处是可以保留deg的顺序
head(deg)
#按照symbol列去重复：用duplicated根据逻辑值去重复，注意如果用unique函数会出现大片NA值
deg <- deg[!duplicated(deg$symbol),] 
#3.加change列,标记上下调基因
# 前两行是定义阈值：这里的改变会直接影响差异基因的数量
logFC_t=1
P.Value_t = 0.01
k1 = (deg$P.Value < P.Value_t)&(deg$logFC < -logFC_t)
k2 = (deg$P.Value < P.Value_t)&(deg$logFC > logFC_t)
change = ifelse(k1,
                "down",
                ifelse(k2,"up","stable"))
#统计一下差异基因的个数：table(change)
deg <- mutate(deg,change)   #把change列加到deg表中

之前是22268行，现在转换ID加上去重复之后变成12403行

image.png

最后再加上ENTREZID:

#4.加ENTREZID列，用于富集分析（symbol转entrezid，然后inner_join）
suppressMessages(library(ggplot2))
suppressMessages(library(clusterProfiler))
suppressMessages(library(org.Hs.eg.db))
s2e <- bitr(deg$symbol,    #bitr是转换的函数
            fromType = "SYMBOL",
            toType = "ENTREZID",
            OrgDb = org.Hs.eg.db)#人类
#class(s2e)  #是一个data.frame
#其他物种http://bioconductor.org/packages/release/BiocViews.html#___OrgDb
deg <- inner_join(deg,s2e,by=c("symbol"="SYMBOL"))

image.png

下面就是可视化环节了，激不激动!!走到这一步真的很不容易啊~~~😄

rm(list = ls()) 
load(file = "data/step1output.Rdata")
load(file = "data/step4output.Rdata")
#1.火山图----
suppressMessages(library(dplyr))
suppressMessages(library(ggplot2))
dat  = deg
for_label <- dat%>% 
  filter(abs(logFC)>logFC_t & P.Value < P.Value_t) %>%  #这一步是取差异基因
  head(10)  #取差异最显著的前10个基因
p <- ggplot(data = dat, 
            aes(x = logFC, 
                y = -log10(P.Value))) +
  geom_point(alpha=0.4, size=3.5, 
             aes(color=change)) +
  ylab("-log10(Pvalue)")+
  scale_color_manual(values=c("blue", "grey","red"))+
  geom_vline(xintercept=c(-logFC_t,logFC_t),lty=4,col="black",lwd=0.8) +
  geom_hline(yintercept = -log10(P.Value_t),lty=4,col="black",lwd=0.8) +
  theme_bw()
p

还是挺好看的呀哈哈哈~

#下面是一些进阶性的需求：在火山图上添加感兴趣的基因名
volcano_plot <- p +
  geom_point(size = 3, shape = 1, data = for_label) +
  ggrepel::geom_label_repel(
    aes(label = symbol),
    data = for_label,  #前面一步取的差异基因最显著的10个基因
    color="black"
  )
volcano_plot
ggsave(plot = volcano_plot,filename = paste0(gse,"volcano.png"))

添加了差异最显著的10个基因，也可以添加自己想添加的基因

#2.差异基因热图----
load(file = 'data/step2output.Rdata')
x=deg$logFC 
names(x)=deg$probe_id 
cg=c(names(head(sort(x),30)),names(tail(sort(x),30)))  #取30个上调和30个下调
#cg = deg$probe_id[deg$change !="stable"]   #这句是不想挑30，想把所有的都画出来，取不等于stable的所有行
n=exp[cg,]
dim(n)
#作热图
library(pheatmap)
annotation_col=data.frame(group=group_list)
rownames(annotation_col)=colnames(n) 
library(ggplotify)
heatmap_plot <- as.ggplot(pheatmap(n,show_colnames =F,
                                   show_rownames = F,
                                   scale = "row",
                                   #cluster_cols = F, 
                                   annotation_col=annotation_col)) 
heatmap_plot
ggsave(heatmap_plot,filename = paste0(gse,"heatmap.png"))

image.png

#下面是把两个图拼到一起
load("data/pca_plot.Rdata")
library(patchwork)
(pca_plot + volcano_plot +heatmap_plot)+ plot_annotation(tag_levels = "A")

#以3:1的比例来保存
ggsave(filename = "figures/test.png",width = 15,height = 5)

image.png

下面是富集分析：

rm(list = ls())  
load(file = 'data/step4output.Rdata')
#富集分析考验网速，因此给大家保存了Rdata
#上课运行示例数据无需修改，在做自己的数据时请注意把本行之后的load()去掉
library(clusterProfiler)
library(dplyr)
library(ggplot2)
source("kegg_plot_function.R")
#source的意思是在不打开某个脚本的情况下运行它，所以该脚本必须在工作目录下才可以运行
#source表示运行整个kegg_plot_function.R脚本，里面是一个function
#以up_kegg和down_kegg为输入数据做图

#1.GO database analysis ----

#(1)输入数据:富集分析的输入数据是EntrezID
gene_up = deg[deg$change == 'up','ENTREZID'] 
gene_down = deg[deg$change == 'down','ENTREZID'] 
gene_diff = c(gene_up,gene_down)
gene_all = deg[,'ENTREZID']
#(2)GO分析，分三部分
#以下步骤耗时很长，实际运行时注意把if后面的括号里F改成T
if(T){
  #细胞组分  CC
  ego_CC <- enrichGO(gene = gene_diff,  #输入数据是差异基因
                     OrgDb= org.Hs.eg.db,   #人类对应的数据库
                     ont = "CC",
                     pAdjustMethod = "BH",
                     minGSSize = 1,
                     pvalueCutoff = 0.01,
                     qvalueCutoff = 0.01,
                     readable = TRUE)
  #生物过程  BP
  ego_BP <- enrichGO(gene = gene_diff,
                     OrgDb= org.Hs.eg.db,
                     ont = "BP",
                     pAdjustMethod = "BH",
                     minGSSize = 1,
                     pvalueCutoff = 0.01,
                     qvalueCutoff = 0.01,
                     readable = TRUE)
  #分子功能：  MM
  ego_MF <- enrichGO(gene = gene_diff,
                     OrgDb= org.Hs.eg.db,
                     ont = "MF",
                     pAdjustMethod = "BH",
                     minGSSize = 1,
                     pvalueCutoff = 0.01,
                     qvalueCutoff = 0.01,
                     readable = TRUE)
  save(ego_CC,ego_BP,ego_MF,file = "data/ego_GSE14520.Rdata")
}
load(file = "data/ego_GSE14520.Rdata")
#class(ego_BP)  #数据类型是enrichResult,出自DOSE包
# tmp = ego_BP@result   #看PPT中的表
# View(tmp)

然后进行可视化：

#(3)可视化
#条带图
barplot(ego_CC,showCategory=20)
#气泡图
dotplot(ego_CC)   #不设置showCategory=几的话就默认展示10个

条带图（barplot）

气泡图（dotplot）

#下面的图需要映射颜色，设置和示例数据一样的geneList
# library(DOSE)
# data(geneList)   #这个是DOSE包中的示例数据，不用运行
geneList = deg$logFC
names(geneList)=deg$ENTREZID
geneList = sort(geneList,decreasing = T)   #排序是从大到小
#(3)展示top5通路的共同基因，要放大看。
#Gene-Concept Network
cnetplot(ego_CC, categorySize="pvalue", foldChange=geneList,colorEdge = TRUE)
cnetplot(ego_CC, foldChange=geneList, circular = TRUE, colorEdge = TRUE)
#Enrichment Map
emapplot(ego_CC)
#(4)展示通路关系
goplot(ego_CC)

cnetplot

image.png

emapplot

展示通路关系

#(5)Heatmap-like functional classification
heatplot(ego_CC,foldChange = geneList)
#太多基因就会糊。可通过调整比例或者减少基因来控制。
pdf("heatplot.pdf",width = 14,height = 5)
heatplot(ego_CC,foldChange = geneList)
dev.off()

因为基因太多了所以底下会糊，要适当调整比例

下面是KEGG:

# 2.KEGG pathway analysis----
#上调、下调、差异、所有基因
#（1）输入数据
gene_up = deg[deg$change == 'up','ENTREZID'] 
gene_down = deg[deg$change == 'down','ENTREZID'] 
gene_diff = c(gene_up,gene_down)
gene_all = deg[,'ENTREZID']
#（2）对上调/下调/所有差异基因进行富集分析
#注意这里又有个F
if(T){
  kk.up <- enrichKEGG(gene         = gene_up,
                      organism     = 'hsa',
                      universe     = gene_all,
                      pvalueCutoff = 0.9,
                      qvalueCutoff = 0.9)
  kk.down <- enrichKEGG(gene         =  gene_down,
                        organism     = 'hsa',
                        universe     = gene_all,
                        pvalueCutoff = 0.9,
                        qvalueCutoff =0.9)
  kk.diff <- enrichKEGG(gene         = gene_diff,
                        organism     = 'hsa',
                        pvalueCutoff = 0.9)
  save(kk.diff,kk.down,kk.up,file = "data/GSE14520kegg.Rdata")
}
load("data/GSE14520kegg.Rdata")
#(3)从富集结果中提取出结果数据框
kegg_diff_dt <- kk.diff@result

#(4)按照pvalue筛选通路
#在enrichkegg时没有设置pvaluecutoff，在此处筛选
down_kegg <- kk.down@result %>%
  filter(pvalue<0.05) %>% #筛选行
  mutate(group=-1) #新增列

up_kegg <- kk.up@result %>%
  filter(pvalue<0.05) %>%
  mutate(group=1)
#(5)可视化
g_kegg <- kegg_plot(up_kegg,down_kegg)

#g_kegg +scale_y_continuous(labels = c(20,15,10,5,0,5)) #手动修改坐标轴的刻度
ggsave(g_kegg,filename = 'figures/kegg_up_down.png')

image.png

#gsea作kegg富集分析，可选----
#(1)查看示例数据
data(geneList, package="DOSE")
#(2)将我们的数据转换成示例数据的格式
geneList=deg$logFC
names(geneList)=deg$ENTREZID
geneList=sort(geneList,decreasing = T)
#(3)富集分析
kk_gse <- gseKEGG(geneList     = geneList,
                  organism     = 'hsa',
                  nPerm        = 1000,
                  minGSSize    = 120,
                  pvalueCutoff = 0.9,
                  verbose      = FALSE)
down_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore < 0,];down_kegg$group=-1
up_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore > 0,];up_kegg$group=1
#(4)可视化
kegg_plot(up_kegg,down_kegg)
ggsave('figures/kegg_up_down_gsea.png')

kegg_up_down_gsea

基本上一篇完整的生信分析，从下载数据，处理数据，差异基因的分析，以及到后面的可视化包括PCA/火山图/热图，还有最后的富集分析的图，就算是一整套流程走下来了。真的很不容易啊！🤭

想想真正开始自己写简书也就这段时间，在宾馆隔离开始，然后5.1号开始告诉自己要每天写，然后现在回到宿舍也一直努力坚持着，希望自己可以一直坚持，写一些自己喜欢的东西，心情啊，想法啊，学习笔记什么的。

然后真正学习生信是放寒假开始，1月10号凌晨的飞机到新桥机场，然后就开始了我的寒假，因为自己课题需要就好好学习生信，学习R语言，学习各种数据库的使用，磕磕绊绊到现在，4个月，从一开始的代码是什么意思，频频报错，不知道怎么办，疯狂问别人，别人也不爱搭理，然后自己画图画到半夜，到现在可以从头到尾顺下来整个流程，虽然代码不是自己写的，但是代码什么意思应该是可以说出一二的，我还记得之前我问一个朋友怎么学的生信，他就说总会有阵痛期的，慢慢熬过去就好了，现在想来真的是这样。很多东西，只有自己经历了，去认真学习了，体会过，焦虑过，印象才最深刻。

刚开始逛CSDN帖子，看别人复现一篇又一篇文章，就觉得哇真的好厉害啊，我要是也可以这样就好了，这是2个月前的想法，由于自己一直不断努力，到现在，我也可以完整复现一篇文章，玩这个走完一套流程，自己学会了搜索，学会解决报错，最后出图，然后在这里记录我努力的过程，画图的结果，这样的过程是快乐的，欣慰的。

我会一直努力下去的❤加油💪猪莎🐷

2020.6.22 后面要更一波通过芯片数据集提取生存分析的数据从而做生存分析，先在这里立一个Flag~