032-熟悉一下Bioconductor

Sep 7, 2018 25 min read cnposts

刘小泽写于18.9.7
各行各业数据呈爆炸式增长，大量的数据等待被处理，R语言就是一个利器，可以说是做数据分析必备的编程语言。当强大的R与包罗万象的生物结合，再一次刺激了R的迅猛发展。随着NGS测序的普及，R语言的生信专业社区Bioconductor诞生，开启了生物信息的R语言时代。好工具，用起来，首先要了解生物数据与R之间的关联

生物知识回顾

基因有三类：第一类编码序列，编码蛋白【转录+翻译功能】；第二类只有转录没有翻译功能【tRNA+rRNA】；第三类不转录基因，调控基因表达【启动子、操纵子】
基因组：单倍体细胞中包含编码序列和非编码序列的全部DNA【核基因组+线粒体基因组+叶绿体基因组】
狭义转录组：mRNA，代表样本整体基因表达水平，一般称为“表达谱”【真正的表达谱是蛋白质组信息，研究手段是质谱】。通常用基因芯片、RNA-seq
广义转录组：
- 非编码ncRNA：有三类【按长度划分】
  - 小于50nt的small RNA【长度单位：单链叫nt，双链叫bp】：miRNA、siRNA、piRNA【small RNA序列短，同源性高，一般采用二代测序检测】
  - 50-500nt：rRNA、tRNA、snRNA、snoRNA
  - 大于500nt：mRNA-like ncRNA、不带polyA尾的ncRNA
- microRNA：也叫miRNAs，20-25nt，初级转录物-》核酸酶剪切加工-〉组装进RNA诱导的沉默复合体-》互补配对识别靶mRNA-〉根据互补程度不同指导沉默复合体降解靶mRNA或者阻遏mRNA翻译
DNA、蛋白互作：Chromatin immuopre-cipitation, ChIP，主要应用于：DNA序列转录因子结合位点（Binding sites）识别，如启动子、增强子等顺式作用文件(Cis-acting element)的识别；DNA甲基化、组蛋白修饰、核小体定位
DNA甲基化：甲基化DNA免疫共沉淀测序(Methylated DNA immunoprecipitation sequencing, MeDIP-seq)、甲基化DNA（蛋白）结合域测序(Methylated DNA binding domain sequencing, MBD-seq)和亚硫酸氢盐测序(Bisulfite sequencing, BS-seq)

基因表达分析

基因表达检测方法

实时荧光定量PCR(Quantitative real time PCR, qRT-PCR)
基因（表达谱）芯片（Microarray）
表达序列标签（Expressed Sequence Tag, EST）
基因表达系列分析(Serial Analysis of Gene Expression, SAGE)
转录组测序

PCR技术应用最为成熟，灵敏度高，特异性强，但其缺点是通量较小；基因芯片方便快捷，适合临床诊断及个体基因组分析；基因测序技术通量高，但周期长、成本高

基因芯片是什么

基因芯片又称DNA微阵列，按照检测物的不同，可分为DNA芯片、RNA芯片等，其中DNA芯片又可分为单核苷酸多肽性（SNP）芯片、比较基因组杂交（CGH）芯片等。

原理：基于A、T；C、G互补理论，将已知序列的核酸探针与未知序列的核酸序列进行杂交检测DNA，并且DNA探针以显微打印的方式大规模集成于芯片（类似于计算机的硅芯片）表面。杂交后通过计算机对杂交信号的检测分析，得出样品的遗传信息(基因序列及表达的信息)。分析单核苷酸变异多态性性价比较高。

主流寡聚核苷酸芯片主要有：Affymetric、Agilen、Illumina公司

基因表达数据

矩阵表示：行名代表一个基因不同条件/样本的表达，列名代表某个条件/样本的所有基因表达。数据代表表达水平。那么一般分析什么？

不同样本/处理中哪些基因表达有显著差异？
基因之间有什么共有的功能，或者参与哪些共同代谢途径？
不同的处理中，哪些基因变化一致，它们受到上游哪些基因的调节，或者它们控制下游哪些基因的表达？
哪些基因表达存在样本特异性，也就是说通过他们的表达可以判断样本的状态（如：细胞的增殖、分化、凋亡、应激、癌变等）

主要的分析

主要有差异显著性分析和时间序列分析，后者主要是测定基因多个时间点的表达量，然后聚类+主成分分析寻找共调控基因

表达显著性分析就是为了找差异基因（DEG）。那么怎样判断基因间是有差异的呢？常用的有3种算法：一是倍数分析（无统计假设），计算每个基因在不同条件/样本的比值，再与阈值比较；二是用统计模型T检验等方法，计算差异表达的置信度p值，以0.05或者0.01作为阈值；三是机器学习方法，利用贝叶斯模型、随机森林等。分析的结果从来不用担心没有差异基因，而是要考虑差异基因可能存在很多，从几十个到上百个不等，那么如何展示他们呢，一张简单粗暴的大表格吗？肯定是不行的！

需要把上游的这些差异基因再进行注释、分组，一个类别就相当于一个GO term，然后看这几大类的区别，肯定比看几十甚至上百个基因或蛋白的差异要更加直观，这就是富集分析，包括GO分析，KEGG分析，GSEA分析等。其中重点研究的基因集叫做前景基因，需要比对的所有基因集叫背景基因，前景是背景的子集。例如转录组数据中的对照组和处理组，处理与对照之间的差异基因就是前景基因，两组所有的表达基因就是背景基因。富集分析的目的就是根据不同功能，把各个分子进行分类，然后使用超几何分布检验进行分析。当然使用不同工具，得到的结果不同，现在clusterProfiler要比DAVID的结果更多。

GO分析（Gene Ontology）

包括GO terms（标签）+GO annotations（注释）。

GO terms存在于由基因本体联合会(Gene Ontology Consortium)建立的数据库中，对基因和蛋白功能进行限定和描述，每个注释信息都有一个GO ID。它由两部分构成，第一部分都是GO，第二部分是以0开头的7位数字，例如GO:0016021 。
GO是一个情报员，他负责调查：包括基因的分子功能：“干啥的”（molecular function，MF），指分子所执行的任务【如与碳水化合物结合或ATP水解酶活性等】、细胞组分：“活动区域”即产物发挥作用的位置（cellular component，CC）【如核仁、端粒和识别起始的复合物】、参与的生物过程：“近期有什么动静”（biological process，BP）【嘌呤代谢、有丝分裂等】GO调查完就给被查对象贴标签term
GO annotations即GO注释，是针对基因产物的而不是基因，表示某些基因的产物是是非编码RNA、蛋白质还是大分子等。这里GO就相当于一个中间媒介，它对基因进行定义GO term，然后其他各个数据库使用GO的定义方法，对它们的基因产物进行标注，例如一个数据库的EntrezID或SYMBOL与GO数据库进行ID对应，或者用一个数据库的序列与GO term进行对应。

在一个GO注释中，例如，一个基因的产物是细胞色素c（cytochrome c），那么这个基因的产物就会被一个分子功能术语（Molecular Function）描述为氧化还原酶活性（oxidoreductase activity ），被生物过程（Biological Process）描述为氧化磷酸化（oxidative phosphorylation ），被细胞成分（Cellular Component ）描述为线性体基质（mitochondrial matrix ）和线粒体内膜（mitochondrial inner membrane ）~引用自“读研笔记”

pathway代谢通路

GO负责分门别类，而pathway负责把每一类对应到具体的代谢网络中。研究pathway的原因是：生物学问题中设定一个“蝴蝶效应”假设：1个Pathway上游基因的改变，会导致下游相关基因改变，从而改变通路中大量基因的表达。现在常用是KEGG，但是它收录的都是是已有的研究结果，而这些信息，还没有完善

熟悉一下Bioconductor

Bioconductor vs R

Bioconductor源于R但高于R，Bioconductor能做到的一般的R编程也可以，但是Bioconductor的优势就在于一个简单的包就能干R几十行代码做出来的事，编写Bioconductor R包的人花费大量精力，将数据-算法-注释紧密结合起来，这个是一般的R编程不能替代的。因此会一个Bioconductor包的使用，比会编程还要重要
Bioconductor拥有上千个扩展包，主要有实验数据包、软件包、注释数据包三大类，例如白血病的ALL包就是利用Affymetrix进行芯片分析的数据包；但最重要的当属软件包

软件包：

注释：GO、Pathway等
微阵列板块（Assay Domains）：处理芯片数据，Bioconductor支持主流的Affymetrix的商业化单色寡聚核苷酸芯片，也支持用户定制的双色cDNA芯片。芯片数据一般流程：数据预处理、差异表达基因筛选、聚类分析。这里的包有以下几部分：
- 比较基因组杂交（Comparative Genomic Hybridization, CGH）
- 细胞水平检测（Cell Based Assays）
- 染色质免疫共沉淀芯片（ChIPchip）
- 拷贝数变异（Copy Number Variants）
- CpG岛（CpGIsland）
- 差异表达（Differential Expression）
- DNA甲基化（DNA Methylation）
- 外显子检测（Exon Assay）
- 基因表达（Gene Expression）
- 遗传变异性（Genetic Variability）
- 单核苷酸多态性（SNP）
- 转录
测序技术（Assay techs）
- 芯片技术（Microassay）
- 微孔板检测（Microtitre Plate Assayå）
- 质谱（Mass Spectrometry）
- 基因表达系列分析（SAGE）
- 流式细胞仪（Flow Cytometry）
- NGS
数据处理：基因芯片数据预处理（背景矫正、归一化、质控）、芯片分析、基因间关系、样本间关系、识别差异基因
聚类分析（Clustering）、分类（Classification）、富集分析（Enrichment）、多组比较（Multiple Comparison）、预处理（Preprocessing）、质控、序列匹配、时间序列分析（Time Course）、可视化、网络分析

应用领域之基因芯片

Bioconductor的起源就是基因芯片分析，因此提供了芯片数据库如GEO、ArrayExpress结口，方便获取数据

芯片提供商：Affymetrix、Illumina、Nimblegen、Agilent

Affemetrix 3‘-biased Arrays芯片数据需要affy、gcrma、affyPLM包，这三个包又依赖于芯片定义（CDF）、探针（Probe）、注释（Annotation）三个包，后三个会随前三个安装而自动装好；另外还需要手动安装ROOT包，这个包可以直接使用Affymetrix的数据文件，支持的格式有CDF、PGF、CLF、CSV
Affymetrix Exon ST Arrays与Affymetrix Gene ST Arrays芯片需要oligo包和xps包。Ologo包又需要调用pdInfoBuilder包来创建一个pdInfoPackage包，pdInfoPackage包的作用就是合并CDF、Probe、Annotation三种数据
Affymetrix SNP Arrays、Affymetrix Tilling Arrays、Nimblegen Arrays芯片数据处理只需要oligo包
Ilumina Expression Microarrays芯片数据需要lumi和beadarray包，这两个包都有各自的映射包和注释包：lumi需要lumiHumanAll.db、lumiHumanIDMapping包；beadarray需要illuminaHumanv1BeadID.db、illuminaHumanv1.db包

芯片类型：基因表达芯片、外显子芯片、拷贝数变异检测芯片、SNP芯片、DNA甲基化芯片等

芯片分析内容：预处理、质量评估、差异基因表达分析、基因集富集分析、遗传基因组学

应用领域之测序数据

R也可以处理多种类型数据：fasta、fastq、SAM、BAM、Gff、Bed、Wig等，可以进行测序结果预处理（去除低质量、序列污染等）、格式转换、序列比对、测序质量评估、RNA-seq、差异表达分析、ChIP-seq等

SRAdb包是SRA数据库的接口
ShortRead包读入测序数据，进行质控、预处理；Rsamtools将数据比对到参考基因组/转录组，完成格式转变和比对结果统计；rtracklayer包将数据导入USCS基因组浏览器(http://genome-asia.ucsc.edu/)进行浏览、操作、输出
IRanges、GenomicRanges、genomeIntervals包是基于DNA区域（如染色体）进行数据操作的拓展包；Biostrings包含了序列比对、模式匹配等基本序列分析函数；BSgenome针对有注释的全基因组数据进行操作；GenomicFeatures是对基因组中序列特征（如非编码区）进行注释
RNA-seq可用的包有：limma（金标准）、edgeR、DESeq、DEXSeq（外显子水平）、baySeq
ChIP-seq：基序发现(Motif discovery)、结合位点检测（Peak calling），可以用CSAR、chipseq、ChIPseqR、ChIPsim、ChIPpeakAnno、DiffBind、iSeq、rGADEM、segmentSeq、BayesPeak、PICS

更多的数据处理，可以看https://www.bioconductor.org/packages/release/BiocViews.html=》softeware=〉technology

应用领域之注释

注释既需要注释工具，有需要链接注释数据库的接口，大体的注释方式有三类

第一类：基于AnnotationDbi包的扩展包

绝大部分注释包都依赖AnnotationDbi包，一旦biocLite安装任何一个“.db”的注释包，AnnotationDbi包都会自动安装。这个包可以创建、操作“.db”，并且可以创建个性化的芯片平台。它又可以分为三类注释包

物种型注释包：org.Mm.eg.db、org.Hs.eg.db、org.Rn.eg.db等包括整个物种相关的基因数据，命名格式：“org.Xx.yy.db”， Xx是种属的缩写（注意首字母大写），“yy“来源数据库ID，这个ID用于把所有数据连在一起。如：“org.Hs.eg.db”中Hs是种属，eg是数据库ID
平台型注释包：如Affymetrix公司的hgu133plus2.db 依赖于具体的实验平台，命名格式：“platformName.db”。“platform”是芯片平台名称，如“hgu95av2.db”就是利用了Affymetrix公司的hgu95av2基因芯片；另外每个“.db”包还需要对应的“.cdf”、“.probe”，名称分别为：“platformName.cdf”、“platformName.probe”
系统生物学注释包：用于下游分析，如基因功能分析（与上游），比如KEGG.db用来获取KEGG数据库的数据；GO.db用来获取GO数据库的数据；PFAM.db用来获取不同蛋白家族的特征和相关性

第二类：基于biomaRt包获取注释

它与AnnotationDbi包虽然注释信息都是来自远程数据库，但是AnnotationDbi是将注释包本地化，而biomaRt只提供了一组数据接口，加载注释内容严重依赖网络，并且无需维护，不占用本地存储，注释信息自动更新。但它的功能最为强大，因为可以连接各个主要的是给你想你数据库来获取注释信息和数据

应用领域之高通量实验

包括流式细胞仪（Flow Cytometry）、定量PCR、质谱（Mass soectrometry）、蛋白质组及其他基于细胞水平的高通量实验数据

定量PCR：HTqPCR、ddCt、qpcrNorm提供如何分析循环阈值（Cycle threshold）的方法

利用Bioconductor分析芯片数据

芯片分析主要就是利用bioconductor，分析的过程也体现了它的设计理念和编程思想

基因芯片背景知识

以Affymetrix为例，一个芯片可以包括上百万探针（通常由25个碱基组成），被整齐印刷在芯片上。

探针组：一个探针组通常由20对个或11个探针对组成，来自一个基因
探针对：匹配探针（Perfect match， PM）+错配探针（Mismatch， MM）。二者仅仅是中间的那个碱基不同。并非所有芯片都有这两个探针，比如外显子芯片，每个探针组只有4个PM探针，没有MM探针
探针序列来源：公共核酸数据库中的参考序列（如NCBI的GenBank或RefSeq）
不同的芯片中，探针组在参考序列中的分布不同，也就是检测范围不同。3‘表达谱芯片探针组排列在参考序列3’末端附近的一至两个外显子上；外显子芯片中，每个长度大于25个碱基的外显子都有对应的探针组；铺瓦芯片（Tilling array）中，探针基本覆盖了所有的外显子和内含子
探针组与基因关系：芯片数据得到的表达矩阵，实际上是以探针组为单位，而不是直接以基因为单位，每一行对应一个探针组的表达量。后期的分析都是先得到探针组的结果，然后根据注释的ID映射才对应到基因。一般是一个基因同时对应多个探针组。通常会把同一个基因对应的探针组表达量求均值，然后找最大的那个探针组作为代表，让它与该基因一一对应
芯片数据获取：一是扫描设备对芯片进行扫描，得到荧光信号图像文件（DAT文件）；二是由系统自带的图像处理软件，经过网格定位（Griding）、杂交点范围确定（Spot identifying）、背景噪音过滤（Noise filtering），从芯片图像中提取数据，得到CEL文件
关于Affymetrix的数据处理过程：https://slideplayer.com/slide/4804237/
CEL文件：CEL文件是Affymetrix最常用的格式，它的主要内容就是每个“cell”的灰度信息，而“cell”就是整个芯片图像划分后得到的小网格，每个小网格中的图像被看作来自一个探针。【补充】当然，如果要得到芯片上每个探针组对应的表达数据，还需要探针的排布信息（哪个探针对应哪个探针组），这部分信息就存储在CDF文件中。要想知道探针对应的序列信息，就需要用Probe文件

完成一个小任务

从CLL数据包载入芯片数据=》预处理=〉获得探针组表达矩阵

先说一下这个CLL数据集，它包含的是慢性淋巴细胞白血病（Chromic lymphocytic leukemia, CLL）的数据，采用Affymetrix的hgu95av2表达谱芯片（含有12625个探针组），共24个样品（“CLL1.CEL”-“CLL24.CEL”），每个样品来自不同的癌症病人。根据健康状态分为两组（对照试验）：稳定组（stable）和恶化组（progressive），每组各12个（平行试验）。得到的e矩阵是一个12625行、24列的探针组表达矩阵

#安装并加载
source("https://bioconductor.org/biocLite.R")
options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")
biocLite("CLL")
library(CLL)
#读入示例数据
data("CLLbatch")
# rma方法进行背景校正【当MM值比PM值还要高时，MM就是杂信号，也就是背景噪声，需要去除】
CLLrma <- rma(CLLbatch)
e_before <- exprs(CLLbatch)
e_after <- exprs(CLLrma)
#对比一下校正前后数据
e_before[1:5,1:5]
e_after[1:5,1:5]

实战～

1 下载芯片数据

使用Fabbrini E于2015年发表在J Clin Invest上的文章：Metabolically normal obese people are protected from adverse effects following weight gain
文章研究了体重适度增加对代谢正常(MNO)和代谢异常(MAO)受试者脂肪组织基因表达的影响，选取11个代谢正常的人、7个代谢异常的人增重前后的皮下脂肪组织，得到36个数据
文章下载链接：https://pdfs.semanticscholar.org/6052/f035dfb296559e19314352506f048e13f02e.pdf
使用的数据集是GSE62832
实验平台是GPL6244， Affymetrix Human Gene 1.0 ST Array

方法一：打开网页https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE62832 【可能会有网速限制】

方法二：使用GEOquery包

if(T){
source("http://bioconductor.org/biocLite.R")
options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")
biocLite("GEOquery")
  
library(GEOquery)
eSet <- getGEO('GSE62832', destdir = '.', getGPL = F, AnnotGPL = F)#一般芯片的注释文件比较大，直接下载不容易成功，后两个选项可以避免下载GPL、Annotation
save(eSet, file = 'GSE62832.eSet.Rdata')
}

下载完，看一下eSet这个对象都包含什么

> eSet
$GSE62832_series_matrix.txt.gz
ExpressionSet (storageMode: lockedEnvironment)
assayData: 33297 features, 36 samples 
  element names: exprs 
protocolData: none
phenoData
  sampleNames: GSM1534034 GSM1534035 ... GSM1534069 (36 total) #36个样本
  varLabels: title geo_accession ... tissue:ch1 (35 total)
  varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
Annotation: GPL6244 #查看对应芯片平台

class(eSet) #先查看数据种类【列表还是数据框】
str(eSet) #再看数据结构【可以看到第二行涉及到了表达矩阵】
#既然是列表，那么从列表中提取信息就是 eSet[[1]]
e <- exprs(eSet[[1]])

2. 将表达矩阵的探针ID转换为gene ID

首先需要知道GSE62832对应平台是GPL6244
然后需要知道GPL6244对应哪个注释包【阅读jimmy之前整理的平台信息https://www.jianshu.com/p/f6906ba703a0】

下载相应的注释包

source("https://bioconductor.org/biocLite.R")
options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")
biocLite("hugene10sttranscriptcluster.db")

开始探索、过滤、整合

library(hugene10sttranscriptcluster.db) 
ls("package:hugene10sttranscriptcluster.db") #查看所有包含的对象
s <- toTable(hugene10sttranscriptclusterSYMBOL) #将SYMBOL对象转为数据框
#【补充：这个注释包是别人上传的，我们只是拿来用。其中的表达矩阵中原来有3w多探针，作者过滤后只留下接近2w的探针。这仅仅是作者过滤后的结果，并非原始数据结果】

简单的探索

s$symbol %>% unique() %>% length() #看一下有多少基因【人类正常蛋白编码基因就2w左右】
s$symbol %>% table() %>% sort() %>% tail() #基因使用探针最多的前6名【发现有两个基因用了10个探针】
s$symbol %>% table() %>% sort() %>% table()

过滤+整合

dim(e) #过滤前
#过滤【只保留和注释文件探针id相同的探针】
efilt <- e[rownames(e)%in%s$probe_id,]
dim(efilt)#过滤后
#整合1【目的：保证一个基因对应一个探针；如果基因和探针一一对应很好说，但如果一个基因对应多个探针：每个探针取一行的均值-》对应同一基因的探针取表达量最大的探针-》按照基因名给他们建索引，因为是按照基因来过滤探针（不用s$probe_id构建索引的原因是，看清楚我们的目的是让注释包的一个基因对应我们自己表达矩阵的一个探针。如果用s$probe_id那么结果就成了让注释包的一个探针对应我们自己表达矩阵的一个探针，当然这样也运行不成功，因为自己表达矩阵的探针过滤后的数量和注释包的探针数量不相等，这样没法一一对应。但基因名数量是不变的，什么是索引？以不变应万变的就是索引）】
maxp = by(efilt,s$symbol,function(x) rownames(x)[which.max(rowMeans(x))]) 
uniprobes = as.character(maxp)
efilt=efilt[rownames(efilt)%in%uniprobes,]
#整合2【目的：将我们表达矩阵的行名换成刚才一对一的基因名，并且match这个函数保证了表达矩阵和注释包的顺序是一致的】
rownames(efilt)=s[match(rownames(efilt),s$probe_id),2]

3. 对表达矩阵进行一些检验

表达矩阵不局限于GEO数据库的芯片分析，转录组及其他涉及基因、样本、分组关系的都会有一个表达量矩阵，就是一个基因在不同样本中（对照、处理；是否患病等）的表达差异。拿到表达矩阵，在进行后续分析之前，首先要检测这个矩阵是不是合理的，比如看管家基因是否表达量突出、一致；样本分组是不是和实验设计一致，用PCA、hclust检验

3.1 检测一些管家基因表达量

boxplot(efilt[,1])#看看第一个样本中总体基因表达量分布，可以看到基本为5左右
efilt['ACTB',] #激动蛋白Beta-actin的基因名是ACTB，管家基因
efilt['GAPDH',] #也是管家基因

3.2 看表达矩阵的整体分布

先把表达矩阵=》tidy data【四列：基因名、样本、表达量、表型分组(看文献按MAO、MNO分组)】

library(reshape2)
pdata=pData(eSet[[1]]) #将样本表型信息从数据框中提取出来【取出来的是表型、样本的数据框】
group_list=as.character(pdata$`metabolic status:ch1`) 
m_efilt = melt(efilt) #先将原来矩阵“融化
colnames(m_efilt)=c('symbol','sample','value') #重新命名三列
m_efilt$group=rep(group_list,each=nrow(efilt))

以图为证

4. 对表达矩阵进行差异分析

只需要提供表达矩阵efilt、分组信息group_list，就能使用limma进行分析详细内容参考jimmy的http://www.bio-info-trainee.com/bioconductor_China/software/limma.html

suppressMessages(library(limma))
#limma需要三个矩阵：表达矩阵（efilt）、分组矩阵（design）、比较矩阵（contrast）
#先做一个分组矩阵～design，说明MAO是哪几个样本，MNO又是哪几个，其中1代表“是”
design <- model.matrix(~0+factor(group_list))
colnames(design) <- levels(factor(group_list))
rownames(design) <- colnames(efilt)
design
#再做一个比较矩阵【一般是case比control】
contrast<-makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)
contrast

4.1 准备就绪，就可以开始差异分析

DEG <- function(efilt,design,contrast){
  ##step1
  fit <- lmFit(efilt,design)
  ##step2
  fit2 <- contrasts.fit(fit, contrast) 
  fit2 <- eBayes(fit2)  
  ##step3
  mtx = topTable(fit2, coef=1, n=Inf)
  deg_mtx = na.omit(mtx) 
  return(deg_mtx)
}
DEG_mtx <- DEG(efilt,design,contrast) #得到全部的差异基因矩阵

4.2 对一个小表达矩阵，如30、50个基因，可以用热图

top30_gene=head(rownames(DEG_mtx),30)
top30_matrix=efilt[top30_gene,] #得到top30的表达量矩阵
top30_matrix=t(scale(t(top30_matrix)))
#这里做个top30 gene heatmap

4.3 对一个大的表达矩阵，如全部的差异基因，可以用火山图

火山图实际上就是根据两列进行作图：logFC、pvalue

plot(DEG_mtx$logFC, -log10(DEG_mtx$P.Value)) #先画一个最简单的图（能说明原理）
#再根据jimmy的代码做一个。网络上也有很多火山图的代码可以参考

5 对结果进行注释～富集分析（基于超几何分布检验）

5.1 得到上调、下调基因

rm(list=ls())
load("GSE62832.DEG.Rdata")
DEG=DEG_mtx
if(T){
  logFC_cutoff <- with(DEG,mean(abs(logFC)) + 2*sd(abs(logFC)) ) 
  DEG$result = as.factor(ifelse(DEG$P.Value < 0.05 & abs(DEG$logFC) > logFC_cutoff,
                                ifelse(DEG$logFC > logFC_cutoff ,'UP','DOWN'),'NOT'))
}
gene <- rownames(DEG[DEG$result != 'NOT', ])

要再具体细分的话也可以

up_gene <- rownames(DEG[DEG$result =='UP',])
down_gene <- rownames(DEG[DEG$result =='DOWN',])

5.2 ID转换

使用clusterProfiler包接受的gene是Entrez ID，因此要先转换SYMBOL=》Entrez

几种主流的ID：

Ensemble id：由欧洲生物信息数据库提供，一般以ENSG开头，后边跟11位数字。如TP53基因：ENSG00000141510 Entrez id：由美国NCBI提供，通常为纯数字。如TP53基因：7157 Symbol id：为我们常在文献中报道的基因名称。如TP53基因的symbol id为TP53 Refseq id：NCBI提供的参考序列数据库：可以是NG、NM、NP开头，代表基因，转录本和蛋白质。如TP53基因的某个转录本信息可为NM_000546

library(clusterProfiler) 
library(org.Hs.eg.db)
if(T){
  gene_tr <- bitr(gene, fromType = "SYMBOL",
                  toType = c("ENSEMBL", "ENTREZID"),
                  OrgDb = org.Hs.eg.db)
  up_gene_tr <- bitr(up_gene, fromType = "SYMBOL",
                     toType = c("ENSEMBL", "ENTREZID"),
                     OrgDb = org.Hs.eg.db)
  down_gene_tr <- bitr(down_gene, fromType = "SYMBOL",
                       toType = c("ENSEMBL", "ENTREZID"),
                       OrgDb = org.Hs.eg.db)
  gene <- gene_tr$ENTREZID
  up_gene <- up_gene_tr$ENTREZID
  down_gene <- down_gene_tr$ENTREZID
}
head(gene_tr)
#结果发现：总体有3.07%没有对应，上调有2.13%没有对应，下调3.76%没有对应

转换完可以去检验一下，比如检查Entrez ID: 23336, https://www.ncbi.nlm.nih.gov/gene/23336。确实得到的是SYNM这个基因

5.3 进行KEGG注释

最常见的就是GO/KEGG数据库注释，根据clusterProfiler，每个注释就对应一个函数

富集分析使用超几何分布检验：比如背景基因这里有18837个，选出来的差异基因（上调235个+下调319个）554个，约占总体3%；目前KEGG 包含了530个pathway，其中一个04110是Cell Cycle通路，其中包含124个基因【怎么统计？打开链接-》https://www.kegg.jp/dbget-bin/www_bget?hsa04110，找到gene栏-》使用linux cat+wc】如果是随机抽样的话，cell cycle在我们这554个基因中，应该能抽到124*3%=3个左右。也就是说，在得到的554个差异基因中，应该只有3个是cell cycle的基因；但结果找到了30个cell cycle的基因，那么就说明差异基因的产生都有cell cycle这个因素。如果是药物处理组的结果发现这样，就说明确实可能是药物起了作用，当然后续还需要结合其他分析并非所有的通路都需要检测：有的通路比较小，只有十几个基因，我们这500多个差异基因中能抽到的理论上最多一个，这样一般就不考虑；还有的通路很大，有上千个基因，我们的背景基因才不到2w，那我们这586个差异基因中在这个pathway中的理论上就得有一百多个，这样理论和实际的变化差异就不明显，说明不了什么问题

enrichKEGG进行候选基因通路分析

if(T){
  kk <- enrichKEGG(gene         = gene,
                   organism     = 'hsa',
                   pvalueCutoff = 0.05)
  kk_gene <- (kk)[,1:6]
}
dotplot(kk,showCategory = 14,color="pvalue",font.size=14)

按照包的说明做一个geneList（标准：Entrez ID对应logFC，然后logFC降序排列）

if(T){
  geneList <- DEG_mtx$logFC #先获得logFC【没有重复的18837个数值】
  names(geneList) <- rownames(DEG_mtx) #获得基因名，但这里是SYMBOL格式【也是没有重复的18837个名称】
  ##问题就出在（下一步）SYMBOL转Entrez名字上，会多出来许多的重复：如果只是SYMBOL转Entrez，那么会得到18841个名称，原因就是一个symbol对应多个entrez ID，用count(geneList_tr,geneList_tr$SYMBOL) %>% arrange(desc(n))这样来查看到底哪些基因，结果查到有4个基因对应了2个entrez ID，其中一个是“HBD”，然后使用filter(geneList_tr, geneList_tr$SYMBOL == "HBD")来看HBD对应了哪个ID，果然发现HBD对应了两个ID；如果将SYMBOL转ENSEMBL+Entrez，那么会得到20198个名称，因为同一个SYNBOL，同一个Entrez，会有不同的Ensembl ID
  #另外还有别人的友情提示：ENTREZID(具有唯一性）：基因编号来自ANNOVAR的注释结果，建议别用SYMBOL，因为这种名称特异性较差，在转成ENTREZID时可能出现不唯一的现象。symbol与entrezID并不是绝对的一一对应的
  geneList_tr <- bitr(names(geneList), 
                      fromType = "SYMBOL",
                      toType = c("ENSEMBL","ENTREZID"),
                      OrgDb = org.Hs.eg.db) 
  #下一步目的：将Entrez ID与logFC联系起来
  #因为转换后的名称已经不是和原来的geneList行名一一对应了，所以不能直接将原来SYMBOL行名变为Entrez ID。但是可以用SYMBOL名作为桥梁，一边连接logFC（做一个数据框），另一边连接Entrez ID（已经有的geneList_tr数据框），再将两个数据库按照SYMBOL联系起来
  new_list <- data.frame(SYMBOL=names(geneList), logFC = as.numeric(geneList)) 
  new_list <- merge(new_list, geneList_tr, by  = "SYMBOL")
  
  geneList <- new_list$logFC #重新生成geneList之logFC
  names(geneList) <- geneList_tr$ENTREZID #重新生成geneList之Entrez ID
  
  geneList <- sort(geneList,decreasing = T) #降序排列
}

gseKEGG GSEA分析

if(T){
  kk2 <- gseKEGG(geneList     = geneList,
                 organism     = 'hsa',
                 nPerm        = 1000,
                 minGSSize    = 120,
                 pvalueCutoff = 0.2,
                 verbose      = FALSE)
}
head(kk2)[,1:6]
gseaplot(kk2, geneSetID = "hsa04310")

查看特定通路图

library(pathview)
hsa01212 <- pathview(gene.data  = geneList,
                     pathway.id = "hsa01212", #上述结果中的hsa01212通路
                     species    = "hsa",
                     limit      = list(gene=max(abs(geneList)), cpd=1))

5.4 进行GO注释

if(T){
  ego_CC <- enrichGO(gene          = gene,
                     universe      = names(geneList),
                     OrgDb         = org.Hs.eg.db,
                     ont           = "CC",
                     pAdjustMethod = "BH",
                     pvalueCutoff  = 0.01,
                     qvalueCutoff  = 0.05,
                     readable      = TRUE)
}
ego_gene <- head(ego_CC)[,1:6]
##可视化--点图
dotplot(ego_CC,title="EnrichmentGO_CC_dot")#点图，按富集的数从大到小的
##可视化--条形图
barplot(ego_CC, showCategory=20,title="EnrichmentGO_CC")#条状图，按p从小到大排，绘制前20个Term
#通路图
library(Rgraphviz)
library(topGO)
plotGOgraph(ego_CC)