数据挖掘流程代码版方便抄袭

　　GEO数据挖掘流程——代码版（方便抄袭）
　　微信公众号：生信小知识
　　关注可了解更多的教程及单细胞知识。问题或建议，请公众号留言;
　　内容目录
　　step0-install_packagesstep1-download_data1.背景知识2.关于GEO中的几个文件说明A. 2种family file(s)B. Series Matrix File(s)C. GSE64634_RAW.tar3.关于下载镜像问题4.关于探针id转换idmap1包——针对bioconductor包附加小功能——对基因名字进行注释（`annoGene`）idmap2包——万能芯片探针ID注释平台包（提取soft文件）idmap3包——idmap1和idmap2都无法注释的平台AnnoProbe包5. 整体代码这里抄step2-checkstep3-DEGstep4-anno-go-kegg真正代码detailed plot综合显示图（更推荐）step5-anno-GSEAstep6-anno-GSVAstep7-visualization致谢
　　step0-install_packages
　　简单安装R包和设置镜像代码：
　　1#镜像设置
　　2if(T){
　　3rm(list=ls())
　　4options()$repos
　　5options()$BioC_mirror
　　6options(BioC_mirror=＂https://mirrors.ustc.edu.cn/bioc/＂)
　　7options(＂repos＂=c(CRAN=＂https://mirrors.tuna.tsinghua.edu.cn/CRAN/＂))
　　8options()$repos
　　9options()$BioC_mirror
　　10}11:12#https://bioconductor.org/packages/release/bioc/html/GEOquery.html
　　13#BiocManager安装的R包
　　14if(T){
　　15if(!requireNamespace(＂BiocManager＂,quietly=TRUE))
　　16install.packages(＂BiocManager＂)17:18pkgs=c(＂KEGG.db＂,＂GSEABase＂,＂GSVA＂,＂clusterProfiler＂,＂GEOquery＂,＂limma＂,＂impute＂,＂genefu＂,＂org.Hs.eg.db＂,＂hgu133plus2.db＂)
　　19for(pkginpkgs){
　　20if(!requireNamespace(pkg,quietly=TRUE))
　　21BiocManager::install(pkg,ask=F,update=F)
　　22}
　　23}24:25#直接通过install.package函数安装的R包
　　26if(T){
　　27options()$repos
　　28pkgs=c(＂FactoMineR＂,＂factoextra＂,＂WGCNA＂,＂ggplot2＂,＂pheatmap＂,＂ggpubr＂)
　　29for(pkginpkgs){
　　30if(!requireNamespace(pkg,quietly=TRUE))
　　31install.packages(pkg)
　　32}
　　33}34:35
　　36#载入R包
　　37if(T){
　　38library(＂FactoMineR＂)
　　39library(＂factoextra＂)
　　40library(GSEABase)
　　41library(GSVA)
　　42library(clusterProfiler)
　　43library(genefu)
　　44library(ggplot2)
　　45library(ggpubr)
　　46library(hgu133plus2.db)
　　47library(limma)
　　48library(org.Hs.eg.db)
　　49library(pheatmap)
　　50}
　　step1-download_data
　　1.背景知识
　　这里通过使用GEOqueryR包来帮助我们下载数据，比手动登录GEO数据库一个个点击要方便很多！
　　而且我发现，最近NCBI更新后，GEO数据库也更新了！
　　点击后：
　　点击后：
　　可以看到官网的代码和我们目前用的代码基本一致。
　　2.关于GEO中的几个文件说明
　　下面一个个的来说吧！
　　A. 2种family file(s)
　　也就是SOFT formatted family file(s)和MINiML formatted family file(s)，通过字面，我们可以理解这是2种不同格式的探针说明文件。为什么我会这样说呢？因为我下载了soft文件来查看：
　　这时首行的内容：
　　就是告诉你这个数据是GPL570这个平台测序得到的：
　　中间有很多行关于GSM的，可能记录的是用到这个平台测序的sample：
　　接着就是一系列探针信息了：
　　我们可以通过在R种提取信息对我们得到的矩阵种探针做注释。
　　B. Series Matrix File(s)
　　这里面包含了3部分资料：数据属性+患者信息+表达矩阵
　　数据属性在最前面的几行，和患者信息之间有一个空行，但是他们都是以＂!＂开头：
　　接下来是患者信息：
　　紧接着患者信息下一行会有个提示：!series_matrix_table_begin，然后下面就是表达矩阵信息了：
　　这个就是我们需要的表达矩阵信息了，矩阵中每一行是一个基因（也就是一个探针），每一列是一个样本（GSM）
　　C. GSE64634_RAW.tar
　　这个我一开始不知道是什么东西，后来问了问jimmy，然后他给我发了个资料你要挖的公共数据集作者上传了错误的表达矩阵肿么办。大致瞅了瞅，知道这个应该是芯片的原始数据，然后也把这个文件给下载下来看了看：
　　根据常识猜了猜想：
　　X和Y应该代表着芯片上的位置，这个可能和探针有对应关系。
　　MEAN是信号的平均值，STDV是信号的标准差。
　　NPIXELS是像素点吧。
　　既然知道了这是原始数据，jimmy又给发了代码，感觉不学习下有点对不起自己，那就跟着代码过一遍吧：
　　1#BiocManager::install(c(＂oligo＂),ask=F,update=F)
　　2library(oligo)
　　3#BiocManager::install(c(＂pd.hg.u133.plus.2＂),ask=F,update=F)
　　4library(pd.hg.u133.plus.2)
　　5　　6#读入cel文件数据
　　7celFiles&lt;-list.celfiles(＂./＂,listGzipped=T)
　　8celFiles
　　9affyRaw&lt;-read.celfiles(celFiles)
　　10#提取矩阵并做normalization
　　11eset&lt;-rma(affyRaw)
　　12eset
　　13#检查数据
　　14dat=exprs(eset)
　　15dim(dat)
　　16boxplot(dat,las=2)
　　可以看到，整个过程非常的简单，就是利用了oligo这个R包而已。读入所有的cel文件后，利用rma()函数将数据进行了normalization，得到的是一个ExpressionSet对象！！
　　后面会比较我们利用rawdata得到的结果和我们直接下载的Series Matrix File(s)文件之间的区别！
　　3.关于下载镜像问题
　　jimmy曾经发表过GEO数据库中国区镜像横空出世推文，因为我每次下载都是用vpn，但是怕以后万一没有vpn了，所以还是学习下，以防万一嘛！再说了，可以学习下新知识，何乐而不为呢？
　　通过学习，发现实际上jimmy是开发了一个R包，或者说包装了一个函数吧，如果不想了解具体原理，那么用法如下：
　　1#安装方法1：
　　2library(devtools)
　　3install_github(＂jmzeng1314/GEOmirror＂)
　　4library(GEOmirror)
　　5#安装方法2：
　　6source(＂http://raw.githubusercontent.com/jmzeng1314/GEOmirror/master/R/geoChina.R＂)
　　7#安装方法3（源码）：
　　8geoChina&lt;-function(gse=＂GSE2546＂,mirror=＂tercent＂){
　　9suppressPackageStartupMessages(library(GEOquery))
　　10options(download.file.method=＂libcurl＂)
　　11#eSet=getGEO(＂GSE2546＂,destdir=＂.＂,AnnotGPL=F,getGPL=F)
　　12#http://49.235.27.111/GEOmirror/GSE2nnn/GSE2546_eSet.Rdata
　　13#gse=＂GSE2546＂;mirror=＂tercent＂
　　14gse=toupper(gse)
　　15down=ifelse(as.numeric(gsub(＂GSE＂,＂＂,gse))&lt;1000,
　　16paste0(＂/GEOmirror/GSEnnn/＂,gse,
　　17＂_eSet.Rdata＂),
　　18paste0(＂/GEOmirror/＂,
　　19gsub(＂[0-9][0-9][0-9]$＂,＂nnn＂,gse),＂/＂,gse,
　　20＂_eSet.Rdata＂))21:22if(mirror==＂tercent＂){
　　23up=＂http://49.235.27.111＂
　　24}
　　25tpf=paste0(gse,＂_eSet.Rdata＂)
　　26download.file(paste0(up,down),tpf,mode=＂wb＂)
　　27suppressWarnings(load(tpf))
　　28return(gset)
　　29}
　　具体用法也非常简单：
　　1eSet=geoChina(＂GSE2345＂)
　　通过安装方法3（源码）我们可以看出这个包的原理：
　　通过访问下载 http://49.235.27.111/GEOmirror/GSE2nnn/GSE2546_eSet.Rdata
　　所以可以猜到，应该是jimmy事先用循环的方式帮我们下载好了很多GEO数据，并做成了Rdata格式的文件！非常的良苦用心了！！
　　4.关于探针id转换
　　因为GEO数据矩阵中，横坐标都是探针的代号，如下图：
　　我们只看这些代号并不能理解具体是什么基因，于是这就需要我们做id转换了：将探针的代号转为基因symbol。
　　idmap1包——针对bioconductor包
　　这里又要提到jimmy＂开发＂的一个R包了：
　　第一个万能芯片探针ID注释平台R包——37个芯片平台
　　老规矩，还是先学习下：
　　一般重要的芯片在R的bioconductor里面都是有包的。但是需要我们单独下载对应的包。具体关系如下：
　　（具体的R包名称是bioc_package后加＂.db＂）
　　1gplbioc_packagetitle
　　2GPL201hgfocus[HG-Focus]AffymetrixHumanHG-FocusTargetArray
　　3GPL96hgu133a[HG-U133A]AffymetrixHumanGenomeU133AArray
　　4GPL571hgu133a2[HG-U133A_2]AffymetrixHumanGenomeU133A2.0Array
　　5GPL97hgu133b[HG-U133B]AffymetrixHumanGenomeU133BArray
　　6GPL570hgu133plus2[HG-U133_Plus_2]AffymetrixHumanGenomeU133Plus2.0Array
　　7GPL13667hgu219[HG-U219]AffymetrixHumanGenomeU219Array
　　8GPL8300hgu95av2[HG_U95Av2]AffymetrixHumanGenomeU95Version2Array
　　9GPL91hgu95av2[HG_U95A]AffymetrixHumanGenomeU95AArray
　　10GPL92hgu95b[HG_U95B]AffymetrixHumanGenomeU95BArray
　　11GPL93hgu95c[HG_U95C]AffymetrixHumanGenomeU95CArray
　　12GPL94hgu95d[HG_U95D]AffymetrixHumanGenomeU95DArray
　　13GPL95hgu95e[HG_U95E]AffymetrixHumanGenomeU95EArray
　　14GPL887hgug4110bAgilent-012097Human1AMicroarray(V2)G4110B(FeatureNumberversion)
　　15GPL886hgug4111aAgilent-011871Human1BMicroarrayG4111A(FeatureNumberversion)
　　16GPL1708hgug4112aAgilent-012391WholeHumanGenomeOligoMicroarrayG4112A(FeatureNumberversion)
　　17GPL13497HsAgilentDesign026652Agilent-026652WholeHumanGenomeMicroarray4x44Kv2(ProbeNameversion)
　　18GPL6244hugene10sttranscriptcluster[HuGene-1_0-st]AffymetrixHumanGene1.0STArray[transcript(gene)version]
　　19GPL11532hugene11sttranscriptcluster[HuGene-1_1-st]AffymetrixHumanGene1.1STArray[transcript(gene)version]
　　20GPL6097illuminaHumanv1Illuminahuman-6v1.0expressionbeadchip
　　21GPL6102illuminaHumanv2Illuminahuman-6v2.0expressionbeadchip
　　22GPL6947illuminaHumanv3IlluminaHumanHT-12V3.0expressionbeadchip
　　23GPL10558illuminaHumanv4IlluminaHumanHT-12V4.0expressionbeadchip
　　24GPL6885illuminaMousev2IlluminaMouseRef-8v2.0expressionbeadchip
　　25GPL81mgu74av2[MG_U74Av2]AffymetrixMurineGenomeU74AVersion2Array
　　26GPL82mgu74bv2[MG_U74Bv2]AffymetrixMurineGenomeU74BVersion2Array
　　27GPL83mgu74cv2[MG_U74Cv2]AffymetrixMurineGenomeU74Version2Array
　　28GPL339moe430a[MOE430A]AffymetrixMouseExpression430AArray
　　29GPL6246mogene10sttranscriptcluster[MoGene-1_0-st]AffymetrixMouseGene1.0STArray[transcript(gene)version]
　　30GPL340mouse4302[MOE430B]AffymetrixMouseExpression430BArray
　　31GPL1261mouse430a2[Mouse430_2]AffymetrixMouseGenome4302.0Array
　　32GPL8321mouse430a2[Mouse430A_2]AffymetrixMouseGenome430A2.0Array
　　因为很多时候，用户是找不到这些GPL平台对应的R包，或者下载安装困难，其实仅仅是需要探针与基因对应关系，没有必要去下载安装这几十个M的包。于是jimmy就开发了idmap1这个R包来帮我们：
　　安装方法（这里好像只能用方法1，因为在载入包的同时附带有一个变量p2s_df加载，如果用方法2和3没办法得到这个变量）：
　　也可以直接下载：https://codeload.github.com/jmzeng1314/idmap1/legacy.tar.gz/master
　　1#安装方法1：
　　2library(devtools)
　　3install_github(＂jmzeng1314/idmap1＂)
　　4library(idmap1)
　　5#安装方法2：
　　6source(＂https://raw.githubusercontent.com/jmzeng1314/idmap1/master/R/getIDs.R＂)
　　7#安装方法3（源码）：
　　8getIDs&lt;-function(gpl){
　　9#gpl=＂gpl570＂
　　10gpl=toupper(gpl)
　　11if(!gpl%in%unique(p2s_df$gpl)){
　　12stop(＂yourgplisnotinourannotationlist.＂)
　　13}
　　14return(p2s_df[p2s_df$gpl==gpl,])
　　15}
　　具体用法也非常简单：
　　1ids=getIDs(＂gpl570＂)
　　2head(ids)
　　关于我们得到的ExpressionSet对象，可以通过gset@annotation得到我们需要的注释平台信息。
　　前面说到了这个变量p2s_df，像我这么喜欢资源的人，当然也要保存一份在本地了。大家有兴趣的可以自己write到本地保存。哈哈哈哈哈哈哈哈哈哈（jimmy应该不会打我吧）！
　　附加小功能——对基因名字进行注释（`annoGene`）
　　同样的，还是idmap1这个R包里的函数，如果安装过这个R包就无须再安装了，如果没有安装又想用这个功能，还真的没有办法，因为这些数据存在于这个包的自带变量中（humanGTF、mouseGTF、ratGTF）:
　　1#安装方法2(无效)：
　　2source(＂https://raw.githubusercontent.com/jmzeng1314/idmap1/master/R/annoGene.R＂)
　　3#安装方法3（源码）(无效)：
　　4annoGene&lt;-function(IDs,ID_type,species=＂human＂,out_file){
　　5if(length(unique(IDs))&lt;1){
　　6stop(＂Youshouldgivemesomegenestobeannotated!!!＂)
　　7}
　　8if(!ID_type%in%c(＂ENSEMBL＂,＂SYMBOL＂)){
　　9stop(＂WeonlyacceptENSEMBLorSYMBOL!!!＂)
　　10}
　　11if(species==＂human＂){
　　12GTF&lt;-humanGTF
　　13}elseif(species==＂mouse＂){
　　14GTF&lt;-mouseGTF
　　15}elseif(species==＂rat＂){
　　16GTF&lt;-ratGTF
　　17}else{
　　18stop(＂Weonlyaccepthumanormouse,orrat,＂)
　　19}
　　20res&lt;-GTF[eval(parse(text=paste0(＂GTF$＂,ID_type)))%in%IDs,]21:22missIds&lt;-IDs[!(IDs%in%eval(parse(text=paste0(＂res$＂,ID_type))))]
　　23missIdsPercentage=round((length(missIds)/length(IDs))*100,2)
　　24if(length(missIds)!=0){
　　25warning(
　　26paste0(missIdsPercentage,＂%ofinputIDsarefailtoannotate...＂)
　　27#5.29%ofinputgeneIDsarefailtomap...
　　28)
　　29}
　　30if(hasArg(out_file)){
　　31results=res
　　32if(grepl(＂.html$＂,out_file)){
　　33Ensembl_prefix&lt;-＂https://asia.ensembl.org/Homo_sapiens/Gene/Summary?g=＂
　　34href=paste0(Ensembl_prefix,results$ENSEMBL)
　　35results$ENSEMBL=paste0(＂<strong><a target=＂＂_black＂＂ href=＂＂,＂>＂,results$ENSEMBL,＂</a></strong>＂)36:37symbol_prefix&lt;-＂http://www.ncbi.nlm.nih.gov/gene?term=＂
　　38href=paste0(symbol_prefix,results$SYMBOL)
　　39results$SYMBOL=paste0(＂<strong><a target=＂＂_black＂＂ href=＂＂,＂>＂,results$SYMBOL,＂</a></strong>＂)40:41y&lt;-DT::datatable(results,escape=F,rownames=F)
　　42DT::saveWidget(y,file=out_file)
　　43}elseif(grepl(＂.csv$＂,out_file)){
　　44write.csv(results,file=out_file)
　　45}else{
　　46stop(＂Weonlyacceptcsvorhtmlformat!!!＂)
　　47}48:49}
　　50return(res)
　　51}
　　这里补充学习了2个函数：
　　R语言parse函数与eval函数的字符串转命令行及执行操作：
　　parse()函数能将字符串转换为表达式expression；eval()函数能对表达式求解
　　idmap2包——万能芯片探针ID注释平台包（提取soft文件）
　　第二个万能芯片探针ID注释平台R包——有122个GPL
　　这次是根据[soft文件](#####A. 2种family file(s))进行提取信息得到的！
　　如果是我自己来处理这样的文件，我应该会分2步：
　　用linux中的grep -v命令将表头以＂#＂和＂！＂开头的行去掉
　　1grep-v^!GPL570.soft|grep-v^#|grep-v^^&gt;GPL570.txt
　　用R读入数据
　　1tmp=read.table(＂data/GPL570.txt＂,header=T,sep=＂ ＂,fill=T)
　　结果证明这样有错误，但是，具体原因有空再去找吧。下面看正确的读入方法——借助GEOqueryR包工具：
　　1library(GEOquery)
　　2gpl&lt;-getGEO(＂GPL570＂,destdir=＂data/＂)
　　3GPL=Table(gpl)
　　一般来说，大家关心的其实就是探针的ID，以及基因的symbol列。有了这个变量后，就可以按照R语言基本操作来提取我们需要的信息了。
　　注意：我检查了得到的结果，里面存在有的探针ID对应2个基因，如下图：
　　虽然不知道这些代表着什么意思，但是，我将这个数据和bioconductor包里的hgu133plus2.db数据做了比较，结果是这样的：自己提取的结果中如果是一个ID对应2个基因，那么这个探针在bioconductor包里基本上找不到数据。而其他一个ID对应2个基因的结果均和bioconductor包一致。
　　当然了，我这只是单独来看一个平台的探针，而在idmap2jimmy已经帮我们整理好了，直接用就行了！！
　　安装方法：
　　（比较慢）
　　也可以直接下载：https://codeload.github.com/jmzeng1314/idmap2/legacy.tar.gz/master
　　1library(devtools)
　　2install_github(＂jmzeng1314/idmap2＂)
　　3library(idmap2)
　　使用方法：
　　1library(idmap2)
　　2ids=get_soft_IDs(＂gpl570＂)
　　查看支持的平台：
　　1gpl_list[,1:4]#gpl_list是R包自带的变量
　　idmap3包——idmap1和idmap2都无法注释的平台
　　第三个万能芯片探针ID注释平台R包——idmap1和idmap2都无法注释的平台
　　如果我们拿到的soft注释文件中是序列信息，那么我们该怎么做呢？
　　应该是先将序列比对到参考基因组上，然后通过提取基因注释文件中的数据得到基因symbol！
　　img
　　而在idmap3包中，jimmy已经帮我们做好了！！说他宠粉也是真的了！！我都懒得做，他居然还写了个循环来完成了这个事。
　　安装方法：
　　（比较慢）
　　也可以直接下载：https://codeload.github.com/jmzeng1314/idmap3/legacy.tar.gz/master
　　1library(devtools)
　　2install_github(＂jmzeng1314/idmap3＂)
　　3library(idmap3)
　　使用方法：
　　1library(idmap3)
　　2ids=idmap3::get_pipe_IDs(＂GPL21827＂)
　　查看支持的平台：
　　1gpl_list[,1:4]#gpl_list是R包自带的变量
　　AnnoProbe包
　　芯片探针序列的基因注释已经无需你自己亲自做了——一个汇总
　　安装方法：
　　https://codeload.github.com/jmzeng1314/AnnoProbe/legacy.tar.gz/master
　　1library(devtools)
　　2install_github(＂jmzeng1314/AnnoProbe＂)
　　3library(AnnoProbe)
　　使用方法：
　　1gpl=＂GPL21827＂
　　2probe2gene=idmap(gpl,type=＂pipe＂)
　　5. 整体代码这里抄
　　1rm(list=ls())##魔幻操作，一键清空~
　　2options(stringsAsFactors=F)
　　3#下面代码只需要修改file的名字
　　4#下载的文件会自动保存到./data/下
　　5#得到的gset是一个ExpressionSet对象
　　6
　　7　　8#只需要修改file的名字，即可下载得到相应的geo数据
　　9#获取患者信息，需要自行修改
　　10file=＂GSE64634＂
　　11#https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE6463412:13if(T){
　　14#数据下载
　　15if(T){
　　16library(GEOquery)
　　17#这个包需要注意两个配置，一般来说自动化的配置是足够的。
　　18#Settingoptions(＂download.file.method.GEOquery＂=＂auto＂)
　　19#Settingoptions(＂GEOquery.inmemory.gpl＂=FALSE)
　　20fdata=paste0(file,＂_eSet.Rdata＂)
　　21fpath=paste0(＂data/＂,fdata)
　　22if(!file.exists(fpath)){
　　23gset&lt;-getGEO(file,destdir=＂data/＂,
　　24AnnotGPL=F,##注释文件
　　25getGPL=F)##平台文件
　　26gset=gset[[1]]
　　27save(gset,file=fpath)##保存到本地
　　28}
　　29load(fpath)
　　30}
　　31gset32:33#获取患者信息，这里需要自行修改
　　34if(T){
　　35pd=pData(gset)#根据diseasestate:ch1一列得知分组信息
　　36group_list=c(rep(＂normal＂,4),rep(＂npc＂,12))#nasopharyngealcarcinomaNPC鼻咽癌
　　37table(group_list)
　　38}39:40#对数据进行normalization
　　41if(T){
　　42dat=exprs(gset)
　　43dim(dat)44:45dat[1:4,1:4]
　　46boxplot(dat,las=2)
　　47dat=dat[apply(dat,1,sd)&gt;0,]#去除都是0的探针
　　48dat[dat&lt;0]=1
　　49boxplot(dat,las=2)50:51#dat=log2(dat+1)
　　52#boxplot(dat,las=2)
　　53library(limma)
　　54dat=normalizeBetweenArrays(dat)
　　55boxplot(dat,las=2)
　　56}57:58
　　59#探针注释2种方法，推荐方法2
　　60#方法1:比较麻烦，而且不方便，一般不用这种方法
　　61if(F){
　　62library(GEOquery)
　　63#DownloadGPLfile,putitinthecurrentdirectory,andloadit:
　　64gpl&lt;-getGEO(＂GPL570＂,destdir=＂data/＂)
　　65GPL=Table(gpl)
　　66colnames(GPL)
　　67head(GPL[,c(1,11)])##youneedtocheckthis,whichcolumndoyouneed
　　68probe2gene=GPL[,c(1,11)]
　　69head(probe2gene)
　　70save(probe2gene,file=＂probe2gene.Rdata＂)
　　71}72:73#方法2:用hgu133plus2.db这个R包比较方便
　　74if(T){
　　75library(hgu133plus2.db)
　　76ids=toTable(hgu133plus2SYMBOL)
　　77head(ids)
　　78ids=ids[ids$symbol!=＂＂,]
　　79ids=ids[ids$probe_id%in%rownames(dat),]#过滤没法注释的探针80:81dat[1:4,1:4]
　　82dat=dat[ids$probe_id,]#调整顺序，让dat的顺序和ids中的一致83:84ids$median=apply(dat,1,median)
　　85ids=ids[order(ids$symbol,ids$median,decreasing=T),]#按照基因名、中位数大小排序
　　86ids=ids[!duplicated(ids$symbol),]#只保留相同symbol中中位数最大的探针
　　87dat=dat[ids$probe_id,]#调整顺序，让dat的顺序和ids中的一致
　　88rownames(dat)=ids$symbol#id转换
　　89dat[1:4,1:4]
　　90}
　　91}92:93save(dat,group_list,file=＂data/step1-output.Rdata＂)
　　step2-check
　　画PCA图时要求是行名时样本名，列名时探针名，因此此时需要转置。格式要求data.frame。
　　关于scale：在基因表达量的数据中，有些基因表达量极低，有些基因表达量极高，因此把每个基因在不同处理和重复中的数据转换为平均值为0，方差为1的数据，所以在这里也需要先转置。
　　1#PCA检查
　　2#PCA，图会保存在pic/all_samples_PCA.png
　　3if(T){
　　4#载入数据
　　5if(T){
　　6rm(list=ls())##魔幻操作，一键清空~
　　7options(stringsAsFactors=F)8:9load(file=＂data/step1-output.Rdata＂)
　　10print(table(group_list))
　　11#每次都要检测数据
　　12dat[1:4,1:4]
　　13}14:15#PCA，图会保存在pic/all_samples_PCA.png
　　16if(T){
　　17exprSet=dat
　　18#过滤标准
　　19if(T){
　　20dim(exprSet)
　　21#过滤标准需要修改,目前是保留每一个基因在&gt;5个人中表达量&gt;1
　　22exprSet=exprSet[apply(exprSet,1,function(x)sum(x&gt;1)&gt;5),]
　　23boxplot(apply(exprSet,1,mean))
　　24dim(exprSet)
　　25#过滤标准需要修改,目前是去除每一个基因在&gt;5个人中表达量&gt;12
　　26exprSet=exprSet[!apply(exprSet,1,function(x)sum(x&gt;12)&gt;5),]
　　27dim(exprSet)
　　28}
　　29#去除文库大小差异normalization，同时利用log做scale
　　30exprSet=log(edgeR::cpm(exprSet)+1)
　　31dat=exprSet
　　32dat=as.data.frame(t(dat))#画PCA图时要求是行名时样本名，列名时探针名，因此此时需要转换。格式要求data.frame
　　33library(＂FactoMineR＂)#计算PCA
　　34library(＂factoextra＂)#画图展示35:36dat.pca&lt;-PCA(dat,graph=F)
　　37#fviz_pca_ind按样本fviz_pca_var按基因
　　38fviz_pca_ind(dat.pca,
　　39geom.ind=＂point＂,#c(＂point＂,＂text)2选1
　　40col.ind=group_list,#colorbygroups
　　41#palette=c(＂#00AFBB＂,＂#E7B800＂),#自定义颜色
　　42addEllipses=T,#加圆圈
　　43legend.title=＂Groups＂#图例名称
　　44)
　　45ggsave(＂pic/all_samples_PCA.png＂)
　　46}
　　47}48:49
　　50#挑选1000个SD最大的基因画表达量热图
　　51#结果图存放在pic/heatmap_top1000_sd.png中
　　52if(T){
　　53#载入数据
　　54if(T){
　　55rm(list=ls())
　　56load(file=＂data/step1-output.Rdata＂)
　　57dat[1:4,1:4]
　　58}59:60#挑选1000个SD最大的基因画表达量热图
　　61#结果图存放在pic/heatmap_top1000_sd.png中
　　62if(T){
　　63cg=names(tail(sort(apply(dat,1,sd)),1000))#找到SD最大的1000个基因
　　64library(pheatmap)
　　65headmap.dat=dat[cg,]
　　66#把每个基因在不同处理和重复中的数据转换为平均值为0，
　　67#方差为1的数据，所以在这里也需要先转置(针对基因)。
　　68headmap.dat=t(scale(t(headmap.dat)))
　　69headmap.dat[headmap.dat&gt;2]=2
　　70headmap.dat[headmap.dat&lt;-2]=-271:72#分组信息设置
　　73ac=data.frame(group=group_list)
　　74rownames(ac)=colnames(headmap.dat)75:76##可以看到TNBC具有一定的异质性，拿它来区分乳腺癌亚型指导临床治疗还是略显粗糙。
　　77pheatmap(headmap.dat,show_colnames=T,show_rownames=F,
　　78annotation_col=ac,
　　79filename=＂pic/heatmap_top1000_sd.png＂)
　　80}
　　81}82:83
　　84#过滤标准需要修改
　　85#利用top500_mad基因画相关性热图热图
　　86#结果图存放在pic/cor_top500_mad.png中
　　87if(T){
　　88#导入数据
　　89if(T){
　　90rm(list=ls())
　　91load(file=＂data/step1-output.Rdata＂)
　　92dat[1:4,1:4]
　　93exprSet=dat
　　94}95:96#利用所有基因画相关性热图
　　97if(T){
　　98ac=data.frame(group=group_list)
　　99rownames(ac)=colnames(exprSet)
　　100pheatmap::pheatmap(cor(exprSet),
　　101annotation_col=ac,
　　102show_rownames=F,
　　103filename=＂pic/cor_all_genes.png＂)
　　104}105:106#过滤标准
　　107if(T){
　　108dim(exprSet)
　　109#过滤标准需要修改,目前是保留每一个基因在&gt;5个人中表达量&gt;1
　　110exprSet=exprSet[apply(exprSet,1,function(x)sum(x&gt;1)&gt;5),]
　　111boxplot(apply(exprSet,1,mean))
　　112dim(exprSet)
　　113#过滤标准需要修改,目前是去除每一个基因在&gt;5个人中表达量&gt;12
　　114exprSet=exprSet[!apply(exprSet,1,function(x)sum(x&gt;12)&gt;5),]
　　115dim(exprSet)
　　116}117:118#数据normalization处理
　　119if(T){
　　120#去除文库大小差异normalization，同时利用log做scale
　　121exprSet=log(edgeR::cpm(exprSet)+1)
　　122#取top500mad的基因画图
　　123exprSet=exprSet[names(sort(apply(exprSet,1,mad),decreasing=T)[1:500]),]
　　124}125:126#利用top500_mad基因画相关性热图热图
　　127if(T){
　　128M=cor(exprSet)
　　129pheatmap::pheatmap(M,
　　130show_rownames=F,
　　131annotation_col=ac,
　　132filename=＂pic/cor_top500_mad.png＂)
　　133}
　　134}
　　all_samples_PCA.png
　　heatmap_top1000_sd.png
　　cor_top500_mad.png
　　step3-DEG
　　关于差异分析是否需要比较矩阵
　　差异分析的统计学本质探究
　　原始版火山图：plot(logFC,-log10(P.Value))
　　原始版MA图：plot(AveExpr,logFC)
　　1#载入数据,检查数据
　　2if(T){
　　3rm(list=ls())
　　4options(stringsAsFactors=F)
　　5library(ggpubr)
　　6load(file=＂data/step1-output.Rdata＂)
　　7#每次都要检测数据
　　8dat[1:4,1:4]
　　9table(group_list)
　　10boxplot(dat[1,]~group_list)#按照group_list分组画箱线图
　　11　　12#boxplot的美化版
　　13bplot=function(g){
　　14df=data.frame(gene=g,stage=group_list)
　　15p&lt;-ggboxplot(df,x=＂stage＂,y=＂gene＂,
　　16color=＂stage＂,palette=＂jco＂,
　　17add=＂jitter＂)
　　18#Addp-value
　　19p+stat_compare_means()
　　20}
　　21}
　　22#利用定义好的函数检查数据
　　23bplot(dat[1,])
　　24bplot(dat[2,])
　　25bplot(dat[3,])
　　26bplot(dat[4,])
　　27
　　28　　29#limma
　　30library(limma)
　　31#方法1：不制作比较矩阵，简单
　　32#但是做不到随心所欲的指定任意两组进行比较
　　33if(T){
　　34design=model.matrix(~factor(group_list))
　　35fit=lmFit(dat,design)
　　36fit2=eBayes(fit)
　　37##上面是limma包用法的一种方式
　　38options(digits=4)#设置全局的数字有效位数为4
　　39topTable(fit2,coef=2,adjust=＂BH＂)
　　40}
　　41　　42#方法2：制作比较矩阵
　　43#可以随心所欲的指定任意两组进行比较
　　44if(T){
　　45design&lt;-model.matrix(~0+factor(group_list))
　　46colnames(design)=levels(factor(group_list))
　　47head(design)
　　48exprSet=dat
　　49rownames(design)=colnames(exprSet)
　　50design
　　51　　52#比较矩阵
　　53#这个矩阵声明，我们要把npc组跟Normal进行差异分析比较
　　54contrast.matrix&lt;-makeContrasts(＂npc-normal＂,
　　55levels=design)
　　56contrast.matrix
　　57　　58deg=function(exprSet,design,contrast.matrix){
　　59##step1
　　60fit&lt;-lmFit(exprSet,design)
　　61##step2
　　62fit2&lt;-contrasts.fit(fit,contrast.matrix)
　　63　　64fit2&lt;-eBayes(fit2)##defaultnotrend!!!
　　65##eBayes()withtrend=TRUE
　　66　　67##step3
　　68#有了比较矩阵后，coef=1，而number=Inf是把所有结果都打印出来
　　69tempOutput=topTable(fit2,coef=1,number=Inf)
　　70nrDEG=na.omit(tempOutput)
　　71#write.csv(nrDEG2,＂limma_notrend.results.csv＂,quote=F)
　　72head(nrDEG)
　　73return(nrDEG)
　　74}
　　75　　76deg=deg(exprSet,design,contrast.matrix)
　　77　　78head(deg)
　　79}
　　80
　　81　　82save(deg,file=＂data/deg.Rdata＂)
　　83　　84load(file=＂data/deg.Rdata＂)
　　85head(deg)
　　86bplot(dat[rownames(deg)[1],])
　　87##forvolcanoandMAplot
　　88#结果存放在pic/volcano.png和pic/MA.png
　　89if(T){
　　90nrDEG=deg
　　91head(nrDEG)
　　92attach(nrDEG)
　　93#原始版火山图
　　94plot(logFC,-log10(P.Value))
　　95library(ggpubr)
　　96df=nrDEG
　　97df$y=-log10(P.Value)
　　98ggscatter(df,x=＂logFC＂,y=＂y＂,size=0.5)
　　99#定义logFC=2为阈值
　　100df$state=ifelse(df$P.Value&gt;0.01,＂stable＂,
　　101ifelse(df$logFC&gt;2,＂up＂,
　　102ifelse(df$logFC&lt;-2,＂down＂,＂stable＂))
　　103)
　　104table(df$state)
　　105df$name=rownames(df)
　　106head(df)
　　107ggscatter(df,x=＂logFC＂,y=＂y＂,size=0.5,color=＂state＂)
　　108ggscatter(df,x=＂logFC＂,y=＂y＂,color=＂state＂,size=0.5,
　　109label=＂name＂,repel=T,
　　110#label.select=rownames(df)[df$state!=＂stable＂],
　　111label.select=c(＂TTC9＂,＂AQP3＂,＂CXCL11＂,＂PTGS2＂),#挑选一些基因在图中显示出来
　　112palette=c(＂#00AFBB＂,＂#E7B800＂,＂#FC4E07＂))
　　113ggsave(＂pic/volcano.png＂)114:115#MA图
　　116ggscatter(df,x=＂AveExpr＂,y=＂logFC＂,size=0.2)
　　117df$p_c=ifelse(df$P.Value&lt;0.001,＂p&lt;0.001＂,
　　118ifelse(df$P.Value&lt;0.01,＂0.001&lt;p&lt;0.01＂,＂p&gt;0.01＂))<!--0.01＂,＂p-->
　　119table(df$p_c)
　　120ggscatter(df,x=＂AveExpr＂,y=＂logFC＂,color=＂p_c＂,size=0.2,
　　121palette=c(＂green＂,＂red＂,＂black＂))
　　122ggsave(＂pic/MA.png＂)
　　123}124:125##forheatmap
　　126if(T){
　　127load(file=＂data/step1-output.Rdata＂)
　　128#每次都要检测数据
　　129dat[1:4,1:4]
　　130table(group_list)
　　131x=deg$logFC
　　132names(x)=rownames(deg)133:134#cg中存放着变化上升和下降的前100个基因名
　　135cg=c(names(head(sort(x),100)),
　　136names(tail(sort(x),100)))
　　137library(pheatmap)
　　138n=t(scale(t(dat[cg,])))139:140n[n&gt;2]=2
　　141n[n&lt;-2]=-2
　　142n[1:4,1:4]
　　143pheatmap(n,show_colnames=F,show_rownames=F)
　　144ac=data.frame(group=group_list)
　　145rownames(ac)=colnames(n)#将ac的行名也就分组信息（是‘noTNBC’还是‘TNBC’）给到n的列名，即热图中位于上方的分组信息
　　146pheatmap(n,show_colnames=F,
　　147show_rownames=F,
　　148cluster_cols=F,
　　149annotation_col=ac,filename=＂pic/heatmap_top200_DEG.png＂)#列名注释信息为ac即分组信息150:151
　　152}153:154write.csv(deg,file=＂data/deg.csv＂)
　　volcano.png
　　MA.png
　　heatmap_top200_DEG.png
　　step4-anno-go-kegg
　　关于基因名和ID之间的各种转换：bitr(geneID, fromType, toType, OrgDb, drop = TRUE)例如：
　　1head(unique(deg$symbol))
　　2#[1]＂LOC388780＂＂SLC6A6＂＂RGS22＂＂AKR1D1＂＂AOC1＂＂FOXJ1＂
　　3bitr(unique(deg$symbol),fromType=＂SYMBOL＂,
　　4toType=c(＂ENTREZID＂),
　　5OrgDb=org.Hs.eg.db)
　　这里将KEGG和GO富集分析函数自定义为3个函数——run_kegg、run_go、kegg_plot
　　每次运行前先运行下面代码：
　　1#KEGGpathwayanalysis
　　2#具体结果数据在data/下，图在pic/下
　　3run_kegg&lt;-function(gene_up,gene_down,geneList=F,pro=＂test＂){
　　4gene_up=unique(gene_up)
　　5gene_down=unique(gene_down)
　　6gene_diff=unique(c(gene_up,gene_down))
　　7###over-representationtest
　　8#下面把3个基因集分开做超几何分布检验
　　9#首先是上调基因集。
　　10kk.up&lt;-enrichKEGG(gene=gene_up,
　　11organism=＂hsa＂,
　　12#universe=gene_all,
　　13pvalueCutoff=0.9,
　　14qvalueCutoff=0.9)
　　15head(kk.up)[,1:6]
　　16kk=kk.up
　　17dotplot(kk)
　　18barplot(kk)
　　19kk=DOSE::setReadable(kk,OrgDb=＂org.Hs.eg.db＂,keyType=＂ENTREZID＂)
　　20head(kk)
　　21write.csv(kk@result,paste0(＂data/＂,pro,＂_kk.up.csv＂))
　　22　　23#首先是下调基因集。
　　24kk.down&lt;-enrichKEGG(gene=gene_down,
　　25organism=＂hsa＂,
　　26#universe=gene_all,
　　27pvalueCutoff=0.9,
　　28qvalueCutoff=0.9)
　　29head(kk.down)[,1:6]
　　30kk=kk.down
　　31dotplot(kk)
　　32barplot(kk)
　　33kk=DOSE::setReadable(kk,OrgDb=＂org.Hs.eg.db＂,keyType=＂ENTREZID＂)
　　34write.csv(kk@result,paste0(＂data/＂,pro,＂_kk.down.csv＂))
　　35　　36#最后是上下调合并后的基因集。
　　37kk.diff&lt;-enrichKEGG(gene=gene_diff,
　　38organism=＂hsa＂,
　　39pvalueCutoff=0.05)
　　40head(kk.diff)[,1:6]
　　41kk=kk.diff
　　42dotplot(kk)
　　43barplot(kk)
　　44kk=DOSE::setReadable(kk,OrgDb=＂org.Hs.eg.db＂,keyType=＂ENTREZID＂)
　　45write.csv(kk@result,paste0(＂data/＂,pro,＂_kk.diff.csv＂))
　　46
　　47　　48kegg_down_dt&lt;-as.data.frame(kk.down)
　　49kegg_up_dt&lt;-as.data.frame(kk.up)
　　50down_kegg&lt;-kegg_down_dt[kegg_down_dt$pvalue&lt;0.01,];down_kegg$group=-1
　　51up_kegg&lt;-kegg_up_dt[kegg_up_dt$pvalue&lt;0.01,];up_kegg$group=1
　　52#画图设置, 这个图很丑，大家可以自行修改。
　　53kegg_plot(up_kegg,down_kegg)
　　54　　55ggsave(g_kegg,filename=paste0(＂pic/＂,pro,＂_kegg_up_down.png＂))
　　56　　57###GSEA
　　58if(geneList){
　　59　　60## GSEA算法跟上面的使用差异基因集做超几何分布检验不一样。
　　61kk_gse&lt;-gseKEGG(geneList=geneList,
　　62organism=＂hsa＂,
　　63nPerm=1000,
　　64minGSSize=20,
　　65pvalueCutoff=0.9,
　　66verbose=FALSE)
　　67head(kk_gse)[,1:6]
　　68gseaplot(kk_gse,geneSetID=rownames(kk_gse[1,]))
　　69gseaplot(kk_gse,＂hsa04110＂,title=＂Cellcycle＂)
　　70kk=DOSE::setReadable(kk_gse,OrgDb=＂org.Hs.eg.db＂,keyType=＂ENTREZID＂)
　　71tmp=kk@result
　　72write.csv(kk@result,paste0(pro,＂_kegg.gsea.csv＂))
　　73
　　74　　75#这里找不到显著下调的通路，可以选择调整阈值，或者其它。
　　76down_kegg&lt;-kk_gse[kk_gse$pvalue&lt;0.05&amp;kk_gse$enrichmentScore&lt;0,];down_kegg$group=-1
　　77up_kegg&lt;-kk_gse[kk_gse$pvalue&lt;0.05 enrichmentscore=＂＂&gt;0,];up_kegg$group=1<!--0.05-->
　　78　　79g_kegg=kegg_plot(up_kegg,down_kegg)
　　80print(g_kegg)
　　81ggsave(g_kegg,filename=paste0(pro,＂_kegg_gsea.png＂))
　　82}
　　83}
　　84　　85#GOdatabaseanalysis
　　86#具体结果数据在data/下，图在pic/下
　　87run_go&lt;-function(gene_up,gene_down,pro=＂test＂){
　　88gene_up=unique(gene_up)
　　89gene_down=unique(gene_down)
　　90gene_diff=unique(c(gene_up,gene_down))
　　91g_list=list(gene_up=gene_up,
　　92gene_down=gene_down,
　　93gene_diff=gene_diff)
　　94#因为go数据库非常多基因集，所以运行速度会很慢。
　　95if(T){
　　96go_enrich_results&lt;-lapply(g_list,function(gene){
　　97lapply(c(＂BP＂,＂MF＂,＂CC＂),function(ont){
　　98cat(paste(＂Nowprocess＂,names(gene),ont))
　　99ego&lt;-enrichGO(gene=gene,
　　100#universe=gene_all,
　　101OrgDb=org.Hs.eg.db,
　　102ont=ont,
　　103pAdjustMethod=＂BH＂,
　　104pvalueCutoff=0.99,
　　105qvalueCutoff=0.99,
　　106readable=TRUE)107:108print(head(ego))
　　109return(ego)
　　110})
　　111})
　　112save(go_enrich_results,file=paste0(＂data/＂,pro,＂_go_enrich_results.Rdata＂))113:114}
　　115#只有第一次需要运行，就保存结果哈，下次需要探索结果，就载入即可。
　　116#load(file=paste0(＂data/＂,pro,＂_go_enrich_results.Rdata＂))117:118n1=c(＂gene_up＂,＂gene_down＂,＂gene_diff＂)
　　119n2=c(＂BP＂,＂MF＂,＂CC＂)
　　120for(iin1:3){
　　121for(jin1:3){
　　122fn=paste0(＂pic/＂,pro,＂_dotplot_＂,n1[i],＂_＂,n2[j],＂.png＂)
　　123cat(paste0(fn,＂＂))
　　124png(fn,res=150,width=1080)
　　125print(dotplot(go_enrich_results[[i]][[j]]))
　　126dev.off()
　　127}
　　128}
　　129}130:131132:133#合并up和down的基因KEGG结果，放在一幅图里展示
　　134kegg_plot&lt;-function(up_kegg,down_kegg){
　　135dat=rbind(up_kegg,down_kegg)
　　136colnames(dat)
　　137dat$pvalue=-log10(dat$pvalue)
　　138dat$pvalue=dat$pvalue*dat$group139:140dat=dat[order(dat$pvalue,decreasing=F),]141:142g_kegg&lt;-ggplot(dat,aes(x=reorder(Description,order(pvalue,decreasing=F)),y=pvalue,fill=group))+
　　143geom_bar(stat=＂identity＂)+
　　144scale_fill_gradient(low=＂blue＂,high=＂red＂,guide=FALSE)+
　　145scale_x_discrete(name=＂Pathwaynames＂)+
　　146scale_y_continuous(name=＂log10P-value＂)+
　　147coord_flip()+theme_bw()+theme(plot.title=element_text(hjust=0.5))+
　　148ggtitle(＂PathwayEnrichment＂)
　　149print(g_kegg)
　　150}
　　真正代码
　　1#载入数据
　　2if(T){
　　3rm(list=ls())
　　4options(stringsAsFactors=F)5:6load(file=＂data/deg.Rdata＂)
　　7head(deg)
　　8}9:10#数据预处理
　　11##不同的阈值，筛选到的差异基因数量就不一样，后面的超几何分布检验结果就大相径庭。
　　12logFC_t=1.5
　　13#预处理1
　　14if(T){
　　15deg$g=ifelse(deg$P.Value&gt;0.05,＂stable＂,
　　16ifelse(deg$logFC&gt;logFC_t,＂UP＂,
　　17ifelse(deg$logFC&lt;-logFC_t,＂DOWN＂,＂stable＂))
　　18)
　　19table(deg$g)
　　20head(deg)
　　21deg$symbol=rownames(deg)
　　22library(ggplot2)
　　23library(clusterProfiler)
　　24library(org.Hs.eg.db)
　　25df&lt;-bitr(unique(deg$symbol),fromType=＂SYMBOL＂,
　　26toType=c(＂ENTREZID＂),
　　27OrgDb=org.Hs.eg.db)
　　28head(df)
　　29DEG=deg
　　30head(DEG)31:32DEG=merge(DEG,df,by.y=＂SYMBOL＂,by.x=＂symbol＂)
　　33head(DEG)
　　34save(DEG,file=＂data/anno_DEG.Rdata＂)
　　35}
　　36#预处理2
　　37if(T){
　　38gene_up=DEG[DEG$g==＂UP＂,＂ENTREZID＂]
　　39gene_down=DEG[DEG$g==＂DOWN＂,＂ENTREZID＂]
　　40gene_diff=c(gene_up,gene_down)
　　41gene_all=as.character(DEG[,＂ENTREZID＂])
　　42data(geneList,package=＂DOSE＂)
　　43head(geneList)
　　44boxplot(geneList)
　　45boxplot(DEG$logFC)46:47geneList=DEG$logFC
　　48names(geneList)=DEG$ENTREZID
　　49geneList=sort(geneList,decreasing=T)
　　50}51:52#detailedplot
　　53if(T){
　　54source(＂kegg_and_go_up_and_down.R＂)
　　55run_kegg(gene_up,gene_down,pro=＂npc_VS_normal＂)
　　56#需要多go数据库的3个条目进行3次富集分析，非常耗时。
　　57run_go(gene_up,gene_down,pro=＂npc_VS_normal＂)
　　58}59:60
　　61#综合显示图
　　62if(F){
　　63go&lt;-enrichGO(gene_up,OrgDb=＂org.Hs.eg.db＂,ont=＂all＂)
　　64library(ggplot2)
　　65library(stringr)
　　66barplot(go,split=＂ONTOLOGY＂)+facet_grid(ONTOLOGY~.,scale=＂free＂)
　　67barplot(go,split=＂ONTOLOGY＂,font.size=10)+
　　68facet_grid(ONTOLOGY~.,scale=＂free＂)+
　　69scale_x_discrete(labels=function(x)str_wrap(x,width=50))+
　　70ggsave(＂pic/gene_up_GO_all_barplot.png＂)71:72go&lt;-enrichGO(gene_down,OrgDb=＂org.Hs.eg.db＂,ont=＂all＂)
　　73barplot(go,split=＂ONTOLOGY＂,font.size=10)+
　　74facet_grid(ONTOLOGY~.,scale=＂free＂)+
　　75scale_x_discrete(labels=function(x)str_wrap(x,width=50))+
　　76ggsave(＂pic/gene_down_GO_all_barplot.png＂)
　　77}
　　detailed plot
　　npc_VS_normal_kegg_up_down.png（kegg结果还可以接受）
　　GO系列结果过于冗余：
　　npc_VS_normal_dotplot_gene_diff_BP.png
　　npc_VS_normal_dotplot_gene_diff_CC.png
　　npc_VS_normal_dotplot_gene_diff_MF.png
　　npc_VS_normal_dotplot_gene_down_BP.png
　　npc_VS_normal_dotplot_gene_down_CC
　　npc_VS_normal_dotplot_gene_down_MF.png
　　npc_VS_normal_dotplot_gene_up_BP.png
　　npc_VS_normal_dotplot_gene_up_CC.png
　　npc_VS_normal_dotplot_gene_up_MF.png
　　综合显示图（更推荐）
　　gene_down_GO_all_barplot.png
　　gene_up_GO_all_barplot.png
　　step5-anno-GSEA
　　GSEA数据下载网页：https://www.gsea-msigdb.org/gsea/downloads.jsp
　　msigdb.v7.0.symbols.gmt：https://www.gsea-msigdb.org/gsea/msigdb/download_file.jsp?filePath=/resources/msigdb/7.0/msigdb.v7.0.symbols.gmt
　　msigdb.v7.0.entrez.gmt：https://www.gsea-msigdb.org/gsea/msigdb/download_file.jsp?filePath=/resources/msigdb/7.0/msigdb.v7.0.entrez.gmt
　　所有打包gmt：https://www.gsea-msigdb.org/gsea/msigdb/download_file.jsp?filePath=/resources/msigdb/7.0/msigdb_v7.0_files_to_download_locally.zip
　　1#载入数据和R包
　　2#DEG为limma得到的差异分析结果
　　3if(T){
　　4rm(list=ls())
　　5options(stringsAsFactors=F)
　　6load(file=＂data/anno_DEG.Rdata＂)
　　7library(ggplot2)
　　8library(clusterProfiler)
　　9library(org.Hs.eg.db)
　　10}11:12
　　13###对 MigDB中的全部基因集做GSEA分析。
　　14###按照FC的值对差异基因进行排序
　　15#http://www.bio-info-trainee.com/2105.html
　　16#http://www.bio-info-trainee.com/2102.html
　　17#自行修改存放gmt文件路径d
　　18#GSEA每个geneset的具体结果保存在gsea_results这个list中
　　19#而最终结果保存在gsea_results_df数据框中
　　20d=＂data/GSEA_gmt/msigdb_v7.0_files_to_download_locally/msigdb_v7.0_GMTs/symbols/＂
　　21if(T){
　　22geneList=DEG$logFC
　　23names(geneList)=DEG$symbol
　　24geneList=sort(geneList,decreasing=T)
　　25#选择gmt文件（MigDB中的全部基因集）26:27gmts=list.files(d,pattern=＂all＂)
　　28gmts29:30#GSEA分析
　　31library(GSEABase)#BiocManager::install(＂GSEABase＂)
　　32##下面使用lapply循环读取每个gmt文件，并且进行GSEA分析
　　33##如果存在之前分析后保存的结果文件，就不需要重复进行GSEA分析。
　　34f=＂data/gsea_results.Rdata＂
　　35if(!file.exists(f)){
　　36gsea_results&lt;-lapply(gmts,function(gmtfile){
　　37#gmtfile=gmts[2]
　　38filepath=paste0(d,gmtfile)
　　39geneset&lt;-read.gmt(filepath)
　　40print(paste0(＂Nowprocessthe＂,gmtfile))
　　41egmt&lt;-GSEA(geneList,TERM2GENE=geneset,verbose=FALSE)
　　42head(egmt)
　　43#gseaplot(egmt,geneSetID=rownames(egmt[1,]))
　　44return(egmt)
　　45})
　　46#上面的代码耗时，所以保存结果到本地文件
　　47save(gsea_results,file=f)
　　48}
　　49load(file=f)
　　50#提取gsea结果，熟悉这个对象
　　51gsea_results_list&lt;-lapply(gsea_results,function(x){
　　52cat(paste(dim(x@result)),＂＂)
　　53x@result
　　54})
　　55}
　　56#随便看几个结果图
　　57dev.off()
　　58gsea_results_df&lt;-do.call(rbind,gsea_results_list)
　　59gseaplot(gsea_results[[2]],geneSetID=＂NIKOLSKY_BREAST_CANCER_7P15_AMPLICON＂)
　　60gseaplot(gsea_results[[2]],＂RICKMAN_HEAD_AND_NECK_CANCER_D＂)
　　step6-anno-GSVA
　　老实说，我对GSVA还不是很理解，但是知道在代码方面，做起来比较简单，有2个输入文件，一个是表达矩阵，一个是gmt文件即可。先把代码谢在这里，日后如果有需要，再仔细研究吧：
　　1###对 MigDB中的全部基因集做GSVA分析。
　　2##还有ssGSEA,PGSEA
　　3#载入数据
　　4if(T){
　　5rm(list=ls())
　　6options(stringsAsFactors=F)7:8library(ggplot2)
　　9library(clusterProfiler)
　　10library(org.Hs.eg.db)
　　11load(file=＂data/step1-output.Rdata＂)
　　12#每次都要检测数据
　　13dat[1:4,1:4]
　　14}15:16#GSVA分析
　　17#存放geneset的文件路径需要具体修改
　　18d=＂data/GSEA_gmt/msigdb_v7.0_files_to_download_locally/msigdb_v7.0_GMTs/symbols/＂
　　19if(T){
　　20X=dat
　　21table(group_list)
　　22##MolecularSignaturesDatabase(MSigDb)
　　23gmts=list.files(d,pattern=＂all＂)
　　24gmts
　　25library(GSVA)#BiocManager::install(＂GSVA＂)
　　26if(!file.exists(＂data/gsva_msigdb.Rdata＂)){
　　27es_max&lt;-lapply(gmts,function(gmtfile){
　　28#gmtfile=gmts[8];gmtfile
　　29geneset&lt;-read.gmt(file.path(d,gmtfile))
　　30es.max&lt;-gsva(X,geneset,
　　31mx.diff=FALSE,verbose=FALSE,
　　32parallel.sz=1)
　　33return(es.max)
　　34})
　　35adjPvalueCutoff&lt;-0.001
　　36logFCcutoff&lt;-log2(2)
　　37es_deg&lt;-lapply(es_max,function(es.max){
　　38#es.max=es_max[[1]]
　　39table(group_list)
　　40dim(es.max)
　　41design&lt;-model.matrix(~0+factor(group_list))
　　42colnames(design)=levels(factor(group_list))
　　43rownames(design)=colnames(es.max)
　　44design
　　45library(limma)
　　46#contrast.matrix&lt;-makeContrasts(paste0(unique(group_list),collapse=＂-＂),levels=design)
　　47contrast.matrix&lt;-makeContrasts(＂npc-normal＂,
　　48levels=design)49:50contrast.matrix##这个矩阵声明，我们要把progres.组跟stable进行差异分析比较51:52deg=function(es.max,design,contrast.matrix){
　　53##step1
　　54fit&lt;-lmFit(es.max,design)
　　55##step2
　　56fit2&lt;-contrasts.fit(fit,contrast.matrix)
　　57##这一步很重要，大家可以自行看看效果58:59fit2&lt;-eBayes(fit2)##defaultnotrend!!!
　　60##eBayes()withtrend=TRUE
　　61##step3
　　62res&lt;-decideTests(fit2,p.value=adjPvalueCutoff)
　　63summary(res)
　　64tempOutput=topTable(fit2,coef=1,n=Inf)
　　65nrDEG=na.omit(tempOutput)
　　66#write.csv(nrDEG2,＂limma_notrend.results.csv＂,quote=F)
　　67head(nrDEG)
　　68return(nrDEG)
　　69}70:71re=deg(es.max,design,contrast.matrix)
　　72nrDEG=re
　　73head(nrDEG)
　　74return(nrDEG)
　　75})
　　76gmts
　　77save(es_max,es_deg,file=＂data/gsva_msigdb.Rdata＂)
　　78}
　　79}80:81#画图展示，结果存放在pic/下
　　82if(T){
　　83load(file=＂data/gsva_msigdb.Rdata＂)
　　84library(pheatmap)
　　85lapply(1:length(es_deg),function(i){
　　86#i=8
　　87print(i)
　　88dat=es_max[[i]]
　　89df=es_deg[[i]]
　　90df=df[df$P.Value&lt;0.01&gt;0.3,]<!--0.01-->
　　91print(dim(df))
　　92if(nrow(df)&gt;5){
　　93n=rownames(df)
　　94dat=dat[match(n,rownames(dat)),]
　　95ac=data.frame(g=group_list)
　　96rownames(ac)=colnames(dat)
　　97rownames(dat)=substring(rownames(dat),1,50)
　　98pheatmap::pheatmap(dat,
　　99fontsize_row=8,height=11,
　　100annotation_col=ac,show_colnames=F,
　　101filename=paste0(＂[pic/gsva_＂,strsplit(gmts[i],＂[.]＂)[[1]][1],＂.pdf＂))102:103}
　　104})105:106adjPvalueCutoff&lt;-0.001
　　107logFCcutoff&lt;-log2(2)
　　108df=do.call(rbind,es_deg)
　　109es_matrix=do.call(rbind,es_max)
　　110df=df[df$P.Value&lt;0.01&gt;0.5,]<!--0.01-->
　　111write.csv(df,file=＂data/GSVA_DEG.csv＂)
　　112}
　　因为用到的这个样本用GSVA没有得到显著性结果，所以没有图出来，具体也没有深究，有需要日后再仔细研究吧
　　step7-visualization
　　http://www.360doc.com/conteant/18/0309/18/33459258_735717104.shtml
　　1###主要是关于KEGG方面的扩展图
　　2###主要是关于KEGG方面的扩展图
　　3　　4#载入数据
　　5if(T){
　　6rm(list=ls())
　　7options(stringsAsFactors=F)
　　8　　9library(ggplot2)
　　10library(clusterProfiler)
　　11library(org.Hs.eg.db)
　　12load(file=＂data/anno_DEG.Rdata＂)13:14head(DEG)
　　15}16:17#pre-processdata
　　18if(T){
　　19gene_up=DEG[DEG$g==＂UP＂,＂ENTREZID＂]
　　20gene_down=DEG[DEG$g==＂DOWN＂,＂ENTREZID＂]
　　21gene_diff=c(gene_up,gene_down)
　　22gene_all=as.character(DEG[,＂ENTREZID＂])
　　23#制作差异基因listL
　　24boxplot(DEG$logFC)25:26geneList=DEG$logFC
　　27names(geneList)=DEG$ENTREZID
　　28geneList=sort(geneList,decreasing=T)
　　29}30:31
　　32#KEGG富集分析得到结果
　　33if(T){
　　34if(!file.exists(＂data/enrichkk.rdata＂)){
　　35gene_down
　　36gene_up
　　37enrichKK&lt;-enrichKEGG(gene=gene_up,
　　38organism=＂hsa＂,
　　39#universe=gene_all,
　　40pvalueCutoff=0.1,
　　41qvalueCutoff=0.1)
　　42save(enrichKK,file=＂data/enrichkk.rdata＂)
　　43}
　　44load(file=＂data/enrichkk.rdata＂)
　　45#查看KEGG结果
　　46head(enrichKK)[,1:6]
　　47#打开网页看相关KEGG通路图
　　48browseKEGG(enrichKK,＂hsa05150＂)49:50#将数据中的entrz-id变成symbol
　　51#更为易读
　　52enrichKK=DOSE::setReadable(enrichKK,OrgDb=＂org.Hs.eg.db＂,keyType=＂ENTREZID＂)
　　53enrichKK
　　54}55:56
　　57##可视化
　　58#条带图
　　59if(T){
　　60#par(mfrow=c(2,1))
　　61barplot(enrichKK,showCategory=20)
　　62ggsave(＂pic/barplot.png＂)
　　63}64:65#气泡图
　　66if(T){
　　67dotplot(enrichKK)
　　68ggsave(＂pic/dotplot.png＂)
　　69}70:71#下面的图需要映射颜色，设置和示例数据一样的geneList72:73#展示top5通路的共同基因，要放大看。
　　74#Gene-ConceptNetwork
　　75if(T){
　　76cnetplot(enrichKK,foldChange=geneList,colorEdge=TRUE,circular=F)
　　77ggsave(＂pic/cnetplot.png＂)
　　78cnetplot(enrichKK,foldChange=geneList,colorEdge=TRUE,circular=T)
　　79ggsave(＂pic/cnetplot_circular.png＂)
　　80}81:82
　　83#EnrichmentMap
　　84if(T){
　　85emapplot(enrichKK)
　　86ggsave(＂pic/Enrichment_Map.png＂)
　　87}88:89#(4)展示通路关系,仅仅是针对于GO数据库结果。
　　90#goplot(enrichKK)
　　91#(5)Heatmap-likefunctionalclassification
　　92if(T){
　　93heatplot(enrichKK,foldChange=geneList)
　　94ggsave(＂pic/Enrichment_Heatmap.png＂)
　　95}
　　条带图：
　　气泡图：
　　Gene-Concept Network图：每一个小蓝圈表示一个基因，其颜色表示FC值，每个KEGGterm圈的大小由里面包含基因的数目决定。
　　成环：更炫酷了，但是感觉图形展示不方便
　　不成环：信息展示更有力吧
　　Enrichment Map图：这里和上面的图类似，只不过不再显示具体的基因，而是直接画出每个term和term之间的关系，每个圆圈代表着一个term，圆圈大小代表着有多少个基因，颜色表示p值。
　　如果term和term之间有共同的基因，那么就会连接起来，聚在一起。
　　Heatmap-like functional classification：
　　和我们常规的热图不太像，这里纵轴是每个KEGG通路，横轴是涉及到的基因。颜色表示FC值。
　　致谢
　　上面所有的代码都来自生信技能树曾老板jimmy的帮助，同时我在测试运行的过程中又进行了部分改进和增补。
　　就用曾老板亲自编辑的感谢词来感谢吧：
　　Fat Yang thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes !
网站目录投稿：安儿