基因表达分析（上） 基因表达

什么是基因表达，如下是来自于维基百科的解释：

Gene expression is the process by whichinformation from agene is used in thesynthesis of a functionalgeneproduct.These products are oftenproteins, but in non-proteincoding genes such astransferRNA (tRNA)orsmallnuclear RNA (snRNA)genes, the product is a functionalRNA. The process of geneexpression is used by all known life—eukaryotes (includingmulticellularorganisms),prokaryotes (bacteria andarchaea), and utilized byviruses—to generate themacromolecular machinery for life.

Flowof genetic information

研究方法

定量PCR

这部分我不太懂，所以就放几段百度百科和维基百科的定义。

百度百科

定量PCR（即时聚合酶链锁反应,Real-time Polymerase Chain Reaction，简称Real-time PCR、即时PCR），

又称定量即时聚合酶链锁反应（Quantitativereal time polymerase chain reaction，简称 Q-PCR/qPCR/rt-qPCR、定量即时PCR、即时定量PCR），是一种在DNA扩增反应中，以萤光染剂侦测每次聚合酶链锁反应（PCR）循环后产物总量的方法技术，有广义概念和狭义概念。

广义概念的定量PCR技术是指以外参或内参为标准，通过对PCR终产物的分析或PCR过程的监测，进行PCR起始模板量的定量。

狭义概念的定量PCR技术（严格意义的定量PCR技术）是指用外标法（荧光杂交探针保证特异性）通过监测PCR过程（监测扩增效率）达到精确定量起始模板数的目的，同时以内对照有效排除假阴性结果（扩增效率为零）。

维基百科

Areal-time polymerase chain reaction (Real-Time PCR), also known as quantitativepolymerase chain reaction (qPCR), is a laboratory technique of molecularbiology based on the polymerase chain reaction (PCR). It monitors theamplification of a targeted DNA molecule during the PCR, i.e. in real-time, andnot at its end, as in conventional PCR. Real-time PCR can be usedquantitatively (quantitative real-time PCR), and semi-quantitatively, i.e.above/below a certain amount of DNA molecules (semi quantitative real-timePCR).

优点：灵敏性高，准确性高，通量也还行。

一般而言，RNA-Seq和microassay分析得到的差异表达基因最终也需要通过这种实验方法进行验证。

但是一般适用于验证实验，而不是用于探索性实验。

microarray<small>基因矩阵</small>

基因芯片的概念在上个世纪80年代就已经提出来了， 被评为1998年度自然科学领域十大进展之一。他的基本原理通过设计专门的短核苷酸作为探针，把这些探针固定在专门的基片表面，然后用样本的cDNA进行杂交，根据杂交信号的强弱来判断基因表达的程序。

microarray

但是microarray检测的基因数量完全取决于你的探针设计的数量，而且难以研究mRNA的可变剪切。

RNA-Seq

RNA-Seq是目前基因表达分析最常用的技术。分为以下几步

1.分离所有mRNA

2.逆转录mRNA成cDNA

3.对cDNA测序

4.比对参考基因组

RNA-Seq实验设计中的“重复”包括：技术重复和生物学重复

重复是为了检测组间和组内的变异，对于假设检验至关重要。

(1)技术重复为了估计测量技术（RNA-Seq）的变异。

(2)生物学重复是为了发现生物组内的变异。

简单的说，两组的基因表达的变化只有比组内变异还大时才能认为时显著的。

RNA-Seq的概率分布

相同基因在不同细胞的表达水平服从log-normal(对数正态)分布，由定量PCR验证。

（注：这与相同细胞不同基因表达的分布不同）

但是大多数基因表达实验都是用一群细胞，几乎没有相应分布提出。

RNA-Seq试验中，抽样得到的raw read counts服从泊松分布。并且同一样本在两次试验中的结果不同，这称为shotnoise。这种变异在RNA-Seq技术重复间成为Possion noise。

生物学上不同的样本间的差异服从负二项(negative binomial)分布,有时称gamma-Poisson分布。

由于RNA-Seq count数据也表现出zeroinflation(大量值为0)的特征，所以很难拟合到负二项分布，所以有文章认为要用Poisson-Tweediefamily建模。

RNA-seq数据和microassay在差异表达分析上的区别：

1.RNA-Seq观察到的数据是抽样过程中产生的离散（discrete）count形式。

也就是说总体是恒定的，表达量越高的基因在抽样结果中所占的比例越大。

表达量低的基因可能即便有也无法被检测出来。

当然，重新对相同文库进行测序，还是有可能找到更多表达的转录本

2.microassay检测的是荧光信号的连续度量。

由于使用固定的核酸序列去杂交，所以不是一种“零和游戏”，只要能杂交，就能被检测。

（但如果没有设计相应的引物，就不能检测到可能的基因）

研究意义

1.在不同背景下比较mRNA水平

①同一物种，不同组织：   研究基因在不同部分的表达情况

②同一物种，同一组织：   研究基因在不同处理下，不同条件下的表达变化

③同一组织，不同物种：   研究基因的进化关系

④时间序列实验：         基因在不同时期的表达情况与发育的关系

2.基因分类：       找到细胞特异，疾病相关，处理相关的基因表达模式，用于诊断疾病和预测等

3.基因网络和通路： 基因在细胞活动中的功能，基因间的相互作用。

公共数据库

当然，如果你要研究一个基因的功能时，不要先急着去花钱找公司测序，

先去一些基因表达公共数据库找找看：

1.http://www.ebi.ac.uk/arrayexpress/

2.https://www.ncbi.nlm.nih.gov/geo/

差异表达(differentialexpression,DE)基因分析

通过研究基因的差异表达，我们可以发现

①细胞特异性的基因；

②发育阶段特异性的基因；

③疾病状态相关的基因；

④环境相关的基因；

...

基本方法就是以生物学意义的方式计算基因表达量，然后通过统计学分析表达量寻找具有统计学显著性差异的基因，从而

①选择合适的基因

②衡量结果的可靠性

分析方法

寻找差异表达基因有三种方式：

第一种是计算Fold change（倍数变化），十分简单粗暴的方法，计算方法如下：

E = mean(group1)

B=mean(group2)

FC = (E-B)/min(E,B)

说人话就是，基因A和基因B的平均值之差与两者中较小的比值。选择2-3倍的基因作为结果

（为什么是2-3倍，就是大家约定俗成）。

foldchange

但是简单粗暴的用2到3倍作为阈值，对于低表达的基因，3倍也是噪音，那些高表达的基因，1.1倍都是生物学显著了。更重要的没有考虑到组内变异，没有统计学意义。所以发文章肯定这个图只能作为附录了。

第二种就是统计检验，写文章的时候总需要给出一个p值告诉主编这个结果可信的（虽然p值也存在争论）。

复习一下:p值指的碰巧是拒绝零假设机会。P值越大假阳性越低，同时真实结果也可能会剔除。

注： 基因表达分析的零假设是：基因在不同处理下的表达量相同。

对于基因芯片的数据而言，由于样本服从正态分布，所以可以用t-test（双处理）或anova分析（多处理以上）。

T检验适用于只有两个处理的实验设计，如植物叶片在相同处理第一天和第二天的基因表达差异。

T-Test 实验设计

进行T-test检验时要注意：是双尾检验（存在差异）还是单尾检验（显著性上调或下降），两个样本的总体是不是等方差（标准T检验还是Welch’s test）

如果存在多于两个处理（条件），就需要用到ANOVA分析了。ANOVA分析能主要是研究结果之间的差异是如何引起的，具体请移步到我写方差分析教程。

对于基因表达而言，研究目标是，对于同一个基因而言，他们之间的差异是处理不同造成，还是因为系统误差造成。

单因素

双因素

当然你可以研究，不同基因的表达差异是由因为处理不同，还是基因不同，还是系统误差，还是其中一些的交互作用。

上面都是针对基因芯片的样本服从正态分布进行的统计检验。现在的RNA-Seq，它的抽样过程是离散的，结果是count，服从泊松分布，样本间的差异是服从负二向分布，显然不能按照上述方法分析。

方差分析（ANOVA)和线性回归分析(regression)都是同一时期发展的两套紧密相连的理论。方差分析考量的是离散型自变量(因子）对连续型应变量（响应变量）的模型分析，而线性回归分析只要求响应变量是连续的，对于自变量无要求。如果响应变量不是连续型分布，就要使用更加一般化的广义线性模型（generalized linear model),通过一个连接函数变换响应变量期望，将响应变量的期望与自变量建立线性关系。

因此，我们可以用广义线性模型去分析RNA-seq前期分析得到的离散型结果（count）

方差分析一般用于分析有计划的实验结果，比如说不同处理下的水稻产量。回归分析一般分析没有计划的数据，比如说你可以找到大量体检的数据，只分析其中性别和身高对体重的影响。所以两者各有侧重，不要拿大炮轰蚊子。

统计检验相对于fold change具有统计意义，不需要参考样本，需要处理随机取样。但是需要重复（ANOVA推荐4-10重复），但由于资金和材料等原因，不一定能够满足。此外，对于1000个基因，就要做1000次ANOVA或t-test，最后的p值会有一定的假阳性，因此要做p值矫正（FDR）筛选。

注：推荐在统计检验前过滤表达量低，也就如果一个基因在所有样本中count均低于某一阀值，请在分析前剔除。这个阀值也是约定俗成，一般设置为3.

第三种：Fold Change + 统计检验。说一比较尴尬的事情，在统计检验中你找到越多的差异表达基因，在p值矫正之后，你反而找不到差异表达基因。也就是说，如果在结果中存在大量滥竽充数的所谓的DE基因，那么在严格的p值矫正筛选后，反而会误删真实的DE基因。

因此在p值矫正之前，你先要手动剔除一部分明显就是假阳性的DE基因。这个步骤就需要用到前面的fold-change分析。

我们可以通过火山图来看看如何确定区间

火山图

实战演练

为了方便理解，我们选用同一组数据进行实际操作。默认你懂基本的R语言操作，如安装R包，查看帮助文件等。

正式学习之前，先感谢一下Bioconductor。 根据他们提供流程，我学习到如何进行RNA-Seq分析。

·        http://www.bioconductor.org/help/workflows/RnaSeqGeneEdgeRQL/

·        http://www.bioconductor.org/help/workflows/rnaseqGene/

·        http://www.bioconductor.org/help/workflows/RNAseq123/

·        http://www.bioconductor.org/help/workflows/ExpressionNormalizationWorkflow/

数据介绍和下载

用到的数据来自于一篇文献“A pooled shRNA screen for regulators of primarymammary stem and progenitor cells identifies roles for Asap1 and Prox1”。

实验目的:

Theaim of this study is to determine the absolute and relative expression levelsof mRNA transcripts across multiple flow cytometrically sorted epithelial celltypes including freshly isolated CD24+CD29hi mammary stem cell-enriched basalcells (MaSC/basal), CD24+CD29loCD61+ luminal progenitor-enriched (LP) and theCD24+CD29loCD61- mature luminal-enriched (ML) cell populations. Additionally, acomparison between these primary cell types and cultured MaSC/Basal-derivedmammosphere cells (mammosphere) and the CommaD-βGeo (CommaDβ) cell line wasperformed.

实验设计：

Total RNA was extracted andpurified from sorted luminal or basal populations from the mammary glands offemale virgin 8- to 10-week-old FVB/N mice (3 independent samples forpopulation), MaSC/Basal cells cultured for 1 week under mammsophere conditionsand CommaDβ cells grown under maintenance conditions (Deugnier et al. 2006) andtheir transcriptomes analysed by RNA-Seq. |

根据文章的介绍，数据是100 bp的单端read，用Rsubread比对到mouse reference genome(mm10), 然后使用featureCounts统计每个基因的count数。然后用TMM进行标准化，转换成log2 counts per million.最后用limma包对每个样本每个基因的平均表达值以观察水平权重的线性模型进行拟合，并用T检验找到不同群体的差异表达基因。以FDR + log2-fold-change对基因排序。

我们的任务是下载他们featureCounts的结果，不用limma,改用DESeq2分析。

#install.packages("R.utils")

library("R.utils")

url <-"https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE63310&format=file"

utils::download.file(url,destfile="GSE63310_RAW.tar", mode="wb")

utils::untar("GSE63310_RAW.tar",exdir = ".")

files<- c("GSM1545535_10_6_5_11.txt","GSM1545536_9_6_5_11.txt", "GSM1545538_purep53.txt","GSM1545539_JMS8-2.txt","GSM1545540_JMS8-3.txt", "GSM1545541_JMS8-4.txt",

"GSM1545542_JMS8-5.txt", "GSM1545544_JMS9-P7c.txt","GSM1545545_JMS9-P8c.txt")

for(i inpaste(files, ".gz", sep=""))

R.utils::gunzip(i, overwrite=TRUE)

数据读取

先定一个包含所有文件的向量，方便后续调用。

files <-c("GSM1545535_10_6_5_11.txt", "GSM1545536_9_6_5_11.txt",

"GSM1545538_purep53.txt",  "GSM1545539_JMS8-2.txt",

"GSM1545540_JMS8-3.txt",  "GSM1545541_JMS8-4.txt",

"GSM1545542_JMS8-5.txt",  "GSM1545544_JMS9-P7c.txt",

"GSM1545545_JMS9-P8c.txt")

查看下文件的数据格式,分别记录了每个基因的EntrezID，基因长度和数量

read.delim(files[1],nrow=5)

#  EntrezID  GeneLength  Count

#1   497097       3634     1

#2 100503874      3259     0

#3100038431       1634     0

#4    19888       9747     0

#5    20671       3130     1

之后我们需要用DESeqDataSetFromMatrix为DESeq2提供数据,命令形式如下：

library("DEseq2")

DESeqDataSetFromMatrix(countData,colData, design, tidy = FALSE,ignoreRank = FALSE, ...)

其中countData存放counts数据，colData存放样本信息的数据，design就是实验设计。

首先从各个count matrix文件中读取count，基因长度部分可以舍弃，因为DESeq2只需要为标准化的count数据，不需要提供基因长度信息。

逻辑就是分别读取每一个文件的count列，然后赋予文件名。

sample1<- read.delim(files[1],header = T)[,c(1,3)]

count.table<- data.frame(sample1)

for( f in files[2:length(files)]){

column <- read.delim(f,header=T,row.names = 1)[,2]

count.table <- cbind(count.table, column)

}

head(count.table)

rownames(count.table)<- count.table[,1]

count.table<- count.table[-1]

samplenames<- substring(files, 12, nchar(files)-4)

colnames(count.table)<- samplenames

然后要提供colData,其中colData存放列名要和countData的行名相同。

#colData

condition<- as.factor(c("LP", "ML", "Basal","Basal", "ML", "LP", "Basal","ML", "LP"))

lane<- as.factor(rep(c("L004","L006","L008"),c(3,4,2)))

colData<- data.frame(condition = condition, lane=lane,row.names = samplenames)

all(rownames(colData)%in% colnames(count.table))

最后导入数据:

dds<- DESeqDataSetFromMatrix(countData = count.table,

colData = colData,

design = ~ condition )

数据预处理

预过滤

虽然DESeq2会自动屏蔽那些低count的基因，但是剔除那些几乎不存在基因的部分能够提高运行速度。

dds<- dds [ rowSums(counts(dds)) > 1, ]

rlog和方差齐性转换

许多常见的多维数据探索性分析的统计分析方法，例如聚类和主成分分析要求，在那些同方差性的数据表现良好。所谓的同方差性就是虽然平均值不同，但是方差相同。

但是对于RNA-Seq count数据而言，当均值增加时，方差期望也会提高。也就说直接对count matrix或标准化count（根据测序深度调整）做PCA分析，由于高count在不同样本间的绝对差值大，也就会对结果有很大影响。简单粗暴的方法就是对count matrix取log后加1。这个1也是约定俗成，看经验了。

随便举个栗子看下效果：

lambda<- 10^seq(from = -1, to = 2, length = 1000)

cts<- matrix(rpois(1000*100, lambda), ncol = 100)

library("vsn")

meanSdPlot(cts,ranks = FALSE)

平滑前

log.cts.one<- log2(cts + 1)

meanSdPlot(log.cts.one,ranks = FALSE)

平滑后

DESeq2为count数据提供了两类变换方法，使得不同均值的方差趋于稳定：

regularized-logarithmtransformation or rlog(Love, Huber, and Anders 2014)和

variancestabilizing transformation(VST)   (Anders and Huber2010)

用于处理含有色散平均趋势负二项数据。

到底用啥：数据集小于30 -> rlog，大数据集 -> VST。还有这个处理过程不是用于差异检验的，在DESeq分析中会自动选择最合适的所以你更不需要纠结了，记得用raw count。

rld<- rlog(dds, blind = FALSE)

head(assay(rld),3)

vsd<- vst(dds, blind = FALSE)

head(assay(vsd),3)

处理的结果如下：

library("dplyr")

library("ggplot2")

dds<- estimateSizeFactors(dds)

df<- bind_rows(

as_data_frame(log2(counts(dds, normalized=TRUE)[, 1:2]+1)) %>%

mutate(transformation = "log2(x + 1)"),

as_data_frame(assay(rld)[, 1:2]) %>% mutate(transformation ="rlog"),

as_data_frame(assay(vsd)[, 1:2]) %>% mutate(transformation ="vst"))

colnames(df)[1:2]<- c("x", "y")

ggplot(df,aes(x = x, y = y)) + geom_hex(bins = 80) +

coord_fixed() + facet_grid( . ~ transformation)

不同转换的分布

结果就是转换后更加集中了。

样本距离

RNA-Seq分析第一步通常是评估样本间的总体相似度。

①那些样本最接近

②那些样本差异较大

③这与实验设计预期符合么

这里使用R内置的 dist 计算 rlog变换数据的Euclidean distance，然后用pheatmap根据结果画热图

sampleDists<- dist(t(assay(rld)))

library("pheatmap")

library("RColorBrewer")

sampleDistMatrix<- as.matrix( sampleDists )

colors<- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)

pheatmap(sampleDistMatrix,

clustering_distance_rows = sampleDists,

clustering_distance_cols = sampleDists,

col = colors)

heatmap

同样的可以用Poisson Distance (Witten 2011)计算距离，计算方式如下：

library("PoiClaClu")

poisd<- PoissonDistance(t(counts(dds)))

PCA分析

还有一种可视化样本-样本距离的方法就是主成分分析。PCA分析我打算找点资料好好理解之后再写，这个说下有这个方法。

DESeq2提供了专门的方法用于作图，还很好看呢！

plotPCA(rld,intgroup=c("condition"))

PCA

能够明显的发现不同处理的距离离得很远。

差异表达基因分析（DEA）

在我们定义好DESeqDataSe 就可以方便的调用DESeq分析.

dds<- DESeq(dds)

在几行输出后信息后，分析就完成了，更多具体参数可以用?DESeq查看手册

构建结果表格

如果直接调用results，只会提取出design matrix最后两个变量（这里是ML vs Basal）的

log2fold changes and p values等结果。

>results(dds)

log2fold change (MAP): condition ML vs Basal

Waldtest p-value: condition ML vs Basal

DataFramewith 27179 rows and 6 columns

baseMean log2FoldChange    lfcSE       stat      pvalue         padj

<numeric>      <numeric> <numeric>  <numeric>    <numeric>   <numeric>

497097   118.7522256     -7.0662875 0.5802473 -12.17806084.068401e-34 8.556128e-33

100503874  1.3308906     -2.6371037 1.3370584  -1.97231744.857338e-02 8.946685e-02

100038431  0.0843771     -0.1541107 1.2443948  -0.12384399.014388e-01           NA

19888      1.4833878      2.9503142 1.3351543  2.2097179 2.712475e-02 5.385925e-02

20671     26.4317752     -1.5256766 0.7754853  -1.96738294.913908e-02 9.034136e-02

...              ...           ...       ...        ...         ...          ...

100861837  368.84237     0.33334299 0.3226665  1.0330883    0.3015626    0.4139846

100861924   22.79272    -1.25169702 0.8307249 -1.5067528    0.1318740    0.2107990

170942     852.61538    -0.07612745 0.2399318 -0.3172879    0.7510252    0.8235139

100861691    0.00000            NA       NA         NA          NA           NA

100504472    0.00000            NA       NA         NA          NA           NA

实际上results有如下众多参数

results(object,contrast, name, lfcThreshold = 0,

altHypothesis = c("greaterAbs", "lessAbs","greater", "less"),

listValues = c(1, -1), cooksCutoff, independentFiltering = TRUE,

alpha = 0.1, filter, theta, pAdjustMethod = "BH", filterFun,

format = c("DataFrame", "GRanges", "GRangesList"),test, addMLE = FALSE,

tidy = FALSE, parallel = FALSE, BPPARAM = bpparam(), ...)

比如可以指定比较对象Basal和PL,可以用mcols查看结果存储的元数据，了解列名的含义。

res<- results(dds,contrast =c("condition","Basal","LP"))

mcols(res,use.names = TRUE)

DataFramewith 6 rows and 2 columns

type                                  description

<character>                                  <character>

baseMean      intermediate     mean of normalized counts for all samples

log2FoldChange     results log2 fold change (MAP): condition Basal vs LP

lfcSE              results         standard error:condition Basal vs LP

stat               results         Wald statistic:condition Basal vs LP

pvalue             results      Wald test p-value: condition Basal vs LP

padj               results                         BH adjusted p-valuess

其中lfcSE是Log2 fold change的标准误。其他项都比较好理解。

最后还可以看一些总结性的内容

summary(res)

outof 22026 with nonzero total read count

adjustedp-value < 0.1

LFC> 0 (up)     : 5528, 25%

LFC< 0 (down)   : 5274, 24%

outliers[1]     : 43, 0.2%

lowcounts [2]   : 2953, 13%

(meancount < 1)

[1]see 'cooksCutoff' argument of ?results

[2]see 'independentFiltering' argument of ?results

在limma分析结果分别是4127，4298，稍微多了那么2000个基因。其他结果也基本上都了2000个。

看起来对于22000多个基因而言，差异不算太大。

当然我们还可以通过

①降低假阳性率（FDR， false discovery rate)阈值

②提高log2 fold change的阈值

>res0.5 <- results(dds, contrast =c("condition","Basal","LP"),alpha=0.05)

>table(res0.5$padj < 0.05)

FALSE TRUE

8860 9750

>resLFC1 <- results(dds, lcontrast =c("condition","Basal","LP"),fcThreshold=1)

>table(resLFC1$padj < 0.1)

FALSE TRUE

891610958

于是乎结果就和limma差不多了。

多重试验

在高通量数据分析中，我们通常不是用p值来拒绝原假设，更多是用来进行多重试验矫正。

为什么要做多重试验矫正？

如果对一个基因，我有99%的把握认为是差异表达的，也就是说1%的可能是错的。那么假设有10000个基因，按照数学期望，有100个是假的。因此为了保证多重试验结果的可靠性要对结果的p-value做矫正。

DESeq2用的Benjamini-Hochberg(BH) adjustment (Benjamini and Hochberg 1995)，计算矫正后的P-value,用于回答如下问题：

如果我们选择了所有小于或等于矫正p-value阈值的所有显著性基因，假阳性比例（ false discovery rate, FDR）是多少？

FDR在高通量试验中是比较有用的统计值，由于我们通常关注一类自己感兴趣的基因，所以我们需要设置一个假阳性上限。

我们从结果中以1%作为显著性，分别找出显著性上调和下调的基因。

regSig<- subset(res, padj < 0.01)

#Up

head(resSig[order(resSig$log2FoldChange), ])

#Down

head(resSig[order(resSig$log2FoldChange, decreasing = TRUE), ])

可视化

人类对图像比较敏感，对文字比较迟钝。所以我们需要一些比较好看的图来放到文章中解释说明。

最简单的就是Counts plot，看看特定基因的count数量。

topGene<- rownames(res)[which.min(res$padj)]

plotCounts(dds,gene = topGene, intgroup=c("condition"))

count plot

可以用MA-plot（也叫mean-differenceplot,Bland-Altman plot）了解模型（如所有基因在不同处理比较的结果）的估计系数的分布

res<- lfcShrink(dds, contrast=c("contion","ML","PL"),res=res)

plotMA(res,ylim = c(-5, 5))

更加喜闻乐见的是基因聚类所提供的热图展示。我们可以找前20个样本件差异比较大，然后看他们在不同样本间的表达情况。

library("genefilter")

topVarGenes<- head(order(rowVars(assay(rld)), decreasing = TRUE), 20)

mat <- assay(rld)[ topVarGenes, ]

mat <- mat - rowMeans(mat)

anno<- as.data.frame(colData(rld)[, c("condition","lane")])

pheatmap(mat,annotation_col= anno)

heatmap

大致可以发现同一组的基因颜色是相同的，也就是说表达量相近。

结语

找到差异表达的基因只是第一步，后续还需要对这些基因进行进一步的分析，如

①基因富集分析

②主成分分析

③聚类分析

而这些内容就是我接下来的学习重点，也是下次更新文章的主题。

基因表达分析（中）-富集分析

回顾

首先对基因表达分析（上）做一个简单的回顾

1.研究基因表达的有如下工具：RNA-Seq，microarray， qRT-PCR等（欢迎补充）

2.RNA-Seq，microarray一般用在探索性阶段，qRT-PCR用于验证

3.RNA-Seq和microarray由于他们的实验方式不同，导致寻找差异表达基因的统计学方法也不同。其中microarray使用寡核苷酸作为探针进行杂交，基因表达量与亮度正相关，而亮度是一个连续型变量，因此大多认为结果是服从正态分布。而RNA-Seq的测序结果是一条条read，是一种离散抽样过程，因此认为是服从泊松分布。

4.ANOVA和简单线性模型都是广义线性模型的特殊情况。ANOVA是研究名义型解释变量和连续型解释变量的关系，简单线性模式是研究连续型解释变量和连续型解释变量的关系。而广义线性模式没特殊要求。

5.在3,4的背景下，microarray一般用t检验（两个条件），ANOVA分析（多个条件），最常用limma（线性模型）进行检验。RNA-Seq有许多基于count的R包，如DESeq，DESeq2,（基于负二向分布广义线性模型）

6.以上要求你每个条件都要有3个重复（目前投稿要求），你要是老板穷，一个重复都不给，那你去Google解决方案吧。

7.用R作差异表达分析大致分为以下几步：

1）根据软件包要求导入数据；

2）数据预处理，把那些只有0或1计数结果的基因去掉，提高效率。

这一步还可以进行探索性数据分析；

3）跑程序，得到结果；

4）对结果进行可视化，看看基因聚类等结果，这一步不是必须的，但却是展示数据最好的手段了。

如果这些内容还有印象,那么我们继续上次没有写完的内容的其中一个部分-基因富集分析

为何要基因富集分析

在基因差异表达分析之后,你得到了好多p值特别小（也就是显著性很高）的基因，那么下一步你想做什么？

①选择一些基因用于验证？

②对其中基因进行后续研究？

③在结果中把这些基因都放在后面？

④尝试着把所有基因相关的文献都都读读看（劝你放弃这个念头）？

欢迎补充

这些想法都是非常顺理成章的，但是不要着急。

首先，差异表达找到的基因往往很多，你简单的粗暴去找每一个基因的详细资料，显然不太现实；

其次，如果我们单纯觉得某一个基因和你研究的课题相关，或者说你其实已经找到了一个有可能的基因（或者你只是希望用一些高大上的实验验证一下）那么这个行为是不是有太多主观性，存在一些偏见。

当然，你觉得基因就是你要找的，可是万一它只是碰巧来打酱油的呢，这不就是很尴尬了。

所以为了让审稿人相信你的结果，你就需要做一个基因富集分析哦。

什么是基因富集分析

基因富集分析(geneset enrichment analysis)是在一组基因或蛋白中找到一类过表达的基因或蛋白。一般是高通量实验，如基因芯片，RNA-Seq，蛋白质组学（质谱结果）的后续步骤。

基因富集分析需要我们提供某一类功能基因的集合用于背景，常用的注释数据库如：

①TheGene Ontology Consortium:                描述基因的层级关系

②KyotoEncyclopedia of Genes and Genomes:      提供了pathway的数据库。

分析方法

在文献TenYears of Pathway Analysis: Current Approaches and Outstanding Challenges（推荐大家看一遍）作者将研究方法归为三种：

其中第三种方法想的很好就是难度很大。而且贴心的把每一种方法有哪些工具都总结出来了：

Over-Repressentation Analysis(ORA)

ORA是目前商业化最多的方法。为了说明他的基本思想，我要举一个喜闻乐见的例子：读书无用论。

这是我百度找到搜狐财经一篇文章《大数据告诉你真正的有钱人是什么样》的有钱人的学历分布情况，高学历人群（本科以上，因为本科生太多了）所占的比例是9.4%，其他都是一般学历占90.6%。这时候，有些公众号就可以开始不带脑子的说了，读书没什么用呀，有钱人中都是一般学历的呀，以后读书读到大学就行了，甚至也可以不上本科呀（34.2%本科和本科以下）。

你每年回家总能回去看到有人炫耀说，虽然我有钱，可是读书太少了，都不能和你们读书人比的。你总感觉哪里不对劲，但是却又不太方便说出来。

实际上，这就是因为没有考虑到背景。因为高学历本身人数就不多，当然在有钱人里面的人数也就相应不多了。我们要证明有钱人更多是富集高学历这一部分。

类别

有钱人

整体

高学历

10

50

一般学历

90

950

100

1000

H0: 是否有钱和学历高无关

Ha: 学历高还是有点用的

然后做一个Fisher精确检验，看看p值。

richer.pop<- matrix(data = c(10,90,50,950),nrow=2)

fisher.test(richer.pop,alternative = "greater")

Fisher's Exact Test for Count Data

data: richer.pop

p-value= 0.03857

alternativehypothesis: true odds ratio is greater than 1

95percent confidence interval:

1.052584     Inf

sampleestimates:

oddsratio

2.109244

p值小于0.05，看来我读个博士让我以后有钱概率变大了。

现在将我们上面的有钱人改成我们找到的基因，整体改成所有基因。高学历表示属于目标注释基因集，一般学历就是非注释基因组.我们就是要判断我们找到的基因更多是在目标注释集中。

所以你需要列出下表，然后再做一个fisher.test()。

类别

gene list

Genome

in anno group

10

50

not in anno group

290

19950

300

20000

上述的基本思想就是统计学的白球黑球实验：

在一个黑箱里，有确定数量的黑白两种球，你随机抽取（不放回）M个球中,其中两种球的比例分别是多少？

除了用Fisher精确检验，还有其他统计方法：

①Hypergeometric (fisher精确检验用的就是超几何检验)http://www.bio-info-trainee.com/1225.html

②Binomial: 二项分布要求是有放回，无放回要求整体足够大大到可以近似。

③Chi-squared chisq.test(counts)

④Z

⑤Kolmogorov-Smirnov

⑥Permutation http://www.bio-info-trainee.com/1237.html

ORA的方法就是如此的简单，但是有一个问题，就是你如何确定哪些基因是差异表达的，你还是需要设置一个人为的cutoff, 主观能动性成分有点大。

Functional Class Scoring(FCS)

FCS认为，“虽然个体基因表达改变之后会更多在通路中体现，但是一些功能相关基因中较弱但协调的变化也有明显的影响。”

Thehypothesis of functional class scoring (FCS) is that although large changes inindividual genes can have significant effects on pathways, weaker butcoordinated changes in sets of functionally related genes (i.e., pathways) canalso have significant effects

FCS分析方法稍微复杂了一点，他要求的输入是一个排序的基因列表和一个基因集合。MIT， Broad Institute 2007年文献就提供了这一方法的软件"GSEA"

the install screen for GSEA

有如下特点：

①计算所有输入基因集合的分数，而不是单个基因

②不需要设置cutoff

③找到一组相关的基因

④提供了更加稳健的统计框架

GSEA是一款图形化的软件，根据他们提供的教程，然后点呀点，就会得到如下结果。

下图就是需要好好理解的部分。

GESA

中间从蓝色到红色的过渡“带”表示基因从上调到下调排列（排序可以按照fold change,也可以是p-value)。黑色像条形码的竖线表示该位置的基因属于某个指定通路。绿色有波动的曲线表示富集分数，从0开始计算，属于基因通路增加，不属于则减少。最后看下黑色的条形码是不是富集在一端。

那如何做统计检验呢？

Thefinal step in FCS is assessing the statistical significance of thepathway-level statistic. When computing statistical significance, the nullhypothesis tested by current pathway analysis approaches can be broadly dividedinto two categories: i) competitive null hypothesis and ii) self-contained nullhypothesis [3], [18], [22], [31]. A self-contained null hypothesis permutesclass labels (i.e., phenotypes) for each sample and compares the set of genesin a given pathway with itself, while ignoring the genes that are not in thepathway. On the other hand, a competitive null hypothesis permutes gene labelsfor each pathway, and compares the set of genes in the pathway with a set ofgenes that are not in the pathway。

我们要检验的目标是基因富集在一端是因为于目标通路相关的基因都在一端富集。那么空假设就是，你把找到的基因随便摆放也能看到富集现象。用比较专业的话说就是先生成一个零假设的数据分布，然后观察实际数据在这个零假设分布下，是不是在尾端。

一些问题

统计检验的能力是有限的，所以还有很多问题存在解决。

①我们希望找到生物学显著的基因，但是生物学显著和统计显著两者并不是完全相关

②无论是ORA还是FCS都对背景（也就是这个物种一共有多少基因）有要求，但是随着我们的研究深入，基因数量会改变。

有些软件会直接设置一个很大的背景数，从而让p值很显著，然后我们就开心地用他们的结果。

③有些基因没有注释，也就是注释缺失，处理方法就是扔（欢迎拍砖）。

④有一些注释项是其他项的子集。

实战：用clusterProfiler做富集分析

clusterProfiler是Y叔良心之作（当然他的良心之作还有很多），目前支持KEGG在线拉数据，支持DAVID，支持BroadiNSTITUTE Molecular Signatures Databases， 支持GSEA。你问我支不支持，我当然是支持的啦。所以这里放上Y叔的公众号了。

安装：

source("https://bioconductor.org/biocLite.R")

biocLite('clusterProfiler')

这里简单举例如何使用，更多内容见https://guangchuangyu.github.io/clusterProfiler/,

这个部分你可以直接过掉，因为我只是跟着Y叔的代码敲了一遍而已。

加载差异表达基因，我这里偷懒就随机挑一些基因名（来自于DOSE包，Disease Ontology Semantic and Enrichment analysis ）出来了。

library(org.Hs.eg.db)

data(geneList)

gene <-names(geneList)[abs(geneList) > 2]

head(gene)

[1]"4312"  "8318"  "10874""55143" "55388" "991"

Y叔为了展示他能够处理不同命名方式的ID，用bitr（来自于clusterProfiler)进行生物学ID转换

gene.df<- bitr(gene, fromType = "ENTREZID",

toType = c("ENSEMBL", "SYMBOL"),

OrgDb = org.Hs.eg.db)

head(gene.df)

ENTREZID         ENSEMBL SYMBOL

1    4312 ENSG00000196611   MMP1

2    8318 ENSG00000093009  CDC45

3   10874 ENSG00000109255    NMU

4   55143 ENSG00000134690  CDCA8

5   55388 ENSG00000065328  MCM10

6     991 ENSG00000117399  CDC20

clusterProfiler: ORA

然后对不同命名的ID都做ORA富集分析。

ego<- enrichGO(gene          =gene,

universe      = names(geneList),

OrgDb         = org.Hs.eg.db,

ont           ="CC",

pAdjustMethod = "BH",

pvalueCutoff  = 0.01,

qvalueCutoff  = 0.05)

ego2<- enrichGO(gene         =gene.df$ENSEMBL,

OrgDb         = org.Hs.eg.db,

keytype       = 'ENSEMBL',

ont           ="CC",

pAdjustMethod = "BH",

pvalueCutoff  = 0.01,

qvalueCutoff  = 0.05)

ego3<- enrichGO(gene         =gene.df$SYMBOL,

OrgDb         = org.Hs.eg.db,

keytype       = 'SYMBOL',

ont           ="CC",

pAdjustMethod = "BH",

pvalueCutoff  = 0.01,

qvalueCutoff  = 0.05)

结果会有如下

ID：          基因本体论的ID

Description： 描述

GeneRatio：   在GO词条所占的比例

BgRatio ：    在背景所占的比例

pvalue：      假设是正确但是被拒绝的概率

p.adjust：    采用BH方法进行多重试验p值校正

qvalue：      Q值=被拒绝但却是正确的概率

geneID：     列出在这个GO词条下的我们提供的基因

count：      数量，如果没有约分，就是GeneRatio的分子了。

所以从GeneRatio：24/199，BgRatio：231/11711可以反推出下表：

类别

gene list

Genome

in anno group

24

199

not in anno group

231

11711

255

11910

最后再画一张美美的点图,根据GeneRatio所做。

dotplot(ego,showCategory=30)

clusterProfiler: GSEC

clusterProfiler支持GSEC，而且很好用。

gsecc<- gseGO(geneList=geneList, ont="CC", OrgDb=org.Hs.eg.db,verbose=F)

head(summary(gsecc)

gseaplot(gsecc,geneSetID="GO:0000779"