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

　　GEO数据挖掘流程代码版（方便抄袭）
　　微信公众号：生信小知识
　　关注可了解更多的教程及单细胞知识。问题或建议，请公众号留言；
　　内容目录
　　step0installpackagesstep1downloaddata1。背景知识2。关于GEO中的几个文件说明A。2种familyfile（s）B。SeriesMatrixFile（s）C。GSE64634RAW。tar3。关于下载镜像问题4。关于探针id转换idmap1包针对bioconductor包附加小功能对基因名字进行注释（annoGene）idmap2包万能芯片探针ID注释平台包（提取soft文件）idmap3包idmap1和idmap2都无法注释的平台AnnoProbe包5。整体代码这里抄step2checkstep3DEGstep4annogokegg真正代码detailedplot综合显示图（更推荐）step5annoGSEAstep6annoGSVAstep7visualization致谢
　　step0installpackages
　　简单安装R包和设置镜像代码：
　　1镜像设置
　　2if（T）｛
　　3rm（listls（））
　　4options（）repos
　　5options（）BioCmirror
　　6options（BioCmirrorhttps：mirrors。ustc。edu。cnbioc）
　　7options（reposc（CRANhttps：mirrors。tuna。tsinghua。edu。cnCRAN））
　　8options（）repos
　　9options（）BioCmirror
　　10｝11：12https：bioconductor。orgpackagesreleasebiochtmlGEOquery。html
　　13BiocManager安装的R包
　　14if（T）｛
　　15if（！requireNamespace（BiocManager，quietlyTRUE））
　　16install。packages（BiocManager）17：18pkgsc（KEGG。db，GSEABase，GSVA，clusterProfiler，GEOquery，limma，impute，genefu，org。Hs。eg。db，hgu133plus2。db）
　　19for（pkginpkgs）｛
　　20if（！requireNamespace（pkg，quietlyTRUE））
　　21BiocManager：：install（pkg，askF，updateF）
　　22｝
　　23｝24：25直接通过install。package函数安装的R包
　　26if（T）｛
　　27options（）repos
　　28pkgsc（FactoMineR，factoextra，WGCNA，ggplot2，pheatmap，ggpubr）
　　29for（pkginpkgs）｛
　　30if（！requireNamespace（pkg，quietlyTRUE））
　　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｝
　　step1downloaddata
　　1。背景知识
　　这里通过使用GEOqueryR包来帮助我们下载数据，比手动登录GEO数据库一个个点击要方便很多！
　　而且我发现，最近NCBI更新后，GEO数据库也更新了！
　　点击后：
　　点击后：
　　可以看到官网的代码和我们目前用的代码基本一致。
　　2。关于GEO中的几个文件说明
　　下面一个个的来说吧！
　　A。2种familyfile（s）
　　也就是SOFTformattedfamilyfile（s）和MINiMLformattedfamilyfile（s），通过字面，我们可以理解这是2种不同格式的探针说明文件。为什么我会这样说呢？因为我下载了soft文件来查看：
　　这时首行的内容：
　　就是告诉你这个数据是GPL570这个平台测序得到的：
　　中间有很多行关于GSM的，可能记录的是用到这个平台测序的sample：
　　接着就是一系列探针信息了：
　　我们可以通过在R种提取信息对我们得到的矩阵种探针做注释。
　　B。SeriesMatrixFile（s）
　　这里面包含了3部分资料：数据属性患者信息表达矩阵
　　数据属性在最前面的几行，和患者信息之间有一个空行，但是他们都是以！开头：
　　接下来是患者信息：
　　紧接着患者信息下一行会有个提示：！seriesmatrixtablebegin，然后下面就是表达矩阵信息了：
　　这个就是我们需要的表达矩阵信息了，矩阵中每一行是一个基因（也就是一个探针），每一列是一个样本（GSM）
　　C。GSE64634RAW。tar
　　这个我一开始不知道是什么东西，后来问了问jimmy，然后他给我发了个资料你要挖的公共数据集作者上传了错误的表达矩阵肿么办。大致瞅了瞅，知道这个应该是芯片的原始数据，然后也把这个文件给下载下来看了看：
　　根据常识猜了猜想：
　　X和Y应该代表着芯片上的位置，这个可能和探针有对应关系。
　　MEAN是信号的平均值，STDV是信号的标准差。
　　NPIXELS是像素点吧。
　　既然知道了这是原始数据，jimmy又给发了代码，感觉不学习下有点对不起自己，那就跟着代码过一遍吧：
　　1BiocManager：：install（c（oligo），askF，updateF）
　　2library（oligo）
　　3BiocManager：：install（c（pd。hg。u133。plus。2），askF，updateF）
　　4library（pd。hg。u133。plus。2）
　　5hr6读入cel文件数据
　　7celFilesamp；lt；list。celfiles（。，listGzippedT）
　　8celFiles
　　9affyRawamp；lt；read。celfiles（celFiles）
　　10提取矩阵并做normalization
　　11esetamp；lt；rma（affyRaw）
　　12eset
　　13检查数据
　　14datexprs（eset）
　　15dim（dat）
　　16boxplot（dat，las2）
　　可以看到，整个过程非常的简单，就是利用了oligo这个R包而已。读入所有的cel文件后，利用rma（）函数将数据进行了normalization，得到的是一个ExpressionSet对象！！
　　后面会比较我们利用rawdata得到的结果和我们直接下载的SeriesMatrixFile（s）文件之间的区别！
　　3。关于下载镜像问题
　　jimmy曾经发表过GEO数据库中国区镜像横空出世推文，因为我每次下载都是用vpn，但是怕以后万一没有vpn了，所以还是学习下，以防万一嘛！再说了，可以学习下新知识，何乐而不为呢？
　　通过学习，发现实际上jimmy是开发了一个R包，或者说包装了一个函数吧，如果不想了解具体原理，那么用法如下：
　　1安装方法1：
　　2library（devtools）
　　3installgithub（jmzeng1314GEOmirror）
　　4library（GEOmirror）
　　5安装方法2：
　　6source（http：raw。githubusercontent。comjmzeng1314GEOmirrormasterRgeoChina。R）
　　7安装方法3（源码）：
　　8geoChinaamp；lt；function（gseGSE2546，mirrortercent）｛
　　9suppressPackageStartupMessages（library（GEOquery））
　　10options（download。file。methodlibcurl）
　　11eSetgetGEO（GSE2546，destdir。，AnnotGPLF，getGPLF）
　　12http：49。235。27。111GEOmirrorGSE2nnnGSE2546eSet。Rdata
　　13gseGSE2546；mirrortercent
　　14gsetoupper（gse）
　　15downifelse（as。numeric（gsub（GSE，，gse））amp；lt；1000，
　　16paste0（GEOmirrorGSEnnn，gse，
　　17eSet。Rdata），
　　18paste0（GEOmirror，
　　19gsub（〔09〕〔09〕〔09〕，nnn，gse），，gse，
　　20eSet。Rdata））21：22if（mirrortercent）｛
　　23uphttp：49。235。27。111
　　24｝
　　25tpfpaste0（gse，eSet。Rdata）
　　26download。file（paste0（up，down），tpf，modewb）
　　27suppressWarnings（load（tpf））
　　28return（gset）
　　29｝
　　具体用法也非常简单：
　　1eSetgeoChina（GSE2345）
　　通过安装方法3（源码）我们可以看出这个包的原理：
　　通过访问下载http：49。235。27。111GEOmirrorGSE2nnnGSE2546eSet。Rdata
　　所以可以猜到，应该是jimmy事先用循环的方式帮我们下载好了很多GEO数据，并做成了Rdata格式的文件！非常的良苦用心了！！
　　4。关于探针id转换
　　因为GEO数据矩阵中，横坐标都是探针的代号，如下图：
　　我们只看这些代号并不能理解具体是什么基因，于是这就需要我们做id转换了：将探针的代号转为基因symbol。
　　idmap1包针对bioconductor包
　　这里又要提到jimmy开发的一个R包了：
　　第一个万能芯片探针ID注释平台R包37个芯片平台
　　老规矩，还是先学习下：
　　一般重要的芯片在R的bioconductor里面都是有包的。但是需要我们单独下载对应的包。具体关系如下：
　　（具体的R包名称是biocpackage后加。db）
　　1gplbiocpackagetitle
　　2GPL201hgfocus〔HGFocus〕AffymetrixHumanHGFocusTargetArray
　　3GPL96hgu133a〔HGU133A〕AffymetrixHumanGenomeU133AArray
　　4GPL571hgu133a2〔HGU133A2〕AffymetrixHumanGenomeU133A2。0Array
　　5GPL97hgu133b〔HGU133B〕AffymetrixHumanGenomeU133BArray
　　6GPL570hgu133plus2〔HGU133Plus2〕AffymetrixHumanGenomeU133Plus2。0Array
　　7GPL13667hgu219〔HGU219〕AffymetrixHumanGenomeU219Array
　　8GPL8300hgu95av2〔HGU95Av2〕AffymetrixHumanGenomeU95Version2Array
　　9GPL91hgu95av2〔HGU95A〕AffymetrixHumanGenomeU95AArray
　　10GPL92hgu95b〔HGU95B〕AffymetrixHumanGenomeU95BArray
　　11GPL93hgu95c〔HGU95C〕AffymetrixHumanGenomeU95CArray
　　12GPL94hgu95d〔HGU95D〕AffymetrixHumanGenomeU95DArray
　　13GPL95hgu95e〔HGU95E〕AffymetrixHumanGenomeU95EArray
　　14GPL887hgug4110bAgilent012097Human1AMicroarray（V2）G4110B（FeatureNumberversion）
　　15GPL886hgug4111aAgilent011871Human1BMicroarrayG4111A（FeatureNumberversion）
　　16GPL1708hgug4112aAgilent012391WholeHumanGenomeOligoMicroarrayG4112A（FeatureNumberversion）
　　17GPL13497HsAgilentDesign026652Agilent026652WholeHumanGenomeMicroarray4x44Kv2（ProbeNameversion）
　　18GPL6244hugene10sttranscriptcluster〔HuGene10st〕AffymetrixHumanGene1。0STArray〔transcript（gene）version〕
　　19GPL11532hugene11sttranscriptcluster〔HuGene11st〕AffymetrixHumanGene1。1STArray〔transcript（gene）version〕
　　20GPL6097illuminaHumanv1Illuminahuman6v1。0expressionbeadchip
　　21GPL6102illuminaHumanv2Illuminahuman6v2。0expressionbeadchip
　　22GPL6947illuminaHumanv3IlluminaHumanHT12V3。0expressionbeadchip
　　23GPL10558illuminaHumanv4IlluminaHumanHT12V4。0expressionbeadchip
　　24GPL6885illuminaMousev2IlluminaMouseRef8v2。0expressionbeadchip
　　25GPL81mgu74av2〔MGU74Av2〕AffymetrixMurineGenomeU74AVersion2Array
　　26GPL82mgu74bv2〔MGU74Bv2〕AffymetrixMurineGenomeU74BVersion2Array
　　27GPL83mgu74cv2〔MGU74Cv2〕AffymetrixMurineGenomeU74Version2Array
　　28GPL339moe430a〔MOE430A〕AffymetrixMouseExpression430AArray
　　29GPL6246mogene10sttranscriptcluster〔MoGene10st〕AffymetrixMouseGene1。0STArray〔transcript（gene）version〕
　　30GPL340mouse4302〔MOE430B〕AffymetrixMouseExpression430BArray
　　31GPL1261mouse430a2〔Mouse4302〕AffymetrixMouseGenome4302。0Array
　　32GPL8321mouse430a2〔Mouse430A2〕AffymetrixMouseGenome430A2。0Array
　　因为很多时候，用户是找不到这些GPL平台对应的R包，或者下载安装困难，其实仅仅是需要探针与基因对应关系，没有必要去下载安装这几十个M的包。于是jimmy就开发了idmap1这个R包来帮我们：
　　安装方法（这里好像只能用方法1，因为在载入包的同时附带有一个变量p2sdf加载，如果用方法2和3没办法得到这个变量）：
　　也可以直接下载：https：codeload。github。comjmzeng1314idmap1legacy。tar。gzmaster
　　1安装方法1：
　　2library（devtools）
　　3installgithub（jmzeng1314idmap1）
　　4library（idmap1）
　　5安装方法2：
　　6source（https：raw。githubusercontent。comjmzeng1314idmap1masterRgetIDs。R）
　　7安装方法3（源码）：
　　8getIDsamp；lt；function（gpl）｛
　　9gplgpl570
　　10gpltoupper（gpl）
　　11if（！gplinunique（p2sdfgpl））｛
　　12stop（yourgplisnotinourannotationlist。）
　　13｝
　　14return（p2sdf〔p2sdfgplgpl，〕）
　　15｝
　　具体用法也非常简单：
　　1idsgetIDs（gpl570）
　　2head（ids）
　　关于我们得到的ExpressionSet对象，可以通过gsetannotation得到我们需要的注释平台信息。
　　前面说到了这个变量p2sdf，像我这么喜欢资源的人，当然也要保存一份在本地了。大家有兴趣的可以自己write到本地保存。哈哈哈哈哈哈哈哈哈哈（jimmy应该不会打我吧）！
　　附加小功能对基因名字进行注释（annoGene）
　　同样的，还是idmap1这个R包里的函数，如果安装过这个R包就无须再安装了，如果没有安装又想用这个功能，还真的没有办法，因为这些数据存在于这个包的自带变量中（humanGTF、mouseGTF、ratGTF）：
　　1安装方法2（无效）：
　　2source（https：raw。githubusercontent。comjmzeng1314idmap1masterRannoGene。R）
　　3安装方法3（源码）（无效）：
　　4annoGeneamp；lt；function（IDs，IDtype，specieshuman，outfile）｛
　　5if（length（unique（IDs））amp；lt；1）｛
　　6stop（Youshouldgivemesomegenestobeannotated！！！）
　　7｝
　　8if（！IDtypeinc（ENSEMBL，SYMBOL））｛
　　9stop（WeonlyacceptENSEMBLorSYMBOL！！！）
　　10｝
　　11if（specieshuman）｛
　　12GTFamp；lt；humanGTF
　　13｝elseif（speciesmouse）｛
　　14GTFamp；lt；mouseGTF
　　15｝elseif（speciesrat）｛
　　16GTFamp；lt；ratGTF
　　17｝else｛
　　18stop（Weonlyaccepthumanormouse，orrat，）
　　19｝
　　20resamp；lt；GTF〔eval（parse（textpaste0（GTF，IDtype）））inIDs，〕21：22missIdsamp；lt；IDs〔！（IDsineval（parse（textpaste0（res，IDtype））））〕
　　23missIdsPercentageround（（length（missIds）length（IDs））100，2）
　　24if（length（missIds）！0）｛
　　25warning（
　　26paste0（missIdsPercentage，ofinputIDsarefailtoannotate。。。）
　　275。29ofinputgeneIDsarefailtomap。。。
　　28）
　　29｝
　　30if（hasArg（outfile））｛
　　31resultsres
　　32if（grepl（。html，outfile））｛
　　33Ensemblprefixamp；lt；https：asia。ensembl。orgHomosapiensGeneSummary？g
　　34hrefpaste0（Ensemblprefix，resultsENSEMBL）
　　35resultsENSEMBLpaste0（lt；stronggt；lt；atargetblackhref，gt；，resultsENSEMBL，lt；agt；lt；stronggt；）36：37symbolprefixamp；lt；http：www。ncbi。nlm。nih。govgene？term
　　38hrefpaste0（symbolprefix，resultsSYMBOL）
　　39resultsSYMBOLpaste0（lt；stronggt；lt；atargetblackhref，gt；，resultsSYMBOL，lt；agt；lt；stronggt；）40：41yamp；lt；DT：：datatable（results，escapeF，rownamesF）
　　42DT：：saveWidget（y，fileoutfile）
　　43｝elseif（grepl（。csv，outfile））｛
　　44write。csv（results，fileoutfile）
　　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种familyfile（s））进行提取信息得到的！
　　如果是我自己来处理这样的文件，我应该会分2步：
　　用linux中的grepv命令将表头以和！开头的行去掉
　　1grepv！GPL570。softgrepvgrepvamp；gt；GPL570。txt
　　用R读入数据
　　1tmpread。table（dataGPL570。txt，headerT，sep，fillT）
　　结果证明这样有错误，但是，具体原因有空再去找吧。下面看正确的读入方法借助GEOqueryR包工具：
　　1library（GEOquery）
　　2gplamp；lt；getGEO（GPL570，destdirdata）
　　3GPLTable（gpl）
　　一般来说，大家关心的其实就是探针的ID，以及基因的symbol列。有了这个变量后，就可以按照R语言基本操作来提取我们需要的信息了。
　　注意：我检查了得到的结果，里面存在有的探针ID对应2个基因，如下图：
　　虽然不知道这些代表着什么意思，但是，我将这个数据和bioconductor包里的hgu133plus2。db数据做了比较，结果是这样的：自己提取的结果中如果是一个ID对应2个基因，那么这个探针在bioconductor包里基本上找不到数据。而其他一个ID对应2个基因的结果均和bioconductor包一致。
　　当然了，我这只是单独来看一个平台的探针，而在idmap2jimmy已经帮我们整理好了，直接用就行了！！
　　安装方法：
　　（比较慢）
　　也可以直接下载：https：codeload。github。comjmzeng1314idmap2legacy。tar。gzmaster
　　1library（devtools）
　　2installgithub（jmzeng1314idmap2）
　　3library（idmap2）
　　使用方法：
　　1library（idmap2）
　　2idsgetsoftIDs（gpl570）
　　查看支持的平台：
　　1gpllist〔，1：4〕gpllist是R包自带的变量
　　idmap3包idmap1和idmap2都无法注释的平台
　　第三个万能芯片探针ID注释平台R包idmap1和idmap2都无法注释的平台
　　如果我们拿到的soft注释文件中是序列信息，那么我们该怎么做呢？
　　应该是先将序列比对到参考基因组上，然后通过提取基因注释文件中的数据得到基因symbol！
　　img
　　而在idmap3包中，jimmy已经帮我们做好了！！说他宠粉也是真的了！！我都懒得做，他居然还写了个循环来完成了这个事。
　　安装方法：
　　（比较慢）
　　也可以直接下载：https：codeload。github。comjmzeng1314idmap3legacy。tar。gzmaster
　　1library（devtools）
　　2installgithub（jmzeng1314idmap3）
　　3library（idmap3）
　　使用方法：
　　1library（idmap3）
　　2idsidmap3：：getpipeIDs（GPL21827）
　　查看支持的平台：
　　1gpllist〔，1：4〕gpllist是R包自带的变量
　　AnnoProbe包
　　芯片探针序列的基因注释已经无需你自己亲自做了一个汇总
　　安装方法：
　　https：codeload。github。comjmzeng1314AnnoProbelegacy。tar。gzmaster
　　1library（devtools）
　　2installgithub（jmzeng1314AnnoProbe）
　　3library（AnnoProbe）
　　使用方法：
　　1gplGPL21827
　　2probe2geneidmap（gpl，typepipe）
　　5。整体代码这里抄
　　1rm（listls（））魔幻操作，一键清空
　　2options（stringsAsFactorsF）
　　3下面代码只需要修改file的名字
　　4下载的文件会自动保存到。data下
　　5得到的gset是一个ExpressionSet对象
　　6hr7hr8只需要修改file的名字，即可下载得到相应的geo数据
　　9获取患者信息，需要自行修改
　　10fileGSE64634
　　11https：www。ncbi。nlm。nih。govgeoqueryacc。cgi？accGSE6463412：13if（T）｛
　　14数据下载
　　15if（T）｛
　　16library（GEOquery）
　　17这个包需要注意两个配置，一般来说自动化的配置是足够的。
　　18Settingoptions（download。file。method。GEOqueryauto）
　　19Settingoptions（GEOquery。inmemory。gplFALSE）
　　20fdatapaste0（file，eSet。Rdata）
　　21fpathpaste0（data，fdata）
　　22if（！file。exists（fpath））｛
　　23gsetamp；lt；getGEO（file，destdirdata，
　　24AnnotGPLF，注释文件
　　25getGPLF）平台文件
　　26gsetgset〔〔1〕〕
　　27save（gset，filefpath）保存到本地
　　28｝
　　29load（fpath）
　　30｝
　　31gset32：33获取患者信息，这里需要自行修改
　　34if（T）｛
　　35pdpData（gset）根据diseasestate：ch1一列得知分组信息
　　36grouplistc（rep（normal，4），rep（npc，12））nasopharyngealcarcinomaNPC鼻咽癌
　　37table（grouplist）
　　38｝39：40对数据进行normalization
　　41if（T）｛
　　42datexprs（gset）
　　43dim（dat）44：45dat〔1：4，1：4〕
　　46boxplot（dat，las2）
　　47datdat〔apply（dat，1，sd）amp；gt；0，〕去除都是0的探针
　　48dat〔datamp；lt；0〕1
　　49boxplot（dat，las2）50：51datlog2（dat1）
　　52boxplot（dat，las2）
　　53library（limma）
　　54datnormalizeBetweenArrays（dat）
　　55boxplot（dat，las2）
　　56｝57：58
　　59探针注释2种方法，推荐方法2
　　60方法1：比较麻烦，而且不方便，一般不用这种方法
　　61if（F）｛
　　62library（GEOquery）
　　63DownloadGPLfile，putitinthecurrentdirectory，andloadit：
　　64gplamp；lt；getGEO（GPL570，destdirdata）
　　65GPLTable（gpl）
　　66colnames（GPL）
　　67head（GPL〔，c（1，11）〕）youneedtocheckthis，whichcolumndoyouneed
　　68probe2geneGPL〔，c（1，11）〕
　　69head（probe2gene）
　　70save（probe2gene，fileprobe2gene。Rdata）
　　71｝72：73方法2：用hgu133plus2。db这个R包比较方便
　　74if（T）｛
　　75library（hgu133plus2。db）
　　76idstoTable（hgu133plus2SYMBOL）
　　77head（ids）
　　78idsids〔idssymbol！，〕
　　79idsids〔idsprobeidinrownames（dat），〕过滤没法注释的探针80：81dat〔1：4，1：4〕
　　82datdat〔idsprobeid，〕调整顺序，让dat的顺序和ids中的一致83：84idsmedianapply（dat，1，median）
　　85idsids〔order（idssymbol，idsmedian，decreasingT），〕按照基因名、中位数大小排序
　　86idsids〔！duplicated（idssymbol），〕只保留相同symbol中中位数最大的探针
　　87datdat〔idsprobeid，〕调整顺序，让dat的顺序和ids中的一致
　　88rownames（dat）idssymbolid转换
　　89dat〔1：4，1：4〕
　　90｝
　　91｝92：93save（dat，grouplist，filedatastep1output。Rdata）
　　step2check
　　画PCA图时要求是行名时样本名，列名时探针名，因此此时需要转置。格式要求data。frame。
　　关于scale：在基因表达量的数据中，有些基因表达量极低，有些基因表达量极高，因此把每个基因在不同处理和重复中的数据转换为平均值为0，方差为1的数据，所以在这里也需要先转置。
　　1PCA检查
　　2PCA，图会保存在picallsamplesPCA。png
　　3if（T）｛
　　4载入数据
　　5if（T）｛
　　6rm（listls（））魔幻操作，一键清空
　　7options（stringsAsFactorsF）8：9load（filedatastep1output。Rdata）
　　10print（table（grouplist））
　　11每次都要检测数据
　　12dat〔1：4，1：4〕
　　13｝14：15PCA，图会保存在picallsamplesPCA。png
　　16if（T）｛
　　17exprSetdat
　　18过滤标准
　　19if（T）｛
　　20dim（exprSet）
　　21过滤标准需要修改，目前是保留每一个基因在amp；gt；5个人中表达量amp；gt；1
　　22exprSetexprSet〔apply（exprSet，1，function（x）sum（xamp；gt；1）amp；gt；5），〕
　　23boxplot（apply（exprSet，1，mean））
　　24dim（exprSet）
　　25过滤标准需要修改，目前是去除每一个基因在amp；gt；5个人中表达量amp；gt；12
　　26exprSetexprSet〔！apply（exprSet，1，function（x）sum（xamp；gt；12）amp；gt；5），〕
　　27dim（exprSet）
　　28｝
　　29去除文库大小差异normalization，同时利用log做scale
　　30exprSetlog（edgeR：：cpm（exprSet）1）
　　31datexprSet
　　32datas。data。frame（t（dat））画PCA图时要求是行名时样本名，列名时探针名，因此此时需要转换。格式要求data。frame
　　33library（FactoMineR）计算PCA
　　34library（factoextra）画图展示35：36dat。pcaamp；lt；PCA（dat，graphF）
　　37fvizpcaind按样本fvizpcavar按基因
　　38fvizpcaind（dat。pca，
　　39geom。indpoint，c（point，text）2选1
　　40col。indgrouplist，colorbygroups
　　41palettec（00AFBB，E7B800），自定义颜色
　　42addEllipsesT，加圆圈
　　43legend。titleGroups图例名称
　　44）
　　45ggsave（picallsamplesPCA。png）
　　46｝
　　47｝48：49
　　50挑选1000个SD最大的基因画表达量热图
　　51结果图存放在picheatmaptop1000sd。png中
　　52if（T）｛
　　53载入数据
　　54if（T）｛
　　55rm（listls（））
　　56load（filedatastep1output。Rdata）
　　57dat〔1：4，1：4〕
　　58｝59：60挑选1000个SD最大的基因画表达量热图
　　61结果图存放在picheatmaptop1000sd。png中
　　62if（T）｛
　　63cgnames（tail（sort（apply（dat，1，sd）），1000））找到SD最大的1000个基因
　　64library（pheatmap）
　　65headmap。datdat〔cg，〕
　　66把每个基因在不同处理和重复中的数据转换为平均值为0，
　　67方差为1的数据，所以在这里也需要先转置（针对基因）。
　　68headmap。datt（scale（t（headmap。dat）））
　　69headmap。dat〔headmap。datamp；gt；2〕2
　　70headmap。dat〔headmap。datamp；lt；2〕271：72分组信息设置
　　73acdata。frame（groupgrouplist）
　　74rownames（ac）colnames（headmap。dat）75：76可以看到TNBC具有一定的异质性，拿它来区分乳腺癌亚型指导临床治疗还是略显粗糙。
　　77pheatmap（headmap。dat，showcolnamesT，showrownamesF，
　　78annotationcolac，
　　79filenamepicheatmaptop1000sd。png）
　　80｝
　　81｝82：83
　　84过滤标准需要修改
　　85利用top500mad基因画相关性热图热图
　　86结果图存放在piccortop500mad。png中
　　87if（T）｛
　　88导入数据
　　89if（T）｛
　　90rm（listls（））
　　91load（filedatastep1output。Rdata）
　　92dat〔1：4，1：4〕
　　93exprSetdat
　　94｝95：96利用所有基因画相关性热图
　　97if（T）｛
　　98acdata。frame（groupgrouplist）
　　99rownames（ac）colnames（exprSet）
　　100pheatmap：：pheatmap（cor（exprSet），
　　101annotationcolac，
　　102showrownamesF，
　　103filenamepiccorallgenes。png）
　　104｝105：106过滤标准
　　107if（T）｛
　　108dim（exprSet）
　　109过滤标准需要修改，目前是保留每一个基因在amp；gt；5个人中表达量amp；gt；1
　　110exprSetexprSet〔apply（exprSet，1，function（x）sum（xamp；gt；1）amp；gt；5），〕
　　111boxplot（apply（exprSet，1，mean））
　　112dim（exprSet）
　　113过滤标准需要修改，目前是去除每一个基因在amp；gt；5个人中表达量amp；gt；12
　　114exprSetexprSet〔！apply（exprSet，1，function（x）sum（xamp；gt；12）amp；gt；5），〕
　　115dim（exprSet）
　　116｝117：118数据normalization处理
　　119if（T）｛
　　120去除文库大小差异normalization，同时利用log做scale
　　121exprSetlog（edgeR：：cpm（exprSet）1）
　　122取top500mad的基因画图
　　123exprSetexprSet〔names（sort（apply（exprSet，1，mad），decreasingT）〔1：500〕），〕
　　124｝125：126利用top500mad基因画相关性热图热图
　　127if（T）｛
　　128Mcor（exprSet）
　　129pheatmap：：pheatmap（M，
　　130showrownamesF，
　　131annotationcolac，
　　132filenamepiccortop500mad。png）
　　133｝
　　134｝
　　allsamplesPCA。png
　　heatmaptop1000sd。png
　　cortop500mad。png
　　step3DEG
　　关于差异分析是否需要比较矩阵
　　差异分析的统计学本质探究
　　原始版火山图：plot（logFC，log10（P。Value））
　　原始版MA图：plot（AveExpr，logFC）
　　1载入数据，检查数据
　　2if（T）｛
　　3rm（listls（））
　　4options（stringsAsFactorsF）
　　5library（ggpubr）
　　6load（filedatastep1output。Rdata）
　　7每次都要检测数据
　　8dat〔1：4，1：4〕
　　9table（grouplist）
　　10boxplot（dat〔1，〕grouplist）按照grouplist分组画箱线图
　　11hr12boxplot的美化版
　　13bplotfunction（g）｛
　　14dfdata。frame（geneg，stagegrouplist）
　　15pamp；lt；ggboxplot（df，xstage，ygene，
　　16colorstage，palettejco，
　　17addjitter）
　　18Addpvalue
　　19pstatcomparemeans（）
　　20｝
　　21｝
　　22利用定义好的函数检查数据
　　23bplot（dat〔1，〕）
　　24bplot（dat〔2，〕）
　　25bplot（dat〔3，〕）
　　26bplot（dat〔4，〕）
　　27hr28hr29limma
　　30library（limma）
　　31方法1：不制作比较矩阵，简单
　　32但是做不到随心所欲的指定任意两组进行比较
　　33if（T）｛
　　34designmodel。matrix（factor（grouplist））
　　35fitlmFit（dat，design）
　　36fit2eBayes（fit）
　　37上面是limma包用法的一种方式
　　38options（digits4）设置全局的数字有效位数为4
　　39topTable（fit2，coef2，adjustBH）
　　40｝
　　41hr42方法2：制作比较矩阵
　　43可以随心所欲的指定任意两组进行比较
　　44if（T）｛
　　45designamp；lt；model。matrix（0factor（grouplist））
　　46colnames（design）levels（factor（grouplist））
　　47head（design）
　　48exprSetdat
　　49rownames（design）colnames（exprSet）
　　50design
　　51hr52比较矩阵
　　53这个矩阵声明，我们要把npc组跟Normal进行差异分析比较
　　54contrast。matrixamp；lt；makeContrasts（npcnormal，
　　55levelsdesign）
　　56contrast。matrix
　　57hr58degfunction（exprSet，design，contrast。matrix）｛
　　59step1
　　60fitamp；lt；lmFit（exprSet，design）
　　61step2
　　62fit2amp；lt；contrasts。fit（fit，contrast。matrix）
　　63hr64fit2amp；lt；eBayes（fit2）defaultnotrend！！！
　　65eBayes（）withtrendTRUE
　　66hr67step3
　　68有了比较矩阵后，coef1，而numberInf是把所有结果都打印出来
　　69tempOutputtopTable（fit2，coef1，numberInf）
　　70nrDEGna。omit（tempOutput）
　　71write。csv（nrDEG2，limmanotrend。results。csv，quoteF）
　　72head（nrDEG）
　　73return（nrDEG）
　　74｝
　　75hr76degdeg（exprSet，design，contrast。matrix）
　　77hr78head（deg）
　　79｝
　　80hr81hr82save（deg，filedatadeg。Rdata）
　　83hr84load（filedatadeg。Rdata）
　　85head（deg）
　　86bplot（dat〔rownames（deg）〔1〕，〕）
　　87forvolcanoandMAplot
　　88结果存放在picvolcano。png和picMA。png
　　89if（T）｛
　　90nrDEGdeg
　　91head（nrDEG）
　　92attach（nrDEG）
　　93原始版火山图
　　94plot（logFC，log10（P。Value））
　　95library（ggpubr）
　　96dfnrDEG
　　97dfylog10（P。Value）
　　98ggscatter（df，xlogFC，yy，size0。5）
　　99定义logFC2为阈值
　　100dfstateifelse（dfP。Valueamp；gt；0。01，stable，
　　101ifelse（dflogFCamp；gt；2，up，
　　102ifelse（dflogFCamp；lt；2，down，stable））
　　103）
　　104table（dfstate）
　　105dfnamerownames（df）
　　106head（df）
　　107ggscatter（df，xlogFC，yy，size0。5，colorstate）
　　108ggscatter（df，xlogFC，yy，colorstate，size0。5，
　　109labelname，repelT，
　　110label。selectrownames（df）〔dfstate！stable〕，
　　111label。selectc（TTC9，AQP3，CXCL11，PTGS2），挑选一些基因在图中显示出来
　　112palettec（00AFBB，E7B800，FC4E07））
　　113ggsave（picvolcano。png）114：115MA图
　　116ggscatter（df，xAveExpr，ylogFC，size0。2）
　　117dfpcifelse（dfP。Valueamp；lt；0。001，pamp；lt；0。001，
　　118ifelse（dfP。Valueamp；lt；0。01，0。001amp；lt；pamp；lt；0。01，pamp；gt；0。01））lt；！0。01，pgt；
　　119table（dfpc）
　　120ggscatter（df，xAveExpr，ylogFC，colorpc，size0。2，
　　121palettec（green，red，black））
　　122ggsave（picMA。png）
　　123｝124：125forheatmap
　　126if（T）｛
　　127load（filedatastep1output。Rdata）
　　128每次都要检测数据
　　129dat〔1：4，1：4〕
　　130table（grouplist）
　　131xdeglogFC
　　132names（x）rownames（deg）133：134cg中存放着变化上升和下降的前100个基因名
　　135cgc（names（head（sort（x），100）），
　　136names（tail（sort（x），100）））
　　137library（pheatmap）
　　138nt（scale（t（dat〔cg，〕）））139：140n〔namp；gt；2〕2
　　141n〔namp；lt；2〕2
　　142n〔1：4，1：4〕
　　143pheatmap（n，showcolnamesF，showrownamesF）
　　144acdata。frame（groupgrouplist）
　　145rownames（ac）colnames（n）将ac的行名也就分组信息（是‘noTNBC’还是‘TNBC’）给到n的列名，即热图中位于上方的分组信息
　　146pheatmap（n，showcolnamesF，
　　147showrownamesF，
　　148clustercolsF，
　　149annotationcolac，filenamepicheatmaptop200DEG。png）列名注释信息为ac即分组信息150：151
　　152｝153：154write。csv（deg，filedatadeg。csv）
　　volcano。png
　　MA。png
　　heatmaptop200DEG。png
　　step4annogokegg
　　关于基因名和ID之间的各种转换：bitr（geneID，fromType，toType，OrgDb，dropTRUE）例如：
　　1head（unique（degsymbol））
　　2〔1〕LOC388780SLC6A6RGS22AKR1D1AOC1FOXJ1
　　3bitr（unique（degsymbol），fromTypeSYMBOL，
　　4toTypec（ENTREZID），
　　5OrgDborg。Hs。eg。db）
　　这里将KEGG和GO富集分析函数自定义为3个函数runkegg、rungo、keggplot
　　每次运行前先运行下面代码：
　　1KEGGpathwayanalysis
　　2具体结果数据在data下，图在pic下
　　3runkeggamp；lt；function（geneup，genedown，geneListF，protest）｛
　　4geneupunique（geneup）
　　5genedownunique（genedown）
　　6genediffunique（c（geneup，genedown））
　　7overrepresentationtest
　　8下面把3个基因集分开做超几何分布检验
　　9首先是上调基因集。
　　10kk。upamp；lt；enrichKEGG（genegeneup，
　　11organismhsa，
　　12universegeneall，
　　13pvalueCutoff0。9，
　　14qvalueCutoff0。9）
　　15head（kk。up）〔，1：6〕
　　16kkkk。up
　　17dotplot（kk）
　　18barplot（kk）
　　19kkDOSE：：setReadable（kk，OrgDborg。Hs。eg。db，keyTypeENTREZID）
　　20head（kk）
　　21write。csv（kkresult，paste0（data，pro，kk。up。csv））
　　22hr23首先是下调基因集。
　　24kk。downamp；lt；enrichKEGG（genegenedown，
　　25organismhsa，
　　26universegeneall，
　　27pvalueCutoff0。9，
　　28qvalueCutoff0。9）
　　29head（kk。down）〔，1：6〕
　　30kkkk。down
　　31dotplot（kk）
　　32barplot（kk）
　　33kkDOSE：：setReadable（kk，OrgDborg。Hs。eg。db，keyTypeENTREZID）
　　34write。csv（kkresult，paste0（data，pro，kk。down。csv））
　　35hr36最后是上下调合并后的基因集。
　　37kk。diffamp；lt；enrichKEGG（genegenediff，
　　38organismhsa，
　　39pvalueCutoff0。05）
　　40head（kk。diff）〔，1：6〕
　　41kkkk。diff
　　42dotplot（kk）
　　43barplot（kk）
　　44kkDOSE：：setReadable（kk，OrgDborg。Hs。eg。db，keyTypeENTREZID）
　　45write。csv（kkresult，paste0（data，pro，kk。diff。csv））
　　46hr47hr48keggdowndtamp；lt；as。data。frame（kk。down）
　　49keggupdtamp；lt；as。data。frame（kk。up）
　　50downkeggamp；lt；keggdowndt〔keggdowndtpvalueamp；lt；0。01，〕；downkegggroup1
　　51upkeggamp；lt；keggupdt〔keggupdtpvalueamp；lt；0。01，〕；upkegggroup1
　　52画图设置，这个图很丑，大家可以自行修改。
　　53keggplot（upkegg，downkegg）
　　54hr55ggsave（gkegg，filenamepaste0（pic，pro，keggupdown。png））
　　56hr57GSEA
　　58if（geneList）｛
　　59hr60GSEA算法跟上面的使用差异基因集做超几何分布检验不一样。
　　61kkgseamp；lt；gseKEGG（geneListgeneList，
　　62organismhsa，
　　63nPerm1000，
　　64minGSSize20，
　　65pvalueCutoff0。9，
　　66verboseFALSE）
　　67head（kkgse）〔，1：6〕
　　68gseaplot（kkgse，geneSetIDrownames（kkgse〔1，〕））
　　69gseaplot（kkgse，hsa04110，titleCellcycle）
　　70kkDOSE：：setReadable（kkgse，OrgDborg。Hs。eg。db，keyTypeENTREZID）
　　71tmpkkresult
　　72write。csv（kkresult，paste0（pro，kegg。gsea。csv））
　　73hr74hr75这里找不到显著下调的通路，可以选择调整阈值，或者其它。
　　76downkeggamp；lt；kkgse〔kkgsepvalueamp；lt；0。05amp；amp；kkgseenrichmentScoreamp；lt；0，〕；downkegggroup1
　　77upkeggamp；lt；kkgse〔kkgsepvalueamp；lt；0。05enrichmentscoreamp；gt；0，〕；upkegggroup1lt；！0。05gt；
　　78hr79gkeggkeggplot（upkegg，downkegg）
　　80print（gkegg）
　　81ggsave（gkegg，filenamepaste0（pro，kegggsea。png））
　　82｝
　　83｝
　　84hr85GOdatabaseanalysis
　　86具体结果数据在data下，图在pic下
　　87rungoamp；lt；function（geneup，genedown，protest）｛
　　88geneupunique（geneup）
　　89genedownunique（genedown）
　　90genediffunique（c（geneup，genedown））
　　91glistlist（geneupgeneup，
　　92genedowngenedown，
　　93genediffgenediff）
　　94因为go数据库非常多基因集，所以运行速度会很慢。
　　95if（T）｛
　　96goenrichresultsamp；lt；lapply（glist，function（gene）｛
　　97lapply（c（BP，MF，CC），function（ont）｛
　　98cat（paste（Nowprocess，names（gene），ont））
　　99egoamp；lt；enrichGO（genegene，
　　100universegeneall，
　　101OrgDborg。Hs。eg。db，
　　102ontont，
　　103pAdjustMethodBH，
　　104pvalueCutoff0。99，
　　105qvalueCutoff0。99，
　　106readableTRUE）107：108print（head（ego））
　　109return（ego）
　　110｝）
　　111｝）
　　112save（goenrichresults，filepaste0（data，pro，goenrichresults。Rdata））113：114｝
　　115只有第一次需要运行，就保存结果哈，下次需要探索结果，就载入即可。
　　116load（filepaste0（data，pro，goenrichresults。Rdata））117：118n1c（geneup，genedown，genediff）
　　119n2c（BP，MF，CC）
　　120for（iin1：3）｛
　　121for（jin1：3）｛
　　122fnpaste0（pic，pro，dotplot，n1〔i〕，，n2〔j〕，。png）
　　123cat（paste0（fn，））
　　124png（fn，res150，width1080）
　　125print（dotplot（goenrichresults〔〔i〕〕〔〔j〕〕））
　　126dev。off（）
　　127｝
　　128｝
　　129｝130：131132：133合并up和down的基因KEGG结果，放在一幅图里展示
　　134keggplotamp；lt；function（upkegg，downkegg）｛
　　135datrbind（upkegg，downkegg）
　　136colnames（dat）
　　137datpvaluelog10（datpvalue）
　　138datpvaluedatpvaluedatgroup139：140datdat〔order（datpvalue，decreasingF），〕141：142gkeggamp；lt；ggplot（dat，aes（xreorder（Description，order（pvalue，decreasingF）），ypvalue，fillgroup））
　　143geombar（statidentity）
　　144scalefillgradient（lowblue，highred，guideFALSE）
　　145scalexdiscrete（namePathwaynames）
　　146scaleycontinuous（namelog10Pvalue）
　　147coordflip（）themebw（）theme（plot。titleelementtext（hjust0。5））
　　148ggtitle（PathwayEnrichment）
　　149print（gkegg）
　　150｝
　　真正代码
　　1载入数据
　　2if（T）｛
　　3rm（listls（））
　　4options（stringsAsFactorsF）5：6load（filedatadeg。Rdata）
　　7head（deg）
　　8｝9：10数据预处理
　　11不同的阈值，筛选到的差异基因数量就不一样，后面的超几何分布检验结果就大相径庭。
　　12logFCt1。5
　　13预处理1
　　14if（T）｛
　　15deggifelse（degP。Valueamp；gt；0。05，stable，
　　16ifelse（deglogFCamp；gt；logFCt，UP，
　　17ifelse（deglogFCamp；lt；logFCt，DOWN，stable））
　　18）
　　19table（degg）
　　20head（deg）
　　21degsymbolrownames（deg）
　　22library（ggplot2）
　　23library（clusterProfiler）
　　24library（org。Hs。eg。db）
　　25dfamp；lt；bitr（unique（degsymbol），fromTypeSYMBOL，
　　26toTypec（ENTREZID），
　　27OrgDborg。Hs。eg。db）
　　28head（df）
　　29DEGdeg
　　30head（DEG）31：32DEGmerge（DEG，df，by。ySYMBOL，by。xsymbol）
　　33head（DEG）
　　34save（DEG，filedataannoDEG。Rdata）
　　35｝
　　36预处理2
　　37if（T）｛
　　38geneupDEG〔DEGgUP，ENTREZID〕
　　39genedownDEG〔DEGgDOWN，ENTREZID〕
　　40genediffc（geneup，genedown）
　　41geneallas。character（DEG〔，ENTREZID〕）
　　42data（geneList，packageDOSE）
　　43head（geneList）
　　44boxplot（geneList）
　　45boxplot（DEGlogFC）46：47geneListDEGlogFC
　　48names（geneList）DEGENTREZID
　　49geneListsort（geneList，decreasingT）
　　50｝51：52detailedplot
　　53if（T）｛
　　54source（keggandgoupanddown。R）
　　55runkegg（geneup，genedown，pronpcVSnormal）
　　56需要多go数据库的3个条目进行3次富集分析，非常耗时。
　　57rungo（geneup，genedown，pronpcVSnormal）
　　58｝59：60
　　61综合显示图
　　62if（F）｛
　　63goamp；lt；enrichGO（geneup，OrgDborg。Hs。eg。db，ontall）
　　64library（ggplot2）
　　65library（stringr）
　　66barplot（go，splitONTOLOGY）facetgrid（ONTOLOGY。，scalefree）
　　67barplot（go，splitONTOLOGY，font。size10）
　　68facetgrid（ONTOLOGY。，scalefree）
　　69scalexdiscrete（labelsfunction（x）strwrap（x，width50））
　　70ggsave（picgeneupGOallbarplot。png）71：72goamp；lt；enrichGO（genedown，OrgDborg。Hs。eg。db，ontall）
　　73barplot（go，splitONTOLOGY，font。size10）
　　74facetgrid（ONTOLOGY。，scalefree）
　　75scalexdiscrete（labelsfunction（x）strwrap（x，width50））
　　76ggsave（picgenedownGOallbarplot。png）
　　77｝
　　detailedplot
　　npcVSnormalkeggupdown。png（kegg结果还可以接受）
　　GO系列结果过于冗余：
　　npcVSnormaldotplotgenediffBP。png
　　npcVSnormaldotplotgenediffCC。png
　　npcVSnormaldotplotgenediffMF。png
　　npcVSnormaldotplotgenedownBP。png
　　npcVSnormaldotplotgenedownCC
　　npcVSnormaldotplotgenedownMF。png
　　npcVSnormaldotplotgeneupBP。png
　　npcVSnormaldotplotgeneupCC。png
　　npcVSnormaldotplotgeneupMF。png
　　综合显示图（更推荐）
　　genedownGOallbarplot。png
　　geneupGOallbarplot。png
　　step5annoGSEA
　　GSEA数据下载网页：https：www。gseamsigdb。orggseadownloads。jsp
　　msigdb。v7。0。symbols。gmt：https：www。gseamsigdb。orggseamsigdbdownloadfile。jsp？filePathresourcesmsigdb7。0msigdb。v7。0。symbols。gmt
　　msigdb。v7。0。entrez。gmt：https：www。gseamsigdb。orggseamsigdbdownloadfile。jsp？filePathresourcesmsigdb7。0msigdb。v7。0。entrez。gmt
　　所有打包gmt：https：www。gseamsigdb。orggseamsigdbdownloadfile。jsp？filePathresourcesmsigdb7。0msigdbv7。0filestodownloadlocally。zip
　　1载入数据和R包
　　2DEG为limma得到的差异分析结果
　　3if（T）｛
　　4rm（listls（））
　　5options（stringsAsFactorsF）
　　6load（filedataannoDEG。Rdata）
　　7library（ggplot2）
　　8library（clusterProfiler）
　　9library（org。Hs。eg。db）
　　10｝11：12
　　13对MigDB中的全部基因集做GSEA分析。
　　14按照FC的值对差异基因进行排序
　　15http：www。bioinfotrainee。com2105。html
　　16http：www。bioinfotrainee。com2102。html
　　17自行修改存放gmt文件路径d
　　18GSEA每个geneset的具体结果保存在gsearesults这个list中
　　19而最终结果保存在gsearesultsdf数据框中
　　20ddataGSEAgmtmsigdbv7。0filestodownloadlocallymsigdbv7。0GMTssymbols
　　21if（T）｛
　　22geneListDEGlogFC
　　23names（geneList）DEGsymbol
　　24geneListsort（geneList，decreasingT）
　　25选择gmt文件（MigDB中的全部基因集）26：27gmtslist。files（d，patternall）
　　28gmts29：30GSEA分析
　　31library（GSEABase）BiocManager：：install（GSEABase）
　　32下面使用lapply循环读取每个gmt文件，并且进行GSEA分析
　　33如果存在之前分析后保存的结果文件，就不需要重复进行GSEA分析。
　　34fdatagsearesults。Rdata
　　35if（！file。exists（f））｛
　　36gsearesultsamp；lt；lapply（gmts，function（gmtfile）｛
　　37gmtfilegmts〔2〕
　　38filepathpaste0（d，gmtfile）
　　39genesetamp；lt；read。gmt（filepath）
　　40print（paste0（Nowprocessthe，gmtfile））
　　41egmtamp；lt；GSEA（geneList，TERM2GENEgeneset，verboseFALSE）
　　42head（egmt）
　　43gseaplot（egmt，geneSetIDrownames（egmt〔1，〕））
　　44return（egmt）
　　45｝）
　　46上面的代码耗时，所以保存结果到本地文件
　　47save（gsearesults，filef）
　　48｝
　　49load（filef）
　　50提取gsea结果，熟悉这个对象
　　51gsearesultslistamp；lt；lapply（gsearesults，function（x）｛
　　52cat（paste（dim（xresult）），）
　　53xresult
　　54｝）
　　55｝
　　56随便看几个结果图
　　57dev。off（）
　　58gsearesultsdfamp；lt；do。call（rbind，gsearesultslist）
　　59gseaplot（gsearesults〔〔2〕〕，geneSetIDNIKOLSKYBREASTCANCER7P15AMPLICON）
　　60gseaplot（gsearesults〔〔2〕〕，RICKMANHEADANDNECKCANCERD）
　　step6annoGSVA
　　老实说，我对GSVA还不是很理解，但是知道在代码方面，做起来比较简单，有2个输入文件，一个是表达矩阵，一个是gmt文件即可。先把代码谢在这里，日后如果有需要，再仔细研究吧：
　　1对MigDB中的全部基因集做GSVA分析。
　　2还有ssGSEA，PGSEA
　　3载入数据
　　4if（T）｛
　　5rm（listls（））
　　6options（stringsAsFactorsF）7：8library（ggplot2）
　　9library（clusterProfiler）
　　10library（org。Hs。eg。db）
　　11load（filedatastep1output。Rdata）
　　12每次都要检测数据
　　13dat〔1：4，1：4〕
　　14｝15：16GSVA分析
　　17存放geneset的文件路径需要具体修改
　　18ddataGSEAgmtmsigdbv7。0filestodownloadlocallymsigdbv7。0GMTssymbols
　　19if（T）｛
　　20Xdat
　　21table（grouplist）
　　22MolecularSignaturesDatabase（MSigDb）
　　23gmtslist。files（d，patternall）
　　24gmts
　　25library（GSVA）BiocManager：：install（GSVA）
　　26if（！file。exists（datagsvamsigdb。Rdata））｛
　　27esmaxamp；lt；lapply（gmts，function（gmtfile）｛
　　28gmtfilegmts〔8〕；gmtfile
　　29genesetamp；lt；read。gmt（file。path（d，gmtfile））
　　30es。maxamp；lt；gsva（X，geneset，
　　31mx。diffFALSE，verboseFALSE，
　　32parallel。sz1）
　　33return（es。max）
　　34｝）
　　35adjPvalueCutoffamp；lt；0。001
　　36logFCcutoffamp；lt；log2（2）
　　37esdegamp；lt；lapply（esmax，function（es。max）｛
　　38es。maxesmax〔〔1〕〕
　　39table（grouplist）
　　40dim（es。max）
　　41designamp；lt；model。matrix（0factor（grouplist））
　　42colnames（design）levels（factor（grouplist））
　　43rownames（design）colnames（es。max）
　　44design
　　45library（limma）
　　46contrast。matrixamp；lt；makeContrasts（paste0（unique（grouplist），collapse），levelsdesign）
　　47contrast。matrixamp；lt；makeContrasts（npcnormal，
　　48levelsdesign）49：50contrast。matrix这个矩阵声明，我们要把progres。组跟stable进行差异分析比较51：52degfunction（es。max，design，contrast。matrix）｛
　　53step1
　　54fitamp；lt；lmFit（es。max，design）
　　55step2
　　56fit2amp；lt；contrasts。fit（fit，contrast。matrix）
　　57这一步很重要，大家可以自行看看效果58：59fit2amp；lt；eBayes（fit2）defaultnotrend！！！
　　60eBayes（）withtrendTRUE
　　61step3
　　62resamp；lt；decideTests（fit2，p。valueadjPvalueCutoff）
　　63summary（res）
　　64tempOutputtopTable（fit2，coef1，nInf）
　　65nrDEGna。omit（tempOutput）
　　66write。csv（nrDEG2，limmanotrend。results。csv，quoteF）
　　67head（nrDEG）
　　68return（nrDEG）
　　69｝70：71redeg（es。max，design，contrast。matrix）
　　72nrDEGre
　　73head（nrDEG）
　　74return（nrDEG）
　　75｝）
　　76gmts
　　77save（esmax，esdeg，filedatagsvamsigdb。Rdata）
　　78｝
　　79｝80：81画图展示，结果存放在pic下
　　82if（T）｛
　　83load（filedatagsvamsigdb。Rdata）
　　84library（pheatmap）
　　85lapply（1：length（esdeg），function（i）｛
　　86i8
　　87print（i）
　　88datesmax〔〔i〕〕
　　89dfesdeg〔〔i〕〕
　　90dfdf〔dfP。Valueamp；lt；0。01amp；gt；0。3，〕lt；！0。01gt；
　　91print（dim（df））
　　92if（nrow（df）amp；gt；5）｛
　　93nrownames（df）
　　94datdat〔match（n，rownames（dat）），〕
　　95acdata。frame（ggrouplist）
　　96rownames（ac）colnames（dat）
　　97rownames（dat）substring（rownames（dat），1，50）
　　98pheatmap：：pheatmap（dat，
　　99fontsizerow8，height11，
　　100annotationcolac，showcolnamesF，
　　101filenamepaste0（〔picgsva，strsplit（gmts〔i〕，〔。〕）〔〔1〕〕〔1〕，。pdf））102：103｝
　　104｝）105：106adjPvalueCutoffamp；lt；0。001
　　107logFCcutoffamp；lt；log2（2）
　　108dfdo。call（rbind，esdeg）
　　109esmatrixdo。call（rbind，esmax）
　　110dfdf〔dfP。Valueamp；lt；0。01amp；gt；0。5，〕lt；！0。01gt；
　　111write。csv（df，filedataGSVADEG。csv）
　　112｝
　　因为用到的这个样本用GSVA没有得到显著性结果，所以没有图出来，具体也没有深究，有需要日后再仔细研究吧
　　step7visualization
　　http：www。360doc。comconteant1803091833459258735717104。shtml
　　1主要是关于KEGG方面的扩展图
　　2主要是关于KEGG方面的扩展图
　　3hr4载入数据
　　5if（T）｛
　　6rm（listls（））
　　7options（stringsAsFactorsF）
　　8hr9library（ggplot2）
　　10library（clusterProfiler）
　　11library（org。Hs。eg。db）
　　12load（filedataannoDEG。Rdata）13：14head（DEG）
　　15｝16：17preprocessdata
　　18if（T）｛
　　19geneupDEG〔DEGgUP，ENTREZID〕
　　20genedownDEG〔DEGgDOWN，ENTREZID〕
　　21genediffc（geneup，genedown）
　　22geneallas。character（DEG〔，ENTREZID〕）
　　23制作差异基因listL
　　24boxplot（DEGlogFC）25：26geneListDEGlogFC
　　27names（geneList）DEGENTREZID
　　28geneListsort（geneList，decreasingT）
　　29｝30：31
　　32KEGG富集分析得到结果
　　33if（T）｛
　　34if（！file。exists（dataenrichkk。rdata））｛
　　35genedown
　　36geneup
　　37enrichKKamp；lt；enrichKEGG（genegeneup，
　　38organismhsa，
　　39universegeneall，
　　40pvalueCutoff0。1，
　　41qvalueCutoff0。1）
　　42save（enrichKK，filedataenrichkk。rdata）
　　43｝
　　44load（filedataenrichkk。rdata）
　　45查看KEGG结果
　　46head（enrichKK）〔，1：6〕
　　47打开网页看相关KEGG通路图
　　48browseKEGG（enrichKK，hsa05150）49：50将数据中的entrzid变成symbol
　　51更为易读
　　52enrichKKDOSE：：setReadable（enrichKK，OrgDborg。Hs。eg。db，keyTypeENTREZID）
　　53enrichKK
　　54｝55：56
　　57可视化
　　58条带图
　　59if（T）｛
　　60par（mfrowc（2，1））
　　61barplot（enrichKK，showCategory20）
　　62ggsave（picbarplot。png）
　　63｝64：65气泡图
　　66if（T）｛
　　67dotplot（enrichKK）
　　68ggsave（picdotplot。png）
　　69｝70：71下面的图需要映射颜色，设置和示例数据一样的geneList72：73展示top5通路的共同基因，要放大看。
　　74GeneConceptNetwork
　　75if（T）｛
　　76cnetplot（enrichKK，foldChangegeneList，colorEdgeTRUE，circularF）
　　77ggsave（piccnetplot。png）
　　78cnetplot（enrichKK，foldChangegeneList，colorEdgeTRUE，circularT）
　　79ggsave（piccnetplotcircular。png）
　　80｝81：82
　　83EnrichmentMap
　　84if（T）｛
　　85emapplot（enrichKK）
　　86ggsave（picEnrichmentMap。png）
　　87｝88：89（4）展示通路关系，仅仅是针对于GO数据库结果。
　　90goplot（enrichKK）
　　91（5）Heatmaplikefunctionalclassification
　　92if（T）｛
　　93heatplot（enrichKK，foldChangegeneList）
　　94ggsave（picEnrichmentHeatmap。png）
　　95｝
　　条带图：
　　气泡图：
　　GeneConceptNetwork图：每一个小蓝圈表示一个基因，其颜色表示FC值，每个KEGGterm圈的大小由里面包含基因的数目决定。
　　成环：更炫酷了，但是感觉图形展示不方便
　　不成环：信息展示更有力吧
　　EnrichmentMap图：这里和上面的图类似，只不过不再显示具体的基因，而是直接画出每个term和term之间的关系，每个圆圈代表着一个term，圆圈大小代表着有多少个基因，颜色表示p值。
　　如果term和term之间有共同的基因，那么就会连接起来，聚在一起。
　　Heatmaplikefunctionalclassification：
　　和我们常规的热图不太像，这里纵轴是每个KEGG通路，横轴是涉及到的基因。颜色表示FC值。
　　致谢
　　上面所有的代码都来自生信技能树曾老板jimmy的帮助，同时我在测试运行的过程中又进行了部分改进和增补。
　　就用曾老板亲自编辑的感谢词来感谢吧：
　　FatYangthankDr。JianmingZeng（UniversityofMacau），andallthemembersofhisbioinformaticsteam，biotrainee，forgenerouslysharingtheirexperienceandcodes！