RNA-seq数据基因水平表达差异分析 129

这篇博文来自于Bioc2014的学习笔记。

使用Bioconductor分析RNA-seq数据,获取基因水平表达差异结果主要分三步:记数(counts) -> 质量控制(diagnostics) -> 统计分析(statistics)

在Bioconductor中,用于RNA-seq数据分析的包很多,这里介绍三种edgeR, DESeq2以及voom。

edgeR

为了避免意外,我们开始一个全新的r session。

rm(list=ls())

载入edgeR。

library(edgeR)
## Loading required package: limma

载入数据

load(url("http://qiubio.com/bioconductor/RNA-seq/ds1.Rdata"))
ls()
## [1] "counts"
class(counts)
## [1] "matrix"
dim(counts)
## [1] 30727    27
head(counts[, 1:7], 3)
##                 ES.07985 DE.07981 GT.66339 FG.08004 PE.07980 FE.66350
## ENSG00000254468        0        0        0        1        0        0
## ENSG00000177951       44       50       24       37       38       41
## ENSG00000188076        0        0        0        0        0        0
##                 ES.66342
## ENSG00000254468        1
## ENSG00000177951      162
## ENSG00000188076        0

为载入的数据分组。这里使用列名提供的信息分组。

grp <- as.factor(substr(colnames(counts), 1, 2)) ##substr截取列名的前两个字母
table(grp)
## grp
## DE ES FE FG GT IS PE PH 
##  3  3  4  3  3  2  6  3
o <- order(grp)
pairs(log2(1+counts[,o[1:7]]), pch=".",lower.panel=NULL) ##这个图给我们样品间关系的一揽印象

1

建立一个DGEList对象,并对数据进行标准化(normalize)

d <- DGEList(counts=counts, group=grp)
d <- calcNormFactors(d)
d$samples
##          group lib.size norm.factors
## ES.07985    ES   769440       0.9721
## DE.07981    DE   947798       1.0560
## GT.66339    GT   548227       0.8668
## FG.08004    FG   726088       1.2161
## PE.07980    PE   688903       1.0602
## FE.66350    FE   521573       1.1357
## ES.66342    ES  2281754       0.9036
## DE.66333    DE  2834208       0.9781
## GT.66341    GT  1596550       0.7663
## FG.66346    FG  2070020       1.1650
## PE.66344    PE  1978504       1.0385
## FE.66331    FE  1517276       1.0618
## PH.66332    PH   840380       1.3760
## PH.66336    PH   852548       1.3735
## PE.66345    PE   915042       1.2394
## PE.66330    PE   938173       1.2296
## FE.66348    FE   629687       1.1131
## FE.66334    FE   647313       1.1129
## ES.66335    ES  2048439       0.7027
## DE.66349    DE  1843746       0.7489
## GT.66337    GT  1761089       0.8636
## FG.66351    FG  1417106       0.9454
## PE.66338    PE  2317354       0.8958
## PH.66340    PH  4087566       1.1389
## PE.66329    PE  3700620       1.0793
## IS.66347    IS 12583285       0.5621
## IS.66343    IS 33286768       0.9297

这里的norm.factors计算方法有3种,分别是TMM, RLE和upperquartile。默认使用TMM方法。

接下来要去除低表达值的数据,这一步可以减少很多计算量,不做也不影响最终的结果。

dim(d)
## [1] 30727    27
cps <- cpm(d) ## count per million
k <- rowSums(cps>=1) > 2
d <- d[k,]
dim(d)
## [1] 21707    27

我们看到数据几乎减少了1/3。

然后进行聚类分析。这一步是质量控制的关键一步。如果同一组的biological replicates不 能相聚在一起,必要时需要考虑重新实验。如果条件允许,这一步就要抛弃一些意外的数据。

cols <- as.numeric(d$samples$group)
plotMDS(d,col=cols)

2

绘制multidimensional scaling plot有很多种手段。LFC就是log2 fold change, BCV就是 biological coefficient of variation. 为了快速绘图,默认情况下只是抽取了top 500个值。

par(mfrow=c(2,2))
plotMDS(d, col=cols, main="500 / lLFC")
plotMDS(d, col=cols, method="bcv", main="500 / BCV")
plotMDS(d, col=cols, top=2000, main="2000 / lLFC")
plotMDS(d, col=cols, top=2000, method="bcv", main="2000 / BCV")

3

而后构建模型。

mm <- model.matrix(~-1+grp)
mm
##    grpDE grpES grpFE grpFG grpGT grpIS grpPE grpPH
## 1      0     1     0     0     0     0     0     0
## 2      1     0     0     0     0     0     0     0
## 3      0     0     0     0     1     0     0     0
## 4      0     0     0     1     0     0     0     0
## 5      0     0     0     0     0     0     1     0
## 6      0     0     1     0     0     0     0     0
## 7      0     1     0     0     0     0     0     0
## 8      1     0     0     0     0     0     0     0
## 9      0     0     0     0     1     0     0     0
## 10     0     0     0     1     0     0     0     0
## 11     0     0     0     0     0     0     1     0
## 12     0     0     1     0     0     0     0     0
## 13     0     0     0     0     0     0     0     1
## 14     0     0     0     0     0     0     0     1
## 15     0     0     0     0     0     0     1     0
## 16     0     0     0     0     0     0     1     0
## 17     0     0     1     0     0     0     0     0
## 18     0     0     1     0     0     0     0     0
## 19     0     1     0     0     0     0     0     0
## 20     1     0     0     0     0     0     0     0
## 21     0     0     0     0     1     0     0     0
## 22     0     0     0     1     0     0     0     0
## 23     0     0     0     0     0     0     1     0
## 24     0     0     0     0     0     0     0     1
## 25     0     0     0     0     0     0     1     0
## 26     0     0     0     0     0     1     0     0
## 27     0     0     0     0     0     1     0     0
## attr(,"assign")
## [1] 1 1 1 1 1 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$grp
## [1] "contr.treatment"

依据模型估算离散度(dispersion)。你一定会问什么是离散度啊?在McCarthy, DJ, Chen, Y, Smyth, GK (2012). Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation. Nucleic Acids Research 40, 4288-4297.一文中定义离散度为生物变异系数(biological coefficient of variation)的平方。这里需要解释一下什么是技术重复(technical replicates),什么是生物重复(biological replicates)。技术重复就是同样的一个生物样品,使用相同的操作技术,两次或者两次以上实验。而生物重复是指不同的生物样品,使用相同的操作技术,两次或者两次以上的实验。而 biological variation在上文中专指排除technical variation以外的因为不同生物样品导致的差异。 \(Total CV^2 = Technical CV^2 + Biological CV^2\)

d <- estimateGLMCommonDisp(d,mm)
d <- estimateGLMTrendedDisp(d,mm)
## Loading required package: splines
d <- estimateGLMTagwiseDisp(d,mm)

其实离散度也可以是质量控制的一部分。我们使用plotBCV来察看同一组内不同样品的不同表达水平值的方差分析。

plotBCV(d)

4

从上图中,可以看到红线所在的Common dispersion有0.34左右。但是当我们调用d$common.dispersion时,它却只有0.12。这是怎么回事呢?这是因为求方根之后才是真正的值。

d$common.dispersion
## [1] 0.1182
sqrt(d$common.dispersion)
## [1] 0.3438

##请无视此段。对于BCV的估计,正常情况下,必须要有很多样品才能比较准确,而生物实验是不可能提供成百上千的重复的,所以通常的手段是将表达值相近的基因等同看待。但是这样还是会对BCV过度估计。对于过度估计,在GLM中可以penalized maximum likelihood estimation来进行修正。

对于plotBCV给出的红线,它在0.2~0.4之间都是容易接受的。但如果超过0.5,则需要考虑实验样品之间的差异是否过大了。

plotMeanVar(d, show.raw=TRUE, show.tagwise=TRUE, show.binned=TRUE)

5

在上图中,很明显地显示出基因水平的差异(variance)随着表达值的升高而升高。这一点会让有些人觉得很困惑,因为很多人会认为当表达值升高时,变异系数应该下降才对。但其实随着表达值升高,变异系数下降只符合技术重复带来的差异变化,而生物重复带来的差异变化却并不是这样。从dispersion的图中我们可以看到,它是一个两头翘的分布。因为生物重复差异在整体的变异系数中占主导地位,所以我们会观察到其基因水平的变异系数反而会随着表达值的升高而升高。

对于基因芯片,我们熟悉了MAplot,对于RNAseq同样可以绘制MAplot。MAplot的y轴的中值应该在0线上。

par(mfrow=c(1,2))
plotSmear(d, pair=c("ES","DE"), ylim=c(-5,5))
plotSmear(d, pair=c("DE","GT"), ylim=c(-5,5))

6

接下来的步骤就是统计分析的步骤了。先做广义线性拟合(Generalized linear models),而后构建比较结构(construct contrast),最后依据比较结构进行比较,给出差异表达的基因。

f <- glmFit(d,mm)
con <- makeContrasts("DE-ES"=grpDE-grpES,levels=colnames(mm))
con
##        Contrasts
## Levels  DE-ES
##   grpDE     1
##   grpES    -1
##   grpFE     0
##   grpFG     0
##   grpGT     0
##   grpIS     0
##   grpPE     0
##   grpPH     0
lrt <- glmLRT(f,contrast=con)
topTags(lrt,20)
## Coefficient:  1*grpDE -1*grpES 
##                  logFC logCPM    LR    PValue       FDR
## ENSG00000095596  8.532  6.658 295.5 3.189e-66 6.923e-62
## ENSG00000076356  7.978  7.297 280.2 6.962e-63 7.556e-59
## ENSG00000164292 10.583  6.440 262.0 6.295e-59 4.555e-55
## ENSG00000185823  7.349  5.685 241.4 1.936e-54 1.051e-50
## ENSG00000125798  9.891  8.044 227.5 2.074e-51 9.003e-48
## ENSG00000132130 11.314  6.398 223.1 1.880e-50 6.801e-47
## ENSG00000126005 -4.501  6.464 218.4 2.059e-49 6.384e-46
## ENSG00000104332  4.707  6.784 217.4 3.278e-49 8.895e-46
## ENSG00000158815 14.932  7.295 202.7 5.407e-46 1.304e-42
## ENSG00000119242  6.266  6.655 198.9 3.606e-45 7.828e-42
## ENSG00000182798  5.062  7.548 185.1 3.694e-42 7.289e-39
## ENSG00000266283 11.071  8.362 181.4 2.338e-41 4.230e-38
## ENSG00000120875  7.807  5.529 178.2 1.193e-40 1.991e-37
## ENSG00000110799  7.639  4.880 177.5 1.693e-40 2.624e-37
## ENSG00000121966  7.529  5.850 174.7 6.890e-40 9.971e-37
## ENSG00000185155  6.432  4.827 171.7 3.206e-39 4.349e-36
## ENSG00000164736 12.848  5.288 171.5 3.560e-39 4.546e-36
## ENSG00000124942  6.854  7.487 170.6 5.400e-39 6.512e-36
## ENSG00000141448 13.186  6.579 166.0 5.478e-38 6.259e-35
## ENSG00000110148  6.790  6.114 162.3 3.517e-37 3.817e-34

最后,得到差异表达基因列表后,我们找出几个来做一个检查,看看结果对不对。

cps <- cpm(d)
o <- order(colnames(counts))
barplot( cps["ENSG00000095596",o], col=cols[o], las=2)

7

输出结果至文件。

tt <- topTags(lrt, n=Inf)$table
write.table(tt, file="LRT1.xls", row.names=FALSE, sep="\t", quote=FALSE)

DESeq2 载入DESeq2

library(DESeq2)

我们载入与上一例相同的数据。

load(url("/bioconductor/RNA-seq/ds1.Rdata"))
ls()
##  [1] "cols"   "con"    "counts" "cps"    "d"      "f"      "grp"   
##  [8] "k"      "lrt"    "mm"     "o"      "tt"
head(counts[, 1:7], 3)
##                 ES.07985 DE.07981 GT.66339 FG.08004 PE.07980 FE.66350
## ENSG00000254468        0        0        0        1        0        0
## ENSG00000177951       44       50       24       37       38       41
## ENSG00000188076        0        0        0        0        0        0
##                 ES.66342
## ENSG00000254468        1
## ENSG00000177951      162
## ENSG00000188076        0
grp <- as.factor(substr(colnames(counts), 1, 2)) ##substr截取列名的前两个字母
dds <- DESeqDataSetFromMatrix(counts, colData=data.frame(grp), design=formula(~-1+grp)) ##这里的colData必须是一个DataFrame或者data.frame。每一行都对应着counts中的一列。design中的公式和limma中的方法一致。
design(dds)
## ~-1 + grp
dds <- DESeq(dds) ##一步到位,不需要太多步骤
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
res <- results(dds) ##得到结果
res
## log2 fold change (MAP): grp PH vs DE 
## Wald test p-value: grp PH vs DE 
## DataFrame with 30727 rows and 6 columns
##                  baseMean log2FoldChange     lfcSE       stat    pvalue
##                 <numeric>      <numeric> <numeric>  <numeric> <numeric>
## ENSG00000254468    0.4287       -1.05357    1.6490  -0.638900  0.522888
## ENSG00000177951  101.9852        0.30692    0.2534   1.211289  0.225785
## ENSG00000188076    0.1621       -0.01532    1.7535  -0.008735  0.993031
## ENSG00000177947    2.6249       -2.84301    1.2344  -2.303244  0.021265
## ENSG00000254559  496.3435       -0.93363    0.3140  -2.973341  0.002946
## ...                   ...            ...       ...        ...       ...
## ENSG00000248275  19.84377      1.595e+00    0.4163  3.830e+00 0.0001281
## ENSG00000248103   0.01374      2.614e-05    0.8715  3.000e-05 0.9999761
## ENSG00000238035 127.55955     -3.819e-01    0.2806 -1.361e+00 0.1736234
## ENSG00000250765   1.86059     -4.989e-01    1.0464 -4.768e-01 0.6335268
## ENSG00000230178   0.00751      5.954e-07    0.3614  1.648e-06 0.9999987
##                      padj
##                 <numeric>
## ENSG00000254468        NA
## ENSG00000177951   0.37016
## ENSG00000188076        NA
## ENSG00000177947   0.05871
## ENSG00000254559   0.01137
## ...                   ...
## ENSG00000248275 0.0007656
## ENSG00000248103        NA
## ENSG00000238035 0.3038560
## ENSG00000250765 0.7615653
## ENSG00000230178        NA
mcols(res, use.names=TRUE) ##结果中各例的意义
## DataFrame with 6 rows and 2 columns
##                        type                          description
##                 <character>                          <character>
## baseMean       intermediate          the base mean over all rows
## log2FoldChange      results log2 fold change (MAP): grp PH vs DE
## lfcSE               results         standard error: grp PH vs DE
## stat                results         Wald statistic: grp PH vs DE
## pvalue              results      Wald test p-value: grp PH vs DE
## padj                results                 BH adjusted p-values
plotMA(res)

8

plotDispEsts(dds)

9

write.csv(as.data.frame(res), file="results.csv") ##输出结果

voom

library(limma)
v <- voom(d, mm, plot=TRUE) ##这里的d就是edgeR中的d,下面就是标准的limma流程了。

10

vf <- lmFit(v,mm)            # 'mm' 前文定义了
cf <- contrasts.fit(vf,con)  # 'con' 前文定义了
cf <- eBayes(cf)
topTable(cf)                 # 
##                  logFC AveExpr      t   P.Value adj.P.Val     B
## ENSG00000104332  4.682  5.3343  16.82 4.643e-17 1.008e-12 28.60
## ENSG00000182798  5.049  2.3421  16.02 1.802e-16 1.781e-12 27.44
## ENSG00000126005 -4.485  5.7958 -15.47 4.669e-16 1.781e-12 26.29
## ENSG00000147869  7.395  1.6008  15.35 5.733e-16 1.781e-12 26.23
## ENSG00000076356  7.788  5.7036  15.35 5.745e-16 1.781e-12 25.22
## ENSG00000158815 11.806  1.6305  15.55 4.016e-16 1.781e-12 24.92
## ENSG00000164292  9.515  5.0669  15.50 4.418e-16 1.781e-12 24.89
## ENSG00000164736  9.722  0.6449  15.09 9.219e-16 2.501e-12 24.24
## ENSG00000003056 -3.454  8.4616 -13.94 7.795e-15 1.538e-11 23.84
## ENSG00000125798  8.925  7.0947  14.11 5.645e-15 1.361e-11 22.79

我们比较不同方法给出的头几行,差别还真是大啊。 问题来了,我们应该相信哪种方法呢?我只能说,都算一下,然后每种方法给出的结果里都挑几个跑一下qRT-PCR。有人问是不是取多种办法相同的的部分会验证率高些?我的经验,真不一定。如果data好,用哪种方法验证率都高,如果data不行,那只能看看哪种方法更合适了。

129 thoughts on “RNA-seq数据基因水平表达差异分析

  1. Reply xiaoma 8月 7,2014 4:18 上午

    “必要时需要考虑重新实验” 是多么不能接受:(

  2. Reply jingming 9月 22,2014 9:21 上午

    博主你好,
    现在我手上有Expression value,RPKM,Ratio of unique to total (transcript reads)这些RNA_SEQ的数据,如何分析啊?

  3. Reply biogg 11月 26,2014 2:08 上午

    大赞,学习了。

  4. Reply sun 12月 18,2014 1:38 上午

    请问手头上只有RPKM值和raw counts 数据。能用edgeR或者DEseq进行分析吗?

  5. Reply mkl 12月 23,2014 2:48 上午

    LZ你好。手头现在只有共7个组织的50多个sample的RPKM值,想用R分析一下这7个组织里差异表达的基因,请问该用什么方法?我现在有点晕的原因是,7个组织分别包括了不同数目的样本,那么该用什么统计方法呢?

  6. Reply funnybomb 12月 24,2014 2:42 上午

    请问构建模型
    mm <- model.matrix(~-1+grp)
    这里为什么要-1+grp

  7. Reply mary zhou 12月 24,2014 2:57 上午

    为什么我每次在做图表相关的语句的时候都会出现R code execution error这个错误呢

  8. Reply Huiying 12月 26,2014 12:50 上午

    您好!想请教一下,ENCODE的标准中提到RPKM的组内皮尔森相关系数应大于0.92,我用HTSeqcount计算出来的raw count计算了组内相关系数,结果基本都符合要求,相关性都很好,但是我用cuffdiff计算出来的FPKM值来计算组内相关系数,却有多个地方相关性很差(<0.5),我感到很奇怪,想问一下其中可能的原因,以及raw count,RPKM,FPKM三者之间差别会如此之大么?谢谢!

    • Reply admin 12月 26,2014 10:32 上午

      数据点多吗?如果足够多,可能还需要再做一个cov.test得一个pvalue比较好。关于选择相信哪一个,并不是100%重要。重要的是你是否有其它的证据相互支撑。

      • Reply Huiying 12月 28,2014 12:33 下午

        count table里大概有将近四万个基因,然后每个组有六个样本。然后弱弱的问一句cov.test是什么呀?听起来很耳熟但是死活想不起来,搜索了一下也没有得到靠谱的结果。另外,您所指的其他证据支撑是指其他生理指标么?还是另有所指?非常感谢您的指导!谢谢!

  9. Reply akun 1月 8,2015 3:30 下午

    刚进入RNA-seq这个圈,什么都还不太懂,刚测的数据返回来一大堆图表,好茫然,好恼火,望大神些多帮忙!

  10. Reply akun 1月 21,2015 2:33 下午

    你好,刚做了RNA-seq,回来了4个DEG的database,分别为材料A-处理 VS 材料B-处理,材料A-处理 VS 材料A-对照,材料B-处理 VS 材料B-对照,材料A-对照 VS 材料B-对照,现在需要找出同时存在于数据库材料A-处理 VS 材料B-处理,材料A-处理 VS 材料A-对照中的DEG,但这些DEG又不能出现在材料B-处理 VS 材料B-对照,材料A-对照 VS 材料B-对照的数据库中,这该怎么办呢?请前辈帮忙指点指点。

    • Reply admin 1月 22,2015 9:14 上午

      有两种方式,一种是使用contrast matrix, 另一种是分别自行处理后再使用set操作对它们做集合处理。后者较为简单,所以大多数人都直接采用后者。你可以搜索setdiff, interset, union等。

  11. Reply Ting 1月 21,2015 7:49 下午

    您好!好不容易翻过来。。。想请教您一下,对于一个基因的各个isoform来说,它们的relative abundance是什么呢?应该是0-1之间的一个box图,是否有准确的定义呢?

    • Reply admin 1月 22,2015 8:55 上午

      relative abundance我没有做过,不好意思。如果按照我的理解,应该是使用一个标准做为1,其余的与之比较即可。所以很有可能会大于1.它只是一个相对比例。当然我的理解可能不一定正确。

  12. Reply akun 1月 26,2015 4:48 下午

    你好,请问如何找出两个数据集A与B的交集,也就是A中上调(下调)的基因同时在B中也上调(下调),请问方便email交流吗?我想请你帮我看看我的EXCEL表,谢谢!

  13. Reply akun 1月 26,2015 5:00 下午

    还有请问如何将excel中的某一列或几列数据取出作为一个独立的数据集或csv文件?

  14. Reply yylin1984 2月 8,2015 8:17 上午

    您好,完全按照您的步骤运行DESeq,但当运行到
    > dds <- DESeq(dds)
    错误于DESeq(dds) :
    betaPrior=TRUE can only be used if the design has an intercept.
    if specifying + 0 in the design formula, use betaPrior=FALSE

    出现错误,请指正

  15. Reply LazyXuan 3月 13,2015 9:47 下午

    博主您好!我想请教一下您,利用要分析lncRNA的差异表达情况,用bowtie建索引并用tophat进行比对可以么?它们能处理非编码的序列么?急求答案,谢谢!谢谢!

    • Reply admin 3月 14,2015 8:23 下午

      对于lncRNA, 分析方法与普通的mRNA是类似的,区别一般在于注释的那一步。所以,你的问题的答案是可以。

      • Reply super 5月 25,2016 6:42 下午

        请问我如果已经把样本map到参考基因组上并且得到了BAM文件,我现在也下载了lncRNA的GTF文件,想分析lncRNA的差异表达。那么我是要把这个lncRNA的GTF注释文件和参考基因组的GTF注释文件合并成一个注释文件重新做tophat?或者我可以直接基于这个组合后的注释文件用cuffdiff或者htseq-count+edgeR做差分分析?还是我可以直接基于lncRNA的GTF注释文件直接用cuffdiff或者htseq-count+edgeR做差分分析?还是上述三种方法都不对? 谢谢!!!

      • Reply super 5月 26,2016 4:37 上午

        接上一个问题,问题2: 此外,除了已有的GTF文件,我需要把lncRNA的fasta文件也加到基因组里的genome.fa里去一起处理么?还是单独处理(不考虑human.fa和human.gtf,只基于ncRNA.fa和ncRNA.gtf做tophat和差分表达分析)就好?

      • Reply super 5月 26,2016 4:40 上午

        还有一个关于lncRNA的问题3,我在ensembl的主页上发现ensembl很多物种都有fasta文件即ncrna.fa,但是很多物种并没有GTF注释下载。。。我只好自己google,我找到一个网站http://www.noncode.org/,但好像也只有人和老鼠的GTF文件。 请问如果这样我如何处理呢? 我的方案是直接把我的fastq样本文件用bowtie对齐到lncRNA的fasta文件上,然后再用RSEM估算counts,最后再用edgeR分析,不知道我的方案可以么?关于lncRNA的差分表达方面连问了三个问题,劳烦老师了。

  16. Reply populus 7月 16,2015 1:28 下午

    博主,您好。在做差异表达分析之前,如何确定有效研究基因呢?我看有些文献说是所有样本的所有重复中都要检测到表达,将这部分基因保留做后续的差异表达分析。不知道这个一般用什么标准合适呢?k =1) > 2,这行的意思是否就是对于1行各样本的cps值中,有2个样本的表达值达到>=1的标准?谢谢

  17. Reply populus 7月 16,2015 1:32 下午

    rowSums(cps>=1) > 2,抱歉没显示全。这个2值一般如何设置较好?谢谢

    • Reply admin 7月 17,2015 8:03 下午

      我想应该是cpm吧。cps是什么?

      • Reply populus 7月 17,2015 9:35 下午

        您好。cps <- cpm(d)。这块我大概了解了。现在有几个问题,(1)我从JGI下载的gff3用在cufflink或者用Rsubread(或HTSeq-counts)做read counts。但是这个gff3中,对于同一个基因下可能有多个mRNA,每个mRNA下又有多个exons,其中存在相同位置的exon在不同转录本下重复列出的情况,这种会否对RPKM值的计算和read counts造成影响?应该如何处理?(2)在通过cpm值计算RPKM时,gene长度应该如何计算呢?是把基因内的所有exons长度相加(包括具有重叠区的exons)还是怎么处理呢?(3)在做read counts技术过程中,往往对gff3的支持不是很好,主要在于这些exons的parent通常都是转录本id,这样计数结果都是转录本的计数,而不是gene的计数,请问如何解决?问题有点多,还请见谅,非常感谢。下面是gff3的一个例子。

        Chr01 phytozomev10 gene 8391 12483 . – . ID=Po.001G000400.v3.0;Name=Po.001G000400
        Chr01 phytozomev10 mRNA 8391 12209 . – . ID=Po.001G000400.3.v3.0;Name=Po.001G000400.3;pacid=27046907;longest=0;Parent=Po.001G000400.v3.0
        Chr01 phytozomev10 exon 12036 12209 . – . ID=Po.001G000400.3.v3.0.exon.1;Parent=Po.001G000400.3.v3.0;pacid=27046907
        Chr01 phytozomev10 exon 11082 11177 . – . ID=Po.001G000400.3.v3.0.exon.2;Parent=Po.001G000400.3.v3.0;pacid=27046907
        Chr01 phytozomev10 exon 10900 10972 . – . ID=Po.001G000400.3.v3.0.exon.3;Parent=Po.001G000400.3.v3.0;pacid=27046907
        Chr01 phytozomev10 exon 8391 8666 . – . ID=Po.001G000400.3.v3.0.exon.4;Parent=Po.001G000400.3.v3.0;pacid=27046907
        Chr01 phytozomev10 mRNA 8391 11349 . – . ID=Po.001G000400.4.v3.0;Name=Po.001G000400.4;pacid=27046908;longest=0;Parent=Po.001G000400.v3.0
        Chr01 phytozomev10 exon 11253 11349 . – . ID=Po.001G000400.4.v3.0.exon.1;Parent=Po.001G000400.4.v3.0;pacid=27046908
        Chr01 phytozomev10 exon 11082 11177 . – . ID=Po.001G000400.4.v3.0.exon.2;Parent=Po.001G000400.4.v3.0;pacid=27046908
        Chr01 phytozomev10 exon 10900 10972 . – . ID=Po.001G000400.4.v3.0.exon.3;Parent=Po.001G000400.4.v3.0;pacid=27046908
        Chr01 phytozomev10 exon 8391 8666 . – . ID=Po.001G000400.4.v3.0.exon.4;Parent=Po.001G000400.4.v3.0;pacid=27046908

        • Reply admin 7月 21,2015 10:19 下午

          问题1, 不影响。其实你的问题3就回答了这个问题。对于这些计数软件,它都是找指定的标签,然后把指定标签下相同名称的exon全部合并起来进行计数,合并形成一个合集,统计落在这个区域内的reads数。
          问题2, 从cpm至rpkm的转换其实并不简单,因为我试过,得到的数据总是和cufflinks给的结果有小小的出入。我的办法就是cpm/coverageLength.这个coverageLength就是将exons的区域reduce之后合并起来,按中文的讲法就是所有的exons长度相加,重叠的部分只数一次。
          问题3,把id attribute改成Parent试试。

  18. Reply super 7月 17,2015 6:31 上午

    欧老师,请问edgeR和DESeq的论文里(http://www.nature.com/nprot/journal/v8/n9/full/nprot.2013.099.html)的simple design和complex design有何区别?你这篇blog用的是complex design,即加了一步mm <- model.matrix(~-1+grp)。

    • Reply admin 7月 17,2015 8:04 下午

      你可以参考limma的说明文档,那里比较明确。使用(~-1+grp)的意思就是默认线性拟合时没有截距。

      • Reply super 1月 18,2016 6:58 上午

        从介绍edgeR/DESeq的paper(http://www.nature.com/nprot/journal/v8/n9/full/nprot.2013.099.html)上看 complex design (GLM模型)主要是多组分析 (i.e., beyond two-group comparisons) 。
        (1)我的理解是如果是simple design,适用于两组分析且组内差异比较小,即[c1,c1,c1]vs[c2,c2,c2]。是么?
        此外
        (2). 如果只有两组样本,complex design适合做两组分析么?
        (3).如果是,是什么情况下做呢?是不是组内方差很大的情况?

      • Reply super 1月 18,2016 9:07 上午

        我在同一个样本集:三个control(C1) 三个case(C2), 3vs3 上试验了,结果发现GLM的(complex design)结果要比simple的多了很多DE genes。所以请专家详细解答。

        • Reply admin 1月 19,2016 10:42 上午

          其实我对edgeR的理解也不是很深。从我了解的来说,simple design在做dispersion估计的时候默认是使用qCML (quantile-adjusted conditional maximum likelihood)的。qCML受限于只能做单一因素的分析。所以在使用GLM的时候,做dispersion估计就换了一种方法,Cox-Reid profile-adjusted likelihood (CR) method。因为使用了不同的方法,所以最后差异分析可能会有不同。但理论上应该差异不大才是。
          无论是几组分析,组内差异如何,GLM模型都是合适的。有些人(edgeR的开发者之一Aaron Lun)甚至只喜欢使用GLM方法,但是他比较年经,对两者差异的理解不一定深刻。以我的经验而言,简单比较我就使用simple design。复杂分析才使用GLM。因为这是没有倾向性的办法。

  19. Reply jiazhongguo 7月 23,2015 7:36 上午

    想请教一个mapping中的问题:我拿到的RNA是双末端read,用tophat进行读段的有参定位。最近,用IGV来看定位结果,发现一对reads中,右端和左端read相聚到几千bp。但是用tophat定位时,-mate-inner-dist也就是-r参数我使用的默认值啊。然后,我与其他人交流,刚开始达成的共识是:1. tophat并不是把插入片段距离作为绝对的map标准,换句话说-r参数主要是参考作用,目的是提高比对准确性和比对速度;2 就是文库构建过程中的随机打断后,待测序片段并不一定都是按要求的大小,可能有长达几KB的,所以中间为测片段很长。。。。这是我最初的想法,可转念一想又觉得不对。。。。
    所以想问你是如何理解的!或者我这个认识你觉得对不

    • Reply admin 7月 24,2015 2:52 下午

      那个r值其实只是一个均值,还有一个值是它的标准差,这个左端reads和右端reads,理论上它们的平均值应该是你提供的r值,或者附近。对于某一具体的例子,可能得考虑具体的mapping的情况。一般你所描述的这种情况,会发生在它确实已经找不到更好的配对的位置了,所以只好把这种相距几千bp的记录下来。这应该算是bug的一种吧。看具体你怎么理解,有的人认为可以接受,有的人认为不能接受。

      • Reply jiazhongguo 7月 26,2015 8:49 下午

        记忆中,你在博客上回复过别人说的是基于tophat的参数设置,是不可能出现一个读段对(pair)的两个read分别定位到两条染色体上的情况的吧!如果这么说,map到同一染色体但是距离远超过设定值的只能算是bug。

  20. Reply jiazhongguo 7月 27,2015 9:19 上午

    不好意思,追加一个问题或疑问。
    1. 参考基因组应该是单倍体基因组吧,而且也只给出了一条链的信息吧!另外,根据参考基因组最基本的组装原理,组装是依赖OVERLAP的片段,一点一点连接起来。那这样子的话,如果不考虑gap,这些片段岂不是位于一条同源染色体上吗。可是如果是这样子的话,双末端数据定位结果为什么会有“+ -”“-+”“——”“++”这四种类型呢。 有点没搞懂!可否解释下!

    2. 你是否做过CNV检测呢。以及RAD测序也就是简化基因组测序!可有比较好的资源推荐!

    • Reply admin 7月 28,2015 2:09 下午

      问题1,如果你的测序不是strand specific的话,那就可能有四种可能啊。否则,应该是以两种可能为主。
      问题2,我没有做过CNV分析。不好意思,帮不上你。

  21. Reply ping 8月 5,2015 1:54 下午

    博主好,
    请问在edgeR分析plotSmear(d, pair=c(“ES”,”DE”), ylim=c(-5,5))中,图中靠近Y轴的黄色点与其他黑色点有什么区别呢?为什么黄色点的分布几乎是靠近Y轴的直线,其他黑点的分布呈鱼尾状?
    DESeq的plotMA(res)图中,红色点与黑色点有什么区别呢?
    topTags(lrt,10)查看的是10个表达差异最大的基因吗?logFC列绝对值最大是表示表达差异最大吗?
    谢谢博主!

    • Reply admin 8月 6,2015 8:44 上午

      您的这个问题还真是有趣。我以前都没想过。经过我初步的研究,可以看到它们都分布在low CPM区域。它的计算是通过两组当中任意一组中的最小值所在的基因以及logFC在5倍quantile以外的值,这些值基本上可以视为outline. 黑点的形状就是真实数据的形状,这是可以想象的,因为表达值越低的基因,其相对的变化值就会越大。

      plotMA的红色的点应该就是由alpha参数控制的,它对数据点通过adjusted p-value做了一个过滤。

      topTags是取出排序最靠前的10条数据,默认是由PValue来排的,所以logFC的绝对值不一定是最大的。

  22. Reply ping 8月 6,2015 12:46 下午

    非常感谢博主的回答。
    根据您的回答,plotSmear 图中,我是不是可以理解为可以视黄色的点为不用管的噪音?

    我还是不太明白plotMA红色点与黑色点的区别。您说红色的点由alpha参数控制?哪里来的alpha参数?
    您说对数据点通过adjusted p-value做了一个过滤。是指红色点符合设定的p-value 的值而黑色的不是吗?您这里说的p-value是FDR吧

    • Reply admin 8月 10,2015 8:05 上午

      你可以读一读plotMA的说明啊。plotMA(object, alpha = 0.1, main = “”, ylim, …),
      alpha
      the significance level for thresholding adjusted p-values

  23. Reply ping 8月 6,2015 1:45 下午

    博主好,我还有问题想请教一下
    聚类分析中,plotMDS的横纵坐标及其值分别是什么意思?横纵坐标有什么区别呢?为什么我的野生型的样品也会落在负值区域?
    pairs(log2(1+counts[,o[1:7]]), pch=”.”,lower.panel=NULL)画出的图应该怎么看呢?例如在您图里的ES66335和FE66350,横纵坐标分别表示什么?其值的大小表示什么呢?在对角线和分散之外的点各有什么区别呢?

    • Reply admin 8月 10,2015 8:13 上午

      其实我也不太明白。你可仔细阅读一下这段话:Distances on the plot can be interpreted as leading log2-fold-change, meaning the typical (root-mean-square) log2-fold-change between the samples for the genes that distinguish those samples.

  24. Reply chenzhnh 9月 16,2015 5:24 下午

    老师您好,关于DESeq2分析我有两个问题:
    1、如果手头是RPKM值(比如从TCGA上下载的RNA seq数据),能够简单粗暴的使用DESeq()这个function来分析吗?不行的话应该如何分析?
    2、您演示的DESeq分析结果中,为什么只有PH vs DE? grp中不是一共有8组吗?应该有很多个两两比较的结果才对呀?
    麻烦您了,谢谢!

    • Reply admin 9月 17,2015 9:03 上午

      感谢你的提问。

      1. 如果手头已经是RPKM值的话,不适合使用DESeq或者edgeR来分析。请使用anova算法分析。
      2. 所以只是演示啊。对于DESeq,你可以设置不同的design来满足复杂的计算要求。请先学习limma包。

      • Reply chenzhnh 9月 17,2015 10:17 上午

        非常感谢您的回答!请问:
        1、如果我手头是已经mapping好,整合成相应gene的read counts数据,可否使用DESeq分析?
        2、如果我手头是两组比较的(比如control&treatment)RPKM,可否把他们当做microarray data一样用CyberT, SAM或者 RankProduct来分析?
        3、对于您演示的文件,如果想得到两两比较的结果,先将limma包的design matrix把8组都包含进去,然后在设计不同的contrast matrix就行了是吗?谢谢您!

        • Reply admin 9月 18,2015 8:57 上午

          1. 对于DESeq和edgeR,它们的输入就是read counts.所以答案是yes
          2. 对于SAM或者RankProduct来说,不可以。对于CyberT,可以使用limma来实现。
          3. 正确。

          • Reply chenzhnh 9月 25,2015 8:11 下午

            老师您好,我仔细研究了一下您这篇文章中关于RNASeq package的教程,发现从头至尾都没有用到contrast matrix,那为何它出的结果只是默认grp PH vs DE呢?而且没有用到contrast matrix的话,我分析时想比较别的组之间的差异该如何设置呢?
            另外一个是关于edgeR教程中画barplot的问题,我搞不清楚barplot()这个函数中col这个参数的用法,为什么令col=cols[o], 其中o <- order(grp), 就可以使颜色那样显示呢?
            多谢您了!

          • Reply admin 9月 28,2015 9:22 上午

            嗯,你可以在文中搜索一下design和makeContrasts。第二个问题,你可能阅读得不太仔细。o已经重新定义了。

          • Reply chenzhnh 9月 28,2015 10:15 上午

            老师您好,我按您说的搜索了一下,的确在整个DESeq2分析过程中没有用到contrast matrix呀,就只用到了design matrix,那这个函数是如何判断比较哪两组呢?谢谢您!

          • Reply admin 10月 2,2015 11:40 上午

            我想有可能就是使用的最后一组减了第一组。

          • Reply Lixuenan 5月 15,2016 10:23 下午

            请问下博主 我要是想比较其他的两两组间比较 该如何设置呢?

          • Reply admin 5月 17,2016 5:59 上午

            你也许注意到了:con <- makeContrasts("DE-ES"=grpDE-grpES,levels=colnames(mm)) 这一句。想比较其它组,你只需要修改这一句即可。比如: con <- makeContrasts(contrasts=c('grpDE-grpES', 'grpFE-grpES'), levels=colnames(mm)) 总之,你可以把想要的比对都放到contrasts来。

  25. Reply Liping 9月 22,2015 2:16 下午

    请教博主,用mRNA-seq找差异表达基因,
    1. 看测序数据量是否足够应该怎么计算啊?是用未mapping前的fasta (或者应该是用fastq?)大小除以genome size还是cds size啊?要求覆盖多少倍呢?10倍吗?或者是用map之后的reads数目来计算呢?

    2. 如果想统计样本中有多少基因表达了,是用map到该基因的reads数,还是RPKM设cutoff呢?一般cutoff设为多少呢?
    谢谢博主

    • Reply admin 9月 25,2015 12:17 下午

      1. 是用raw reads数乘以read length阶以genome size来做的。hi-seq是要求10倍。但很多都达不到。所以有时候,我们只是简单地通过总共有多少reads来判断。对于人,小鼠等,有很多研究的文献,你可以找到相应的经验数据。
      2. 用map到的基因reads数。当然,有些时候可能会依据不同的实验目的来过滤一下低表达的基因。

  26. Reply liling 10月 16,2015 1:10 上午

    老师您好。
    edgeR的输入是raw count,用TMM方法进行样本间的normalization, 那么是否做了样本内的标准化呢?如果没有怎么做呢?可以用rpkm吗?

    • Reply liling 10月 16,2015 1:23 上午

      看到很多文章都将计算了RPKM,这个不是只是样本内的标准化吗?不需要样本间标准化吗?

      • Reply admin 10月 16,2015 12:17 下午

        再强调一下,RPKM或者FPKM是做了样品间标准化的。它是通过library size来达到样品间标准化的。

        • Reply liling 10月 19,2015 3:30 上午

          谢谢博主。不过还有些问题请教:
          edgeR如何导出TMM标准化后的表达矩阵呢?
          另外,我看到有些帖子里写,可以再TMM后再计算一下cpm或rpkm。这是不是对一套数据用了两种标准化方法呢?这种可行吗?

          • Reply admin 10月 20,2015 1:53 下午

            严格的来讲,cpm并不是一种normalization的方法,但是广意上来讲它又是的。在TMM之后再计算cpm和之前计算的cpm是不同的。因为TMM之后对library size进行了调整,有理由相信,这种调整是有意义的,所以之后再求得的cpm应该更理想一些。

    • Reply admin 10月 16,2015 12:16 下午

      样品内标准化的概念是由microarray带来的。一般而言,基因芯片因为染色不均和杂交不均,所以同一芯片内不同区间需要进行标准化。而这种技术上的障碍并不存在于RNA-seq上,所以在RNA-seq的结果比较中,不讨论样品内标准化的问题。使用RPKM并不是很好的选择。RPKM的标准化方法是在library size的基因上再对基因长度做了标准化,以达到不同基因间也可以相互比较的目的。但是,它的false positive rate会不如TMM以及DESeq2的好。

  27. Reply ping 1月 19,2016 11:46 上午

    博主好,
    使用edgeR的时候,是用classic还是glm approach,应该如何选择呢?例如我的数据为0h,3h,6h,12h;想分别得到3hvs0h,6hvs0h和12hvs0h的差异表达基因。
    另外,我还有复杂的实验设计,对照组水处理0h,3h,6h,12h,实验组药处理0h,3h,6h,12h。我想得到每个时间点(即3h,6h,12h)药处理特异表达的基因,即排除掉每个时间点对照组水处理的影响。我应该用makeContrasts,glmfit 和glmLRT来做,还是别的方法呢?

    • Reply admin 1月 19,2016 12:28 下午

      我推荐使用edgeR、或者DESeq2的GLM方法。无论多么复杂的设计,你都可以按照其说明来进行操作。在设计matrix的时候,你可以参考limma的说明文档。

      • Reply ping 1月 19,2016 7:54 下午

        多谢博主!
        例如我想比较wildtype和mutant时,推荐用edgeR的classic 还是glm方法呢?
        因为我有同学做上述case时(比较wildtype和mutant)总是用edgeR的classic approach。
        不知道博主比较过classic 和glm方法没有?结果差别大吗?

  28. Reply boya888 1月 28,2016 9:25 下午

    欧老师好,

    请问差异分析还要看log2FC的z-score,依据|z-score|>1.5确定某个基因是否差异表达是什么节奏?

    请问单个基因的z-score是咋算的呢?谢谢

    • Reply admin 1月 29,2016 12:01 下午

      使用z-score的我很少看到。不知道你是如何得到”依据|z-score|>1.5确定某个基因是否差异表达”这句的。我一般就只看p.value和FDR两个。

  29. Reply Ting 5月 3,2016 1:16 下午

    load(url(“http://qiuworld.com/bioconductor/RNA-seq/ds1.Rdata”))
    找不到这个文件了呢!qiuworld.com/bioconductor/RNA-seq/ds1.Rdata,求大神重新上传!

  30. Reply Ting 5月 9,2016 10:54 上午

    欧老师,还有个问题,如果在实验设计中,同一个样本有多个属性,例如性别和不同的温度处理,那在做差异分析的时候,怎么考虑交互作用呢?是否首先应该检验interaction,并除去具有interaction效应的这类基因呢?

  31. Reply Ting 5月 20,2016 1:32 下午

    谢谢欧老师,还想问一个问题,在过滤了低表达的基因之后,lib size会改变么?需要重新定义 d$samples$lib.size么?看这里是没有重新定义的。另外,每一个样本的lib.size是这个样本所有基因上落入的reads加起来的总和么?

    • Reply admin 5月 20,2016 5:25 下午

      按照手册来说,是需要重新定义library size的。但是通常低表达的基因对library size的贡献有限,是否重新定义library size对最后的结果的影响几乎可以忽略。但是如果你要是过滤掉了tRNA,rRNA的污染的话,library size还是最好重新定义。它的lib.size就是reads总和。

  32. Reply super 5月 31,2016 2:14 上午

    欧老师,我上次提的一个问题好像您还没解答。我做蝙蝠的基因组RNA-Seq,从Ensembl上下载蝙蝠的基因组fasta (Myotis_lucifugus.Myoluc2.0.dna.toplevel.fa)和GTF。但是一旦运行cufflinks和cuffmerge,在cuffmerge一步就报错!!!
    我的命令:
    cufflinks -p 12 -m 140 -u -g genes.gtf -o Bat_clout aligned.bam
    cuffmerge -g genes.gtf -s genome.fa assemblies.txt
    运行后就会报错:
    Error (GFaSeqGet): end coordinate (18260) cannot be larger than sequence length 18208
    Error (GFaSeqGet): end coordinate (15959) cannot be larger than sequence length 15869
    Error (GFaSeqGet): subsequence cannot be larger than 14458
    Error getting subseq for CUFF.28575.2 (5..14586)!
    [FAILED]
    我有一个解决方法,就是运行cuffmerge不要-s genome.fa选项,就能通过并生成merged.gtf,但运行cuffdiff还会报错!!!
    OK,如果我把cuffdiff的 -b genome.fa选项也省略,就能通过了……
    通过是通过了,但是我不知道
    (1)为什么会出错?
    (2)我这样的处理方法是否妥当(把cufflinks,cuffmerge和cuffdiff所有和genome.fa相关的选项都禁用)
    谢谢!

  33. Reply super 5月 31,2016 10:05 上午

    绝对的同一版本啊。比如大蝙蝠 GTF:ftp://ftp.ensembl.org/pub/release-84/gtf/pteropus_vampyrus/ Pteropus_vampyrus.pteVam1.84.gtf.gz 基因组版本: ftp://ftp.ensembl.org/pub/release-84/fasta/pteropus_vampyrus/dna/Pteropus_vampyrus.pteVam1.dna.toplevel.fa.gz。只要cuffmerge有 -s genome.fa就报错,必须把cufflinks,cuffmerge和cuffdiff所有和genome.fa相关的选项都禁用才能运行成功,只是不知道这样的处理方法是否妥当。

  34. Reply Lee 6月 6,2016 4:47 上午

    老师问一下,如果我的数据没有做生物学重复,只有对照组和实验组的分别对应5个时间点的数据,用什么方法进行差异分析比较好?

    • Reply admin 6月 13,2016 9:24 上午

      如果要上杂志的话,没有好的办法。唯一的解决之路是至少再做一组重复。如果只是预实验,可以试试edgeR或者DEseq2。

  35. Reply liu 7月 2,2016 5:46 上午

    请问版主,如何对小RNA测序数据进行差异表达分析?

  36. Reply 阳光 8月 9,2016 6:55 上午

    请问版主,已知FPKM值怎样筛选差异表达基因

  37. Reply 阳光 8月 10,2016 8:28 上午

    请问版主这几步载入数据load(url(“/bioconductor/RNA-seq/ds1.Rdata”))
    ls()
    head(counts[, 1:7], 3)
    grp <- as.factor(substr(colnames(counts), 1, 2))
    可不可以换成
    countData = read.csv( "1.csv", header=TRUE, row.names=1 )
    colData <- data.frame(condition=factor(c("R_mock_3dpi","R_mock_3dpi","R_mock_3dpi","R_mock_3dpi","R_mock_3dpi","R_mock_3dpi","R_mock_3dpi","R_mock_3dpi","R_mock_3dpi","R_mock_3dpi","R_mock_3dpi","R_mock_3dpi","R_fusarium_3dpi","R_fusarium_3dpi","R_fusarium_3dpi","R_fusarium_3dpi","R_fusarium_3dpi","R_fusarium_3dpi","R_fusarium_3dpi","R_fusarium_3dpi","R_fusarium_3dpi","R_fusarium_3dpi","R_fusarium_3dpi","R_fusarium_3dpi")),libType = c( "paired-end", "paired-end", "paired-end","paired-end","paired-end","paired-end","paired-end","paired-end","paired-end","paired-end","paired-end","paired-end","paired-end","paired-end","paired-end","paired-end","paired-end","paired-end","paired-end","paired-end","paired-end","paired-end","paired-end","paired-end") )
    因为没太看懂read counts数据是怎样读入和处理,所以不知道要不要写上condition和type

    • Reply admin 8月 11,2016 8:34 上午

      嗯,恕我直言,你的R还需要加强学习与练习,否则直接上手这样大的数据操作并不是一个非常好的主意。countData是matrix或者data.frame都可以,但mode应该保证是numeric或者integer的。colData中的condition没看出来你是从哪里得到这些的,你要用这个condition来作为DESeqDataSetFromMatrix的一部分吗?是你实验设计的一部分吗?我写成load(url(“/bioconductor/RNA-seq/ds1.Rdata”))是一个笔误,你可以写成load(url(“http://qiubio.com:8080/bioconductor/RNA-seq/ds1.Rdata”))。然后再看看你的counts是否和我的里的数据结构一致。

  38. Reply 阳光 8月 12,2016 7:43 上午

    好的,谢谢版主。我还有一个问题困扰好久了,现在手上只有RPKM/FPKM值,怎样求差异表达基因。DESeq和edgeR处理的是read counts 参照的是FC值和 P值,FDR。那么RPKM/FPKM应该怎样处理呢

  39. Reply 流星北 10月 13,2016 1:04 上午

    博主:谢谢您的分享,请教两个问题。
    1.当 dds = DESeqDataSetFromMatrix(counts, colData=data.frame(grp), design=formula(~-1+grp)),其中参数design=formula(~-1+grp)时
    dds2 <- DESeq(dds) 会报错
    错误于DESeq(dds) :
    betaPrior=TRUE can only be used if the design has an intercept.
    if specifying + 0 in the design formula, use betaPrior=FALSE
    如果修改为design=~grp时,运行正常,
    请问出错原因

    2.差异分析结果
    res
    log2 fold change (MAP): grp PH vs DE
    Wald test p-value: grp PH vs DE
    DataFrame with 30727 rows and 6 columns
    在这里是不是说差异表达的比较组是在 PH和DE两组之间的样本做比较
    那其他组(Levels: DE ES FE FG GT IS PE PH)之间做了两两比较了吗,还是差异分析是做的多组之间的分析。

    • Reply admin 10月 13,2016 2:55 下午

      1. betaPrior=TRUE时要求线性拟合时有截距,所以不能使用无截距的公式。

      2. 好问题。你可以试着自己分步算一下,看看是不是只考虑了第一组。

  40. Reply orange 11月 10,2016 5:43 上午

    欧老师,你好! 我想比较4种条件下的RNA-seq, 4中条件互相比较? 有什么好的方法做多模型之间比较吗? 过滤基因,主要是过滤掉低表达基因吗?一般如果称为低表达呢?RPKM<1 ?期待您的回答!

    • Reply admin 11月 10,2016 9:51 上午

      你有真值和假值的结果吗?如果有,可以做ROC curve来比较多个模型。
      对于过滤基因,通常来讲有两种类型,一种是做比较前的过滤,这种主要是过滤低表达值的基因。但其实它们并不影响最终的结果。对于过滤低表达值,主要是影响计算量。但对于现在的计算机来说,这点计算量不算什么。如果要过滤,低实际情况而定,不一定非是RPKM<1. 这里也没有一个定数。我一般把超过半数是零的过滤掉。一种是比较之后的过滤,这时就是使用p.value之类的过滤。也有把两头比较极端的值过滤掉的。

  41. Reply 蓝色未央 11月 14,2016 7:34 上午

    哈喽,老师您好,想请教一下,一般差异基因分析,采用DEseq软件,上下调基因的数目大概多少比较恰当,比如我卡完p值0.05,发现还有好几千,就会再卡0.01,但是这次做的时候我发现上下调都只有110左右,是不是太少了,这个数值我没有什么概念。
    O(∩_∩)O谢谢老师,可以答疑解惑

    • Reply admin 11月 14,2016 9:13 上午

      RNA-seq通常只是筛选性质的实验,最终还是看验证效果。一般情况下,使用adjust.p.value筛选得到三四千以下个差异基因是比较常见的,但并不意味着它就有实际的生物学意义。

  42. Reply 徐璐 3月 29,2017 1:26 上午

    博主你好,我想请问一下,如果利用GEO上的数据是否可以使用edgeR进行基因表达差异分析?现在我想找到肺癌中一些表达水平较高的基因,但是看了edgeR好像都需要count数据,我不知道count数据是怎么来的?每一个基因针对不同癌症是否有不同的count?我要到哪里去找count数据?
    谢谢博主

    • Reply admin 4月 7,2017 9:33 下午

      没有count number的话,edgeR和DEseq2都不合适。你可以试试limma,但是不推荐。count数据可以通过重新map然后再计数。如果你只是希望比较不同的癌症你可以试着在http://www.cbioportal.org/试试。

  43. Reply 蔡勋超 4月 20,2017 9:24 下午

    请教博主:
    博主您好,您的文章使我受益匪浅。不过我在计算离散度这一部分有一个疑问。我的理解是离散度应该主要用来衡量生物学重复间的变异情况,本博文是不是将所有重复放在一起画图了呢?R语言小白,望博主不吝解答,我要是想要做每组样品的离散度分析该怎么办呢?

    • Reply admin 4月 24,2017 7:02 上午

      是的。它的前提假设是大部分基因的表达水平都是不变的。如果是每组的话,你可以只使用它的subset来做,但是我个人认为并没有太大的意义。

  44. Reply 蔡勋超 4月 25,2017 5:20 上午

    好的,非常感谢博主

  45. Reply zsq 6月 15,2017 3:32 上午

    博主,您好。我有四个组织,每个组织三个重复。做差异分析的时候,一次性就可以输出他们之间的所有两两比较吗?design还是选择formula(~-1+grp))吗?您给的这个公式的链接我打不开,麻烦您再详细解释一下这个公式。谢谢了

    • Reply admin 6月 16,2017 6:07 下午

      我修改了链接,请再试一下吧。
      至于你的问题,在design上没有什么变化,但是,在contrast matrix上有变化,你可以问一下makeContrasts这个函数,在limma的手册里有。

  46. Reply ss 7月 14,2017 8:17 上午

    你好,我想问下,readcount列表里应该是一个基因的所有转录本的count和,还是挑里边转录本count最大的一个count作为这个基因的readcount

Leave a Reply

  

  

  

%d 博主赞过: