差异表达分析之FDR

差异表达分析之FDR
随着测序成本的不断降低，转录组测序分析已逐渐成为一种很常用的分析手段。但对于转录组分析当中的一些概念，很多人还不是很清楚。今天，小编就来谈谈在转录组分析中，经常会遇到的一个概念FDR，那什么是FDR?为什么要用FDR呢?一起来学习吧！

什么是FDR

FDR (false discovery rate)，中文一般译作错误发现率。在转录组分析中，主要用在差异表达基因的分析中，控制最终分析结果中，假阳性结果的比例。

为什么要用FDR

在转录组分析中，如何确定某个转录本在不同的样品中表达量是否有差异是分析的核心内容之一。一般来说，我们认为，不同样品中，表达量差异在两倍以上的转录本，是具有表达差异的转录本。为了判断两个样品之间的表达量差异究竟是由于各种误差导致的还是本质差异，我们需要根据所有基因在这两个样本中的表达量数据进行假设检验。常用的假设检验方法有t-检验、卡方检验等。很多刚接触转录组分析的人可能会有这样一个疑问，一个转录本是不是差异表达，做完假设检验看P-value不就可以了么？为什么会有FDR这样一个新的概念出现？这是因为转录组分析并不是针对一个或几个转录本进行分析，转录组分析的是一个样品中所转录表达的所有转录本。所以，一个样品当中有多少转录本，就需要对多少转录本进行假设检验。这会导致一个很严重的问题，在单次假设检验中较低的假阳性比例会累积到一个非常惊人的程度。举个不太严谨的例子。

假设现在有这样一个项目：

● 包含两个样品，共得到10000条转录本的表达量数据，

● 其中有100条转录本的表达量在两个样品中是有差异的。

● 针对单个基因的差异表达分析有1%的假阳性。

由于存在1%假阳性的结果，在我们分析完这10000个基因后，我们会得到100个假阳性导致的错误结果，加上100条真实存在的结果，共计200个结果。在这个例子中，一次分析得到的200个差异表达基因中，有50%都是假阳性导致的错误结果，这显然是不可接受的。为了解决这个问题，FDR这个概念被引入，以控制最终得到的分析结果中假阳性的比例。

如何计算FDR

FDR的计算是根据假设检验的P-value进行校正而得到的。一般来说，FDR的计算采用Benjamini-Hochberg方法（简称BH法），计算方法如下：

1. 将所有P-value升序排列.P-value记为P，P-value的序号记为i，P-value的总数记为m

2. FDR(i)=P(i)*m/i

3. 根据i的取值从大到小，依次执行FDR(i)=min{FDR(i),FDR(i+1)}

注：实际上，BH法的原始算法是找到一个最大的i，满足P≤i/m*FDR阈值，此时，所有小于i的数据就都可以认为是显著的。在实践中，为了能够在比较方便的用不同的FDR阈值对数据进行分析，采用了步骤3里的方法。这个方法可以保证，不论FDR阈值选择多少，都可以直接根据FDR的数值来直接找到所有显著的数据。

下面我们以一个包含10个数据的例子来看一下FDR计算的过程差异表达分析之FDR

在这个例子中，第一列是原始的P-value，第二列是排序后的序号，第三列是根据P-value校正得到的初始FDR，第四列是最终用于筛选数据的FDR数值。如果我们设定FDR<0.05，那么绿色高亮的两个数据就是最终分析认为显著的数据。

FDR的阈值选择在转录组分析中是非常重要的一个环节，常用的阈值包括0.01、0.05、0.1等。实践中也可以根据实际的需要来灵活选择。例如，在做真菌或者原核生物的转录组分析时，由于这些物种转录本数量较少，假阳性累积的程度较低，所以可以适当将FDR阈值设置的较高一些，这样可以获得较多的差异表达结果，有利于后续的分析。

浅谈多重检验校正FDR
Posted: 四月 12, 2017 Under: Basic By Kai no Comments

例如，在我们对鉴定到的差异蛋白做GO功能注释后，通常会计算一个p值。当某个蛋白的p值小于0.05（5%）时，我们通常认为这个蛋白在两个样本中的表达是有差异的。但是仍旧有5%的概率，这个蛋白并不是差异蛋白。那么我们就错误地否认了原假设（在两个样本中没有差异表达），导致了假阳性的产生（犯错的概率为5%）。

如果检验一次，犯错的概率是5%；检测10000次，犯错的次数就是500次，即额外多出了500次差异的结论（即使实际没有差异）。 为了控制假阳性的次数，于是我们需要对p值进行多重检验校正，提高阈值。

第一种方法Bonferroni，最简单严厉的方法。
例如，如果检验1000次，我们就讲阈值设定为5% / 1000 = 0.00005；即使检验1000次，犯错误的概率还是保持在N×1000 = 5%。最终使得预期犯错误的次数不到1次，抹杀了一切假阳性的概率。但是该方法虽然简单，但是检验过于严格，导致最后找不到显著表达的蛋白（假阴性）。

第二种方法FDR（False Discovery Rate）
相对Bonferroni来说，FDR用比较温和的方法对p值进行了校正。其试图在假阳性和假阴性间达到平衡，将假/真阳性比例控制到一定范围之内。例如，如果检验1000次，我们设定的阈值为0.05（5%），那么无论我们得到多少个差异蛋白，这些差异蛋白出现假阳性的概率保持在5%之内，这就叫FDR＜5%。

那么我们怎么从p value 来估算FDR呢，人们设计了几种不同的估算模型。其中使用最多的是Benjamini and Hochberg方法，简称BH法。虽然这个估算公式并不够完美，但是也能解决大部分的问题，主要还是简单好用！

FDR的计算方法
除了可以使用excel的BH计算方法外，对于较大的数据，我们推荐使用R命令p.adjust。 p.adjust(p, method = p.adjust.methods, n = length§)

p.adjust.methods

c(“holm”, “hochberg”, “hommel”, “bonferroni”, “BH”, “BY”,

“fdr”, “none”)

我们还可以从R命令p.adjust的源代码，了解其运行的机制是什么。

p.adjust
function (p, method = p.adjust.methods, n = length§){
method <- match.arg(method)
if (method == “fdr”)
method <- “BH”
nm <- names§
p <- as.numeric§
……
BH = {
i <- lp:1L
o <- order(p, decreasing = TRUE)
ro <- order(o)
pmin(1, cummin(n/i * p[o]))[ro]
}
……
p0
}
其实该函数表达的意思是这样的：

我们将一系列p值、校正方法（BH）以及所有p值的个数（length§）输入到p.adjust函数中。
将一系列的p值按照从大到小排序，然后利用下述公式计算每个p值所对应的FDR值。 公式：p * (n/i)， p是这一次检验的p value，n是检验的次数，i是排序后的位置ID（如最大的P值的i值肯定为n，第二大则是n-1，依次至最小为1）。
将计算出来的FDR值赋予给排序后的p值，如果某一个p值所对应的FDR值大于前一位p值（排序的前一位）所对应的FDR值，则放弃公式计算出来的FDR值，选用与它前一位相同的值。因此会产生连续相同FDR值的现象；反之则保留计算的FDR值。
将FDR值按照最初始的p值的顺序进行重新排序，返回结果。

p值还是 FDR ?
差异分析
如何筛选显著性差异基因，p value, FDR 如何选
经常有同学询问如何筛选差异的基因（蛋白）。已经计算了表达量和p value值，差异的基因（蛋白）太多了，如何筛选。其中最为关键的是需要对p value进行校正。

基本概念：

零假设：在随机条件下的分布。

p值：在零假设下，观测到某一特定实验结果的概率称为p值。

假阳性：得到了阳性结果，但这个阳性结果是假的。

假阴性：得到了阴性结果，但这个阴性结果是假的。

单次检验：

针对单个基因（蛋白），采用统计检验，假设采用的p值为小于0.05，我们通常认为这个基因在两个（组）样本中的表达是有显著差异的，但是仍旧有5%的概率，这个基因并不是差异基因。

单多次检验：

当两个（组）样本中有10000个基因采用同样的检验方式进行统计检验时，这个时候就有一个问题，单次犯错的概率为0.05， 进行10000次检验的话，那么就有0.05*10000=500 个基因的差异被错误估计了。

多重检验矫正：

为了解决多次检验带来的问题，我们需要对多次检验进行校正。那如何校正呢？在此介绍两种方法：

Bonferroni 校正法
Bonferroni校正法：如果进行N次检验，那么p值的筛选的阈值设定为p/N。 比如，进行10000次检验的话，如果p值选择为0.05， 那么校正的p值筛选为0.000005。 p值低于此的基因才是显著性差异基因。
该方法虽然简单，但是过于严格，导致最后找的差异基因很少，甚至找不到差异的基因。

FDR(False Discovery Rate) 校正法
FDR错误控制法是Benjamini于1995年提出的一种方法，基本原理是通过控制FDR值来决定p值的值域。相对Bonferroni来说，FDR用比较温和的方法对p值进行了校正。其试图在假阳性和假阴性间达到平衡，将假/真阳性比例控制到一定范围之内。
那么怎么从p值来估算FDR呢，人们设计了几种不同的估算模型。其中使用最多的是Benjamini and Hochberg方法，简称BH法。该方法分两步完成，具体如下：
2.1 假设总共有m个候选基因，每个基因对应的p值从小到大排列分别是p(1),p(2),…,p(m)
2.2 若想控制FDR不能超过q，则只需找到最大的正整数i，使得 p(i)<= (i*q)/m . 然后，挑选对应p(1),p(2),…,p(i)的基因做为差异表达基因，这样就能从统计学上保证FDR不超过q。

如何实现多重检验：

如果你了解R语言的话，那么采用p.adjust方法就可以了。

简单使用DESeq做差异分析
Posted: 五月 06, 2017 Under: Transcriptomics By Kai no Comments

DESeq这个R包主要针对count data，其数据来源可以是RNA-Seq或者其他高通量测序数据。类似地，对于CHIP-Seq数据或者质谱肽段数据也是使用的。

由于DESeq是一个R包，因此使用它需要一点点R基础语法。

首先需要读入一个数据框，列代表每个sample，行代表每个gene

database_all <- read.table(file = “readcount”, sep = “\t”, header = T, row.names = 1)
database <- database_all[,1:6]
这里主要对于两两比较的数据，因此我取了数据的前6列，分别是两组样品，每组3个生物学重复

设定分组信息，也就是样本分组的名称

type <- factor(c(rep(“LC_1”,3), rep(“LC_2”,3)))
我这里是样品1是LC_1，样品2是LC_2

由于DESeq包要求接下来的count data必须要整数型，因此我们需要对数据进行取整,然后将数据database和分组信息type读入到cds对象中

database <- round(as.matrix(database))
cds <- newCountDataSet(database,type)
接下里对于不同类型的数据要进行不同的处理，可以粗略分为有生物学重复数据、有部分生物学重复数据以及无生物学重复数据

4.1 有生物学重复

cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds)
res <- nbinomTest(cds,“LC_1”,“LC_2”)
其通过estimate the dispersion并对count data进行标准化，然后得到每个gene做T test检验

4.2 对于部分样品有生物学重复

cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds)
res <- nbinomTest(cds,“LC_1”,“LC_2”)
其步骤跟有上述的一样的，DESeq会根据有生物学重复的样品来estimate the dispersion，当然要保证unreplicated condition does not have larger variation than the replicated one

4.3 对于没有生物学重复

cds <- estimateDispersions(cds, method=“blind”, sharingMode=“fit-only” )
res <- nbinomTest(cds,“LC_1”,“LC_2”)
注意参数method=”blind” 和 sharingMode=”fit-only”即可

最后就是查看符合阈值的差异基因有多少个即可，然后将结果输出到csv文件中方便查看

table(res p a d j < 0.05 ) r e s < − r e s [ o r d e r ( r e s padj <0.05) res <- res[order(res padj<0.05)res<−res[order(respadj),]
sum(res$padj<=0.01,na.rm = T)
write.csv(resdata,file = “LC_1_vs_LC_2_DESeq.csv”)
Summary
DESeq在前几年的文章中经常被使用，但是现在有了其升级版DESeq2，后者相比前者对于犯第一类错误卡的并不是那么严格了，所以在同样的padj的阈值下，筛选到的差异基因的数目也会多一点。

并且上述中，我都只使用nbinomTest来获得由显著差异的gene，其实DESeq包还提供了其他检验方法，比如nbinomGLMTest函数，具体可以查看DESeq文档，还有其他对应对因素分析的使用方法。

DESeq2和EdgeR都可用于做基因差异表达分析，主要也是用于RNA-Seq数据，同样也可以处理类似的ChIP-Seq,shRNA以及质谱数据。

这两个都属于R包，其相同点在于都是对count data数据进行处理，都是基于负二项分布模型。因此会发现，用两者处理同一组数据，最后在相同阈值下筛选出的大部分基因都是一样的，但是有一部分不同应该是由于其估计离散度的不同方法所导致的。

DESeq2的使用方法：
输入矩阵数据，行名为sample，列名为gene；DESeq2不支持无生物学重复的数据，因此我选择了2个样本，3个生物学重复的数据；并对count data取整（经大神指点，这里需要说明下，我的测试数据readcount是RSEM定量的结果，并不是常见的htseq-count的结果，所以count值会有小数点，而DESeq2包不支持count数有小数点，所以这里需要round取整）。

database_all <- read.table(file = “readcount”, sep = “\t”, header = T, row.names = 1)
database <- database_all[,1:6]
#type <- factor(c(rep(“LC_1”,3), rep(“LC_2”,3)))
database <- round(as.matrix(database))
设置分组信息以及构建dds对象

condition <- factor(c(rep(“LC_1”,3), rep(“LC_2”,3)))
coldata <- data.frame(row.names = colnames(database), condition)
dds <- DESeqDataSetFromMatrix(countData=database, colData=coldata, design=~condition)
使用DESeq函数进行估计离散度，然后进行标准的差异表达分析，得到res对象结果

dds <- DESeq(dds)
res <- results(dds)
最后设定阈值，筛选差异基因，导出数据

table(res p a d j < 0.05 ) r e s < − r e s [ o r d e r ( r e s padj <0.05) res <- res[order(res padj<0.05)res<−res[order(respadj),]
resdata <- merge(as.data.frame(res), as.data.frame(counts(dds, normalized=TRUE)),by=“row.names”,sort=FALSE)
write.csv(resdata,file = “LC_1_vs_LC_2.csv”)
EdgeR的使用方法：
跟DESeq2一样，EdgeR输入矩阵数据，行名为sample，列名为gene；DESeq2不支持无生物学重复的数据，因此我选择了2个样本，3个生物学重复的数据。

exprSet_all <- read.table(file = “readcount”, sep = “\t”, header = TRUE, row.names = 1, stringsAsFactors = FALSE)
exprSet <- exprSet_all[,1:6]
group_list <- factor(c(rep(“LC_1”,3), rep(“LC_2”,3)))
设置分组信息，去除低表达量的gene以及做TMM标准化

exprSet <- exprSet[rowSums(cpm(exprSet) > 1) >= 2,]
exprSet <- DGEList(counts = exprSet, group = group_list)
exprSet <- calcNormFactors(exprSet)
使用qCML（quantile-adjusted conditional maximum likelihood）估计离散度（只针对单因素实验设计）

exprSet <- estimateCommonDisp(exprSet)
exprSet <- estimateTagwiseDisp(exprSet)
寻找差异gene(这里的exactTest函数还是基于qCML并且只针对单因素实验设计)，然后按照阈值进行筛选即可

et <- exactTest(exprSet)
tTag <- topTags(et, n=nrow(exprSet))
tTag <- as.data.frame(tTag)
write.csv(tTag,file = “LC_1_vs_LC_2_edgeR.csv”)
Summary
以上我主要针对单因素两两比较组进行差异分析，其实DESeq2和EdgeR两个R包都可以对多因素进行差异分析。

DESeq2修改以上代码的分组信息design参数以及在差异分析results函数中添加所选定的分组因素，其他代码基本一样，具体参照DESeq2手册

EdgeR则需要用Cox-Reid profile-adjusted likelihood (CR)方法来估算离散度，y <- estimateDisp(y, design)或者分别使用三个函数（y <- estimateGLMCommonDisp(y, design)，y <- estimateGLMTrendedDisp(y, design)， ）y <- estimateGLMTagwiseDisp(y, design)；然后差异表达分析也跟单因素分析不同，主要使用generalized linear model (GLM) likelihood ratio test 或者 quasi-likelihood(QL) F-test，具体代码可以参照EdgeR手册。

简单使用limma做差异分析
Posted: 五月 12, 2017 Under: Transcriptomics By Kai no Comments

首先需要说明的是，limma是一个非常全面的用于分析芯片以及RNA-Seq的差异分析，按照其文章所说：

limma is an R/Bioconductor software package that provides an integrated solution for analysing data from gene expression experiments.

在这我只是对其中的一种情况进行简单的总结，比如这个包可以处理RNA-Seq数据，我简单的以两个比较组进行分组为例，至于其他分组情况，请看limma说明文档，有非常详细的说明，非常亲民。

首先我们还是输入count矩阵，这里也跟其他差异分析R包一样，不要输入已经标准化的数据。顺便也加载下edgeR这个R包

library(limma)
library(edgeR)
counts <- read.table(file = “conut_all.txt”, sep = “\t”, header = TRUE, row.names = 1, stringsAsFactors = FALSE)
接着按照文档的说明以及limma包的习惯，我们需要对count进行标准化以及转化为log2的值，这里标准化的方法为TMM，使用edgeR里面的calcNormFactors函数即可

dge <- DGEList(counts = counts)
dge <- calcNormFactors(dge)
logCPM <- cpm(dge, log=TRUE, prior.count=3)
这里prior.count值我粗略理解为是为了防止log2()遇到过于小的值

然后跟其他包一样，设置分组信息

group_list <- factor(c(rep(“control”,2), rep(“siSUZ12”,2)))
design <- model.matrix(~group_list)
colnames(design) <- levels(group_list)
rownames(design) <- colnames(counts)
接下来就是常规的差异分析

fit <- lmFit(logCPM, design)
fit <- eBayes(fit, trend=TRUE)
output <- topTable(fit, coef=2,n=Inf)
sum(output$adj.P.Val<0.05)
到这里为止，我们主要是用了limma包里对RNA-Seq差异分析的limma-trend方法，该方法主要适用于样本间测序深度基本保持一致的情况，也就是说两个样本的文库(reads数目)大小相差的不悬殊（说明文档中是默认3倍以内？）

当文库大小在样本间变化幅度相当大的话，可以使用limma的voom方法，可按照下面的代码流程：

count数据的输入以及数据标准化还是跟之前的一样

counts <- read.table(file = “conut_all.txt”, sep = “\t”, header = TRUE, row.names = 1, stringsAsFactors = FALSE)
dge <- DGEList(counts = counts)
dge <- calcNormFactors(dge)
还是一样需要分组信息

group_list <- factor(c(rep(“control”,2), rep(“siSUZ12”,2)))
design <- model.matrix(~group_list)
colnames(design) <- levels(group_list)
rownames(design) <- colnames(counts)
接下来进行voom转化

v <- voom(dge, design, plot=TRUE)
其实可以不进行TMM标准化，直接用count数据进行voom转化，如：

v <- voom(counts, design, plot=TRUE)
最后就是普通的差异分析过程

fit <- lmFit(v, design)
fit <- eBayes(fit)
output <- topTable(fit, coef=2,n=Inf)
sum(output$adj.P.Val<0.05)
Summary
Limma长久以来就是一个非常流行的差异分析R包，其内容涉及的非常广泛，用于RNA-Seq只是其内容的一小部分，并且使其处理RNA-Seq数据也使用芯片类似线性模型下，并且按照其说法，limma包比其他基于负二项式分布模型的差异分析R包更加的优秀。

其实差异分析不外乎数据的标准化以及统计模型分析差异两个方面的作用，每个差异分析R包都有其自身的优点，个人理解，取舍在于自己的理解以及想法即可。

其实，自己对于limma包的理解还是比较粗浅的

转录组差异表达分析小实战（一）
Posted: 七月 28, 2017 Under: Transcriptomics By Kai no Comments

读文献获取数据
文献名称：AKAP95 regulates splicing through scaffolding
RNAs and RNA processing factors

查找数据：Data availability
The RIP-seq an RNA-seq data have been deposited in the Gene
Expression Omnibus database, with accession code GSE81916. All other data is
available from the author upon reasonable request.

获得GSE号：GSE81916

下载测序数据
https://www.ncbi.nlm./geo/query/acc.cgi?acc=GSE81916获取数据信息，并点击网址下方的ftp，下载测序数据

从https://trace.ncbi.nlm./Traces/study/?acc=PRJNA323422可知我们需要的mRNA测序编号为SRR3589956到SRR3589962

通过Apera下载SRR数据，这里以SRR3589956为例：

ascp -T -i /home/anlan/.aspera/connect/etc/asperaweb_id_dsa.openssh [e-mail].:sra/sra-instant/reads/ByRun/sra/SRR/SRR358/SRR3589956/SRR3589956.sra ./
转化fastq测序数据
通过sratoolkit工具将SRR文件转化为fastq格式的测序数据（写了个shell循环）

for i in ( s e q 5662 ) ; d o n o h u p f a s t q − d u m p − − s p l i t − 3 S R R 35899 (seq 56 62);do nohup fastq-dump --split-3 SRR35899 (seq5662);donohupfastq−dump−−split−3SRR35899{i} &;done
通过fastqc对每个fastq文件进行质检，用multiqc查看整体质检报告（对当前目录下的fastq测序结果进行质检，生成每个fq文件的质检报告总multiqc整合后统计查看）

fastqc *.fastq
multiqc ./
点击这个url可以查看我这个multiqc报告：http://www./data/multiqc_report.html

如果有接头或者质量值不达标的需要进行过滤，这次的数据质量都不错，因此直接进行比对即可
序列比对
安装hisat2软件，下载人类的hiast2索引文件

hisat2下载并安装：

ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/downloads/hisat2-2.1.0-Linux_x86_64.zip
unzip hisat2-2.1.0-Linux_x86_64.zip
下载hisat2的human索引

ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gz
tar zxvf hg19.tar.gz
用hisat2进行比对，测序数据放在data目录下，索引文件放在reference/index/hisat2/hg19目录下，SRR3589956-SRR3589958为人的测序数据

for i in KaTeX parse error: Undefined control sequence: \ at position 28: …do hisat2 -p 4 \̲ ̲-x ~/reference/…{i}_1.fastq -2 ./data/SRR35899KaTeX parse error: Undefined control sequence: \ at position 13: {i}_2.fastq \̲ ̲-S SRR35899i.sam >SRR35899${i}.log;done
用samtools将sam文件转化为bam文件，并使用默认排序

for i in ( s e q 5658 ) ; d o s a m t o o l s s o r t − @ 5 − o S R R 35899 (seq 56 58);do samtools sort -@ 5 -o SRR35899 (seq5658);dosamtoolssort−@5−oSRR35899{i}.bam SRR35899${i}.sam;done
reads计数
用htseq对比对产生的bam进行count计数

htseq安装，使用miniconda，省事！唯一的问题是htseq版本不是最新的，是0.7.2。想要最新版还是要正常安装，可参考http://www./thread-1847-1-2.html

conda install -c bioconda htseq
用htseq将对比后的结果进行计数

for i in KaTeX parse error: Undefined control sequence: \ at position 48: …m -r pos -s no \̲ ̲SRR35899{i}.bam ~/reference/genome/hg19/gencode.v26lift37.annotation.gtf
1>SRR35899 i . c o u n t 2 > S R R 35899 {i}.count 2>SRR35899 i.count2>SRR35899{i}_htseq.log;done
将3个count文件（SRR3589956.count，SRR3589957.count，SRR3589958.count）合并成一个count矩阵，这是就需要脚本来解决这个问题，不然其他方法会稍微麻烦点

#!/usr/bin/perl -w
use strict;

my $path = shift @ARGV;
opendir DIR, $path or die;
my @dir = readdir DIR;

my $header;
my @sample;
my %hash;
foreach my KaTeX parse error: Expected '}', got 'EOF' at end of input: …dir) { if (file =~ /^\w+.*.count/) {
push @sample, $file;
KaTeX parse error: Undefined control sequence: \t at position 12: header .= "\̲t̲file";
open my $fh, f i l e o r d i e ; w h i l e ( < file or die; while (< fileordie;while(<fh>) {
chomp;
next if ($_ =~ /^\W+/);
my @array = split /\t/, $_;
KaTeX parse error: Expected '}', got 'EOF' at end of input: hash{array[0]} -> {$file} = $array[1];
}
close KaTeX parse error: Expected 'EOF', got '}' at position 9: fh; }̲ } print "header\n";
map{
my $gene = ; p r i n t " _; print " ;print"gene";
foreach my KaTeX parse error: Undefined control sequence: \t at position 33: … print "\̲t̲".hash{KaTeX parse error: Expected 'EOF', got '}' at position 5: gene}̲ -> {file};
}
print “\n”;
}keys %hash;
按照接下来的剧本，应该讲count_matrix文件导入DESeq进行差异表达分析。但是从这篇文章的Bioinformatic analyses部分可以发现，作者的control组的2组数据是来自2个不同的批次（一个是SRR3589956，另外一个来源GSM1095127 in GSE44976），treat组倒是同一个批次（SRR3589957和SRR3589958）。但是对于Mouse cells来说，倒是满足2个control和2个treat都正常来自同个批次，因此打算重新用SRR3589959-SRR3589962重新做个一个count_matrix进行后续差异分析

转录组差异表达分析小实战（二）
Posted: 八月 14, 2017 Under: Transcriptomics By Kai no Comments

差异基因表达分析
我按照前面的流程转录组差异表达分析小实战（一），将小鼠的4个样本又重新跑了一遍，从而获得一个新的count文件：mouse_all_count.txt，有需要的话，可以下载下来进行后续的差异分析。

一般来说，由于普遍认为高通量的read count符合泊松分布，所以一些差异分析的R包都是基于负二项式分布模型的，比如DESeq、DESeq2和edgeR等，所以这3个R包从整体上来说是类似的（但各自标准化算法是不一样的）。

当然还有一个常用的R包则是Limma包，其中的limma-trend和limma-voom都能用来处理RNA-Seq数据（但对应适用的情况不一样）

下面准备适用DESeq2和edgeR两个R包分别对小鼠的count数据进行差异表达分析，然后取两者的结果的交集的基因作为差异表达基因。

DEseq2

library(DESeq2)
##数据预处理
database <- read.table(file = “mouse_all_count.txt”, sep = “\t”, header = T, row.names = 1)
database <- round(as.matrix(database))

##设置分组信息并构建dds对象
condition <- factor(c(rep(“control”,2),rep(“Akap95”,2)), levels = c(“control”, “Akap95”))
coldata <- data.frame(row.names = colnames(database), condition)
dds <- DESeqDataSetFromMatrix(countData=database, colData=coldata, design=~condition)
dds <- dds[ rowSums(counts(dds)) > 1, ]

##使用DESeq函数估计离散度，然后差异分析获得res对象
dds <- DESeq(dds)
res <- results(dds)

#最后设定阈值，筛选差异基因，导出数据(全部数据。包括标准化后的count数)
res <- res[order(res$padj),]
diff_gene <- subset(res, padj < 0.05 & (log2FoldChange > 1 | log2FoldChange < -1))
diff_gene_DESeq2 <- row.names(diff_gene)
resdata <- merge(as.data.frame(res), as.data.frame(counts(dds, normalized=TRUE)),by=“row.names”,sort=FALSE)
write.csv(resdata,file = “control_vs_Akap95.csv”,row.names = F)
最终获得572个差异基因（筛选标准为padj < 0.05, |log2FoldChange| > 1）

edgeR

library(edgeR)
##跟DESeq2一样，导入数据，预处理（用了cpm函数）
exprSet<- read.table(file = “mouse_all_count.txt”, sep = “\t”, header = TRUE, row.names = 1, stringsAsFactors = FALSE)
group_list <- factor(c(rep(“control”,2),rep(“Akap95”,2)))
exprSet <- exprSet[rowSums(cpm(exprSet) > 1) >= 2,]

##设置分组信息，并做TMM标准化
exprSet <- DGEList(counts = exprSet, group = group_list)
exprSet <- calcNormFactors(exprSet)

##使用qCML（quantile-adjusted conditional maximum likelihood）估计离散度（只针对单因素实验设计）
exprSet <- estimateCommonDisp(exprSet)
exprSet <- estimateTagwiseDisp(exprSet)

##寻找差异gene(这里的exactTest函数还是基于qCML并且只针对单因素实验设计)，然后按照阈值进行筛选即可
et <- exactTest(exprSet)
tTag <- topTags(et, n=nrow(exprSet))
diff_gene_edgeR <- subset(tTagKaTeX parse error: Expected 'EOF', got '&' at position 19: …le, FDR < 0.05 &̲ (logFC > 1 | l…table,file = “control_vs_Akap95_edgeR.csv”)
最终获得688个差异基因（筛选标准为FDR < 0.05, |log2FC| > 1）

取DESeq2和edgeR两者结果的交集

diff_gene <- diff_gene_DESeq2[diff_gene_DESeq2 %in% diff_gene_edgeR]
最终的差异基因数目为545个

head(diff_gene)
[1] “ENSMUSG00000003309.14” “ENSMUSG00000046323.8” “ENSMUSG00000001123.15”
[4] “ENSMUSG00000023906.2” “ENSMUSG00000044279.15” “ENSMUSG00000018569.12”
其他两个R包（DESeq和limma）就不在这尝试了，我之前做过对于这4个R包的简单使用笔记，可以参考下：
简单使用DESeq做差异分析
简单使用DESeq2/EdgeR做差异分析
简单使用limma做差异分析

GO&&KEGG富集分析
以前一直没有机会用Y叔写的clusterProfiler包，这次正好看说明用一下。

GO富集，加载clusterProfiler包和org.Mm.eg.db包（小鼠嘛），然后将ENSEMBL ID后面的版本号去掉，不然后面不识别这个ID，然后按照clusterProfiler包的教程说明使用函数即可。

library(clusterProfiler)
library(org.Mm.eg.db)

##去除ID的版本号
diff_gene_ENSEMBL <- unlist(lapply(diff_gene, function(x){strsplit(x, “\.”)[[1]][1]}))
##GOid mapping + GO富集
ego <- enrichGO(gene = diff_gene_ENSEMBL, OrgDb = org.Mm.eg.db,
keytype = “ENSEMBL”, ont = “BP”, pAdjustMethod = “BH”,
pvalueCutoff = 0.01, qvalueCutoff = 0.05)
##查看富集结果数据
enrich_go <- as.data.frame(ego)
##作图
barplot(ego, showCategory=10)
dotplot(ego)
enrichMap(ego)
plotGOgraph(ego)
KEGG富集，首先需要将ENSEMBL ID转化为ENTREZ ID，然后使用ENTREZ ID作为kegg id，从而通过enrichKEGG函数从online KEGG上抓取信息，并做富集

library(clusterProfiler)
library(org.Mm.eg.db)

##ID转化
ids <- bitr(diff_gene_ENSEMBL, fromType = “ENSEMBL”, toType = “ENTREZID”, OrgDb = “org.Mm.eg.db”)
kk <- enrichKEGG(gene = ids[,2], organism = “mmu”, keyType = “kegg”,
pvalueCutoff = 0.05, pAdjustMethod = “BH”, qvalueCutoff = 0.1)
##查看富集结果数据
enrich_kegg <- as.data.frame(kk)
##作图
dotplot(kk)
到这里为止，转录组的差异表达分析算是做完了，简单的来说，这个过程就是将reads mapping 到reference上，然后计数获得count数，然后做差异分析，最后来个GO KEGG，over了。。。

对于mapping和计数这两部还有其实还有好多软件，具体可见文献：Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis，有时间可以都尝试下。

至于GO && KEGG这两步，对于人、小鼠等模式物种来说，不考虑方便因素来说，完全可以自己写脚本来完成，数据可以从gene ontology官网下载，然后就是GO id与gene id相互转化了。KEGG 也是一样，也可以用脚本去KEGG网站上自行抓取，先知道URL规则，然后爬数据即可。