01

自学miRNA-seq分析第七讲~miRNA样本配对mRNA表达量获取

这一讲其实算不上是自学miRNA-seq分析,本质就是affymetrix的mRNA表达芯片数据分析,而且还是最常用的那种GPL570    HG-U133_Plus_2,但是因为是跟miRNA样本配对检测的,而且后面会利用到这两个数据分析结果来做共表达网络分析等等,所以就贴出对该芯片数据的分析结果。文章里面也提到了 Messenger RNA expression analysis identified 731 probe sets with significant differential expression,作者挑选的差异分析结果的显著基因列表如下: Continue reading

01

自学miRNA-seq分析第六讲~miRNA表达量差异分析

这一讲是miRNA-seq数据分析的分水岭,前面的5讲说的是读文献下载数据比对然后计算表达量,属于常规的流程分析,一般在公司测序之后都可以拿到分析结果,或者文献也会给出下载结果。但是单纯的分析一个样本意义不大,一般来说,我们做研究都是针对于不同状态下的miRNA表达量差异分析,然后做注释,功能分析,网络分析,这才是重点,也是难点。我这里就直接拿文献处理好的miRNA表达量来展示如何做下游分析,首先就是差异分析啦: Continue reading

05

用samr包对芯片数据做差异分析

本来搞差异分析的工具和包就一大堆了,而且limma那个包已经非常完善了,我是不准备再讲这个的,正好有个同学问了一下这个包,我就随手测试了一下,顺便看看它跟limma有什么差异没有!手痒了就记录了测试流程!

学习一个包其实非常简单,就是找到包的官网看看说明书即可!说明书链接

 

Continue reading

09

差异分析是否需要比较矩阵

最流行的差异分析软件就是limma了,它现在更新了一个voom的算法,所以既可以对芯片数据,也可以对转录组高通量测序数据进行分析,其它所有的差异分析软件其实都是模仿这个的。

我以前讲到过做差异分析,需要三个数据:

  • 表达矩阵
  • 分组矩阵
  • 差异比较矩阵

前面两个肯定是必须的,有表达矩阵,样本必须进行分组,才能分析,但是我看到过好几种例子,有的有差异比较矩阵,有的没有。

后来我仔细研究了一下limma包的说明书,发现这其实是一个很简单的问题。

大家仔细观察下面的两个代码

首先是不需要差异比较矩阵的

    library(CLL)
    data(sCLLex)
    library(limma)
    design=model.matrix(~factor(sCLLex$Disease))
    fit=lmFit(sCLLex,design)
    fit=eBayes(fit)
    options(digits = 4)
    #topTable(fit,coef=2,adjust='BH') 
    > topTable(fit,coef=2,adjust='BH')
               logFC AveExpr      t   P.Value adj.P.Val     B
    39400_at  1.0285   5.621  5.836 8.341e-06   0.03344 3.234
    36131_at -0.9888   9.954 -5.772 9.668e-06   0.03344 3.117
    33791_at -1.8302   6.951 -5.736 1.049e-05   0.03344 3.052
    1303_at   1.3836   4.463  5.732 1.060e-05   0.03344 3.044
    36122_at -0.7801   7.260 -5.141 4.206e-05   0.10619 1.935
    36939_at -2.5472   6.915 -5.038 5.362e-05   0.11283 1.737
    41398_at  0.5187   7.602  4.879 7.824e-05   0.11520 1.428
    32599_at  0.8544   5.746  4.859 8.207e-05   0.11520 1.389
    36129_at  0.9161   8.209  4.859 8.212e-05   0.11520 1.389
    37636_at -1.6868   5.697 -4.804 9.355e-05   0.11811 1.282

然后是需要差异比较矩阵的

    library(CLL)
    data(sCLLex)
    library(limma)
    design=model.matrix(~0+factor(sCLLex$Disease))
    colnames(design)=c('progres','stable')
    fit=lmFit(sCLLex,design)
    cont.matrix=makeContrasts('progres-stable',levels = design)
    fit2=contrasts.fit(fit,cont.matrix)
    fit2=eBayes(fit2)
    options(digits = 4)
    topTable(fit2,adjust='BH')

               logFC AveExpr      t   P.Value adj.P.Val     B
    39400_at -1.0285   5.621 -5.836 8.341e-06   0.03344 3.234
    36131_at  0.9888   9.954  5.772 9.668e-06   0.03344 3.117
    33791_at  1.8302   6.951  5.736 1.049e-05   0.03344 3.052
    1303_at  -1.3836   4.463 -5.732 1.060e-05   0.03344 3.044
    36122_at  0.7801   7.260  5.141 4.206e-05   0.10619 1.935
    36939_at  2.5472   6.915  5.038 5.362e-05   0.11283 1.737
    41398_at -0.5187   7.602 -4.879 7.824e-05   0.11520 1.428
    32599_at -0.8544   5.746 -4.859 8.207e-05   0.11520 1.389
    36129_at -0.9161   8.209 -4.859 8.212e-05   0.11520 1.389
    37636_at  1.6868   5.697  4.804 9.355e-05   0.11811 1.282

大家运行一下这些代码就知道,两者结果是一模一样的。

而差异比较矩阵的需要与否,主要看分组矩阵如何制作的!

design=model.matrix(~factor(sCLLex$Disease))

design=model.matrix(~0+factor(sCLLex$Disease))

有本质的区别!!!

前面那种方法已经把需要比较的组做出到了一列,需要比较多次,就有多少列,第一列是截距不需要考虑,第二列开始往后用coef这个参数可以把差异分析结果一个个提取出来。

而后面那种方法,仅仅是分组而已,组之间需要如何比较,需要自己再制作差异比较矩阵,通过makeContrasts函数来控制如何比较!

22

用RankComp的思想来做差异基因分析

是福建医科大学的学者开发的,
文章里面详细讲解了他们的这个差异分析的统计学原理
大意就是找到同一组织的normal样本的表达量数据,几百个,这样就可以分析2万基因之间的互相配对,检测表达量是否在几百个样本里面稳定的不一样!

我现在还不是很确定这个方法,只是试一试,欢迎与我交流对该方法的讨论!

文章是:

Wang H, Sun Q, Zhao W, et al. Individual-level analysis of differential expression of genes and pathways for personalized medicine[J]. Bioinformatics, 2014: btu522.

他们把它写成了一个R包,可以下载使用,但是必须用R2.15.2版本,我用了一下,不好用!

We can download the R code for in http://bioinformatics.oxfordjournals.org/content/31/1/62/suppl/DC1

他们这个程序真心不好用,但是很容易看懂算法,可以自己用R语言写一个来实现同样的过程!

比如A基因在几百个样本里面表达量都是3左右,而B基因都是5左右,而且满足99%的A表达量高于B,那么这就是一个稳定的基因对!
一般2万基因之间可以配成2亿个基因对,其中稳定的大概有10%~40%
然后我们对每个疾病样本都可以进行检验,看看这样稳定的基因对是否被改变!
比如,我们拿到一个疾病样本的2万个基因的表达量,我们挑取一个,如果它有100个稳定的up的基因对,100个稳定的down的基因对
那么,我看看这些基因对是否被改变了,如果这样还有70基因对在该疾病样本里面仍然是up,60个是down,那么我用Fisher精确检验的结果是
4
这个基因在该疾病样本,相对于normal pool并没有差异表达!当然检验得到的P值最后可以做FDR校验。
依次这样,把所有的gene都分析完,就知道这个样本有哪些差异的gene了。
介绍完原理,我们拿一个具体的例子来看看吧:
首先我们下载一个2008年的一个人的肝脏表达数据(Gene Expression in Human Liver),都是正常组织,共427个样本。
不过这个芯片比较小众,是默克医药公司定制化的, 需要下载探针对应gene的文件!
我们读取GSE9588这个数据,得到表达矩阵,然后计算rank矩阵,然后计算得到comp矩阵
5
> table(rank_comp)
rank_comp
     down        no        up
    58479 465752098     58479
>
不知道为什么这个数据,stable的那么少,不知道是不是出了什么问题!
其实我的程序都是对的了,只是因为这个数据集已经不是纯粹的表达量数据了,而是这427个样本的数据都减去了某个样本的表达量。
这样每个个体的基因之间的表达量排序就会被干扰,导致得到的稳定基因对非常少!!!
但是,我后来下载了GTEx的表达数据,拿那里面的normal组织样本表达量来做,可以得到非常多的稳定基因对。
实际代码大概是:
得到正常组织的表达矩阵:
然后计算表达矩阵的rank,得到各个样本自己的基因排序情况,得到排序矩阵!
处理排序矩阵,每个基因对之间都算一下是否稳定,得到稳定性描述矩阵!
然后根据每个疾病个体的基因表达情况,来循环每个基因, 看看该基因是否差异!

 

十二 11

关于limma包差异分析结果的logFC解释

首先,我们要明白,limma接受的输入参数就是一个表达矩阵,而且是log后的表达矩阵(以2为底)。

那么最后计算得到的logFC这一列的值,其实就是输入的表达矩阵中case一组的平均表达量减去control一组的平均表达量的值,那么就会有正负之分,代表了case相当于control组来说,该基因是上调还是下调。

我之前总是有疑问,明明是case一组的平均表达量和control一组的平均表达量差值呀,跟log foldchange没有什么关系呀。

后来,我终于想通了,因为我们输入的是log后的表达矩阵,那么case一组的平均表达量和control一组的平均表达量都是log了的,那么它们的差值其实就是log的foldchange

首先,我们要理解foldchange的意义,如果case是平均表达量是8,control是2,那么foldchange就是4,logFC就是2咯

那么在limma包里面,输入的时候case的平均表达量被log后是3,control是1,那么差值是2,就是说logFC就是2。

这不是巧合,只是一个很简单的数学公式log(x/y)=log(x)-log(y)

 

十二 11

用excel表格做差异分析

其实主要要讲的不是用excel来做差异分析,只是想讲清楚差异分析的原理,用excel可视化的操作可能会更方便理解,而且想告诉大家,其实生物信息学分析,本来就很简单的,那么多软件,只有你理解了原理,你自己就能写出来的!

首先,还是得到表达矩阵,下面绿色的样本是NASH组,蓝色的样本是normal组

image001

我们进行差异分析,很简单,就是看两组的表达值,是否差异,而检验的方法就是T检验。

=AVERAGE(D2:L2)    ##求NASH组的平均表达量

=AVERAGE(M2:S2)    ###求normal的平均表达量

=T2-U2             ##计算得到logFOLDchange值

=AVERAGE(D2:S2)    ###得到所有样本的平均表达量

=T.TEST(D2:L2,M2:T2,2,3)  ###用T检验得到两个组的表达量的差异显著程度。

简单检查几个值就可以看到跟limma包得到的结果差不多。

image002

 

十二 11

用limma包对芯片数据做差异分析

下载该R语言包,然后看说明书,需要自己做好三个数据(表达矩阵,分组矩阵,差异比较矩阵),总共三个步骤(lmFit,eBayes,topTable)就可以啦

image001

首先做第一个数据,基因表达矩阵!

自己在NCBI里面可以查到下载地址,然后用R语言读取即可

exprSet=read.table("GSE63067_series_matrix.txt.gz",comment.char = "!",stringsAsFactors=F,header=T)

rownames(exprSet)=exprSet[,1]

exprSet=exprSet[,-1]

image002

然后做好分组矩阵,如下

image003

然后做好,差异比较矩阵,就是说明你想把那些组拿起来做差异分析,如下

image004

最后输出结果:

我进行了6次比较,所以会输出6次比较结果

image005

最后打开差异结果,解读,说明书如下!

 

 

image006

在我的github有完整代码