对于芯片数据中fold change方法和t统计量方法得比较
A comparison of fold-change and the t-statistic for microarray data analysis，这文章是大神Robert Tibshirani和Daniela M. Witten（是Robert的一个学生，现在华盛顿大学生统专业做PI）写的。整篇文章每句话都非常重要！
对于芯片数据中control和treatment样本间信号强度的比较，筛选差异大的基因，会用到fold-change或者t-statisitic。
文章作者比较了这两种情况的计算出的差异基因之间的差别。
说明在没有背景噪声的情况下用fc比较好，在有背景噪声的情况下用modified t-statisitic比较好。
作者还指出可重复性
同精确性
之间也不是完全一致。
传统的t检验在可重复性和精确性上都不如fold-change和改良的t检验，所以在芯片数据分析时，不要用传统的t检验。
Introduction
这部分，文章重点介绍了他们的实验方法，用真实的microarry数据和模拟的数据来做分析，采用的筛选基因的方法包括传统t检验，改良的t检验，两种fold-change方法。 模拟后结果说明用究竟是用fold-change或者改良的t检验取决于我们的研究是对基因表达量的绝对变化感兴趣还是对基因表达量相对于其噪声的变化感兴趣。
Therefore, a researcher’s decision to use fold-change or a modified t-statistic should be based on
biological, rather than statistical, considerations.
Statistical measures of differential expression
方法1. 传统t检验的统计量
\[T_i=\frac{\bar{x_i}-\bar{y_i}}{s_i}\]其中\(s_i\)是对于基因i的组内样本重复间的标准差。
方法2. 改良的t统计量
\[T^{"}_{i}= \frac{\bar{x_i}-\bar{y_i}}{s_i+s_0}\]其中\(s_0\)是为了让\(T^{"}_{i}\)的变异系数最小的一个常量，在本文中用的是Significance Analysis of Microarrays (Tusher et al. 2001)中的计算方法。 我在wiki1上查了一下，上面说\(s_0\)是根据\(\alpha\)分位数来定的。另外有些文章2中直接将\(s_0\)定义成\(s_i\)的中位数。
方法3. 标准fold-change
\[FC_i=\frac{\bar{x}^{'}_{i\cdot}}{\bar{y}^{'}_{i\cdot}}\]其中\(\bar{x^{'}_{ij}}\)和\(\bar{y^{'}_{ij}}\)是基因i在组内样本j的原始表达量（control除以treatment）。
方法4. 另一种fold-change算法
\[FC_i=\bar{x}^{'}_{i\cdot}-\bar{y}^{'}_{i\cdot}\]其中我们标记方法3的fold-change为\(FC_{ratio}\)，方法4中的为\(FC_{difference}\)。
另外，我们可以看出在增加\(s_0\)的大小后，方法2的排序结果会逼近方法4的排序结果。 \($s_0\)这个常数越大，t统计量的分母越一致，那么对基因排序起到关键作用的就是分子，分子又通方法4的一样…..所以方法2和方法4的结果会比较一致。
为了证明高度可重复的统计量并不是十分准确的，他们在文中构建了一个人工统计量：
\[P_i=(\bar{x_{i\cdot}})^3-(\bar{y_{i\cdot}})^3\]P代表power，文中说不建议在实际中使用这个人工统计量。
Overview of the genes selected using the different measures
文中分别用模拟的数据，按上面的4个方法来计算差异基因。
从下图中可以看出，按普通的t统计量挑选的基因的标准差都很小， fold-change方法找到的基因在control和treatment里有很大的差异， 改良后的t统计量找到的基因有较小的标准差和组间较大的差异。 由于普通的t统计量和fold-change的结果完全不一样，对于研究者来说，要考虑某个基因在研究中是否重要的根据是表达量的偏移还是标准差的偏移。
Simulated and real data Concordance
这节主要介绍了他们用的两种模拟策略：一种是基于正态分布的偏移，另一种是基于样本分组间均值和方差的偏移。
关于一致(Concordance)的定义：在一个基因（挑选后的）集合中的基因也在另一个基因集合中的比例。
Analysis of the accuracy of the different measures
下图就是他们分析的一个结果，可以看到样本量在红线左侧时，用\(FC_{difference}\)比较好，样本量很大时就用改良的t统计量较好。
其实从我们用的方法来寻找差异基因，一共就度量两类差异： 第一类是\(\mid \mu_{control} - \mu_{treatment} \mid\) ，另一类是\(\mid \frac{\mu_{control} - \mu_{treatment}}{\sigma}\mid\)。 后者主要是标准化基因间的噪声（方差），因为基因均值之间的差异可能是源于基因间表达量的方差有差异。
接下来，他们又考虑了一个比较极端的例子。
| Control | Treatment | FCdifference | FCratio | T | |
|---|---|---|---|---|---|
| Gene1 | 150, 200, 250 | 1, 50, 100 | 3.51 | 3.97 | 1.69 |
| Gene2 | 101.1, 101.2, 101.3 | 100.1, 100.2, 100.3 | 0.014 | 1.01 | 12.25 |
这个例子中Gene1表达量在两组间差异很大，Gene2表达量在两组间差异很小。但是如果用t统计量来衡量的话，Gene2远大于Gene1。
我自己计算的结果同论文上的不太一样，Gene1的t统计量是3.6844，而不是1.69。
t.test(c(101.1,101.2,101.3),c(100.1,100.2,100.3))
Welch Two Sample t-test
data: c(101.1, 101.2, 101.3) and c(100.1, 100.2, 100.3)
t = 12.247, df = 4, p-value = 0.0002552
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.7733042 1.2266958
sample estimates:
mean of x mean of y
101.2 100.2
t.test(c(150,200,250),c(1,50,100))
Welch Two Sample t-test
data: c(150, 200, 250) and c(1, 50, 100)
t = 3.6844, df = 3.9996, p-value = 0.02113
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
36.87866 262.45467
sample estimates:
mean of x mean of y
200.00000 50.33333
那么在生物学家的眼中，究竟哪个基因是需要关注的基因呢？我觉得应该是Gene1吧。
From this perspective, the question of whether the fold-changes or a modified t-statistic results in more accurate gene orderings is really a biological one, rather than a statistical one, as it depends on what types of expression differences between control and treatment have biological relevance.
Analysis of the reproducibility of the different measures
接下来他们又检测了各种方法的可重复性。 从下面的两张图中可以看出前面自己定义的P统计量有很好的重复性（下面第一个图），但是没有很好的精确性（下面第二个图）。 因而说明了，重复性和精确性在某些检验指标中不可间得的问题。
Conclusions
- \(FC_{difference}\)和改良的t检验都是t检验的一种改良形式，也就是有不同的\(s_0\)。一些\(s_0\)的选取可以提高精确性。
- 别用普通的t检验。
- 可重复性高并不暗示着精确性高，The issues of reproducibility and accuracy should be kept separate when evaluating the performance of a statistic.
- 在实际分析中并没有FC和t统计量谁好谁坏的说法，都要看生物学意义。有噪声干扰就用改良的t检验，没有噪声干扰就用FC。
参考资料
新博客主题：书签
在网站首页加了一个新的连接：书签，是我做的类似于Product Hunt的博客主题， 还没有做外观的设计，主要实现的是每个post的连接都指向一个外部网站。 这个子页面主要用来记录我随手在网上看到的有用信息，以及我对这个信息的简要评论， 如果今后有时间对这个信息写了长篇博客，那么这个连接就会在书签中消失。
度量距离时对数据做不做比例归一化(scale)
当我们处理数据时时常会遇到是否要对数据进行比例调整(scale)的问题，那么究竟应不应该做比例调整要取决于数据的实际含义。
前一段时间为了弄明白Gower Distance在网上查找了一个slides1，里面介绍了很多常用距离，并且简单解释了数据比例尺度调整的问题。
- 对于4个人的年龄和身高，有数据
| Person | Age [years] | Height [cm] |
|---|---|---|
| A | 35 | 190 |
| B | 40 | 190 |
| C | 35 | 160 |
| D | 40 | 160 |
画图可以看到A、B比较近，C、D比较近。
有些地区喜欢用feet来作为身高的度量单位，如果换成feet，数据会变成
| Person | Age [years] | Height [feet] |
|---|---|---|
| A | 35 | 6.232 |
| B | 40 | 6.232 |
| C | 35 | 5.248 |
| D | 40 | 5.248 |
画图的话可以看到此时A、C比较近，B、D比较近。
那么究竟哪两个人的数据比较接近呢？
我们来做一下scale,
| Person | Age [scaled] | Height [scaled] |
|---|---|---|
| A | -0.87 | 0.87 |
| B | 0.87 | 0.87 |
| C | -0.87 | -0.87 |
| D | 0.87 | -0.87 |
结果发现这四个人距离差不多，分不出子类。
- 在来看另一种情况
| Object | x1 | x2 |
|---|---|---|
| A | 13.3 | 38.0 |
| B | 12.4 | 45.4 |
| C | -122.7 | 45.6 |
| D | -122.4 | 37.7 |
有四个观测，分别知道它们的变量x1和变量x2数值，在R中scale(dat) 会发现四个观测分散很远，如果直接画图，就发现其实A、B距离近，C、D距离远。
如果x1和x2分别代表经度和纬度，那么这个数据就不应该标准化，A、B两个地点本来就是距离近，标准化后它本身的特点就不存在了。
到底用不用归一化呢？
1.做不做归一化，要知道
- 变量取值范围大，这个变量就在计算距离时权重大
- 距离的远近是由归一化后的数值决定的，不同的归一化，最后求出的距离也不一样
- 归一化对每个变量赋予同样大的权重
- 另一种可行方法是重赋值权重
2.这些情况下必须归一化
- 变量单位不同
- 我们自己期望属于要有相同的权重
3.这些情况下不要归一化
- 变量单位相同
4.一般情况下
- 请归一化
Reference
-
https://stat.ethz.ch/education/semesters/ss2012/ams/slides/v4.2.pdf ↩
看过且收藏的一些电影
整理硬盘，要删除点东西，一看电影都这么占空间，所以就决定删了电影，在删之前先小结一下。
说收藏其实不太准确，因为在国内看电影绝大多数是网络盗版，没有付费，收藏也不是正经的买碟，下载了别人上传的内容。
美丽心灵
今年最唏嘘的事情就是纳什在领完奖回家时做出租车被撞死了。
在看美丽心灵这部片子时，感受最大的特点是——完全没写什么事实。我是先看完纪录片和别人的记事报道，才看的这部电影，看了两三次，才完全看完。
在光辉的外表下，纳什只是一个（自愈）的精神病人。当然，更痛苦的是他的二儿子，也得了这种病，并且完全没有好的迹象。
具体的可以去看玑衡写的《我所认识的约翰纳什》，以及纳什的纪录片。
我特别有感触的是在纪录片里，他妻子说自己的儿子也有这种病（需要人照顾，不能完全自理），但是现在她还能照顾孩子，如果纳什和她不在了，她的孩子会怎么样？
结果，她和纳什就一起突然离世。
美丽心灵是一部奥斯卡获奖影片，但是它内容描写的太美好，而现实往往是非常残酷的。一人得精神疾病，会影响整个家庭，不可逆的。
并且不是说数学家都是mad and crazy，中国宣传了陈景润之后，数学家的形象明显都变味了，结果美国宣传了一个“疯子”，这数学家的形象什么时后才能在人们的心中变得正常起来。-_-|||
还有就是，没机会在他活着时见面了，克里克也早没了，现在还想见的就只剩下沃森和Le Guin，不知道今后有没有机会。
我想亲自见面的前辈贤者都在不断变成先贤，而我也在不断变老，就是这样，完成于2015年还剩下2个月的10月31日。
黑客帝国
硬盘里的黑客帝国三部曲我看了又看，感觉这是我看过的最好看的片子（说道这里，我想到了前几年特火的阿凡达，我完全没觉得它做的又多好，各种外星生物和殖民的创意感觉游戏中都已经出现过了）。
黑客帝国的逻辑构思是这部片子最大的亮点，第一次看之前，我只看过它的几个动画版，所以对整部剧没有任何了解，看了正片之后觉得这真是厉害：人被当作燃料，思维进入计算机的世界继续存活，所有人都变成了一串代码，由更高级的计算机人工智能系统控制，程序里有bug，会导致异常的事件和人物的出现，计算机人工智能系统为了维护统治，而不断的捕杀这些想知道真相的人。
这部片子引出了一句经典的网络用语：脑后插管
。看完后估计大家都会思考是不是我们现在所存在的空间就是一个Matrix。
也就是宇宙。当然这个说法是毫无科学根据的。
说个体外话，我觉得逃出太阳系没准是未来人们必须要完成的生存任务。
那么宇宙的边界，宇宙大爆炸前没有时间没有空间，爆炸后延伸出来的物质究竟是向哪里扩散？
扩散的边界会是什么样子？理解这些问题会帮助人类向太阳系外前进。
另外这个三部曲，我实在是舍不得删除。;)
但是我对导演兄弟/兄妹/姐妹 实在是理解不能_(:з」∠)_。
降世神通：最后的气宗
我没看过动画，偶然间看的这部片子，特效还不错（我看的少，所以评价一般都很高），里面的武术动作也比较精彩。 但是网上整体评价不高，我也不知道怎么回事，我就看个热闹。
基因表达与spiked-in
Revisiting Global Gene Expression Analysis这篇文章在刚发表的时候我们实验室就研读过，但是当时由于自己关注不够多，并且自己处理的数据进行了比例转换所以不涉及这个问题。
但是前一段时间又看了一遍这篇文章，我们都知道在single cell RNA表达分析中一般都要用spiked-in来解决normalization的问题。
但是也有很多单细胞表达量的分析没有做这个分析。
由于同第一次看这篇文章已经有很久了，我早已忘记spiked-in主要应用在肿瘤的研究中，所以现在觉得只要是肿瘤的研究都要做这个分析，才能准确确定表达量。
今天我开始有时间好好读一遍这篇文章，发现文章中提到的表达量normalization用到的两个样本high c-Myc和low c-Myc都是同样的细胞数量，所以对于bulk数据，这个方法好像没法用。
Encode已经把spiked-in当作标准流程了，新上传的RNA表达量数据有些已经带有spiked-in了，我需要在去检查一下究竟那些数据是单细胞（定量细胞的）还是bulk（组织）的大致看了一下，好像tissue的就有。
现在简单复述一下这篇文章的内容。
- 首先以往研究中发现c-Myc高的细胞中很多基因的表达水平整体提高，整体的RNA水平要比c-Myc低的细胞高出2-3倍。 这个现象引出问题：以往的标准化方法没有考虑到整体表达水平的升高和抑制，这样会导致解读RNA表达量数据出现问题。
-
如上图所示，AB图中两个细胞的RNA水平一致，所以标准化后可以看到基因B和E表达量是明显升高的。 对于CD图中的两个细胞中的一个细胞表达了比另一个细胞多1到2倍的RNA，如果还用正常的标准化方法可以看到A、G、I这三个基因表达量上升。 而D、E、F这三个基因表达量下降。但实际情况呢，其实这些基因基本上表达量都升高了。
这个标准化产生的问题是基于这样一个假设，即我们认为每个细胞的所有mRNA表达水平是一致的。
-
进一步说明我们常用的分析方法将这种表达量的差异看作技术误差，也就是理解为噪声， 并希望在研究中对于不同的样本或实验之间的
表达水平
要有同样的中位数，或均值，或在一个范围内的表达量分布要基本一致。 -
为了得到可靠的基因表达量，文章中采用spiked-in 标准，加入表达量基本确定的RNA作为参照。 他们分别在Microarrays，RNA-seq以及Nanostring中做了实验，加入spiked-in RNA的表达谱变化更接近真实情况。 另外，做实验时严格统计了c-Myc高和c-Myc低的细胞数目1在同等细胞数目的条件下进行的RNA表达情况的分析。
-
当无法具体统计细胞数目时怎么办？ When cell counting may be problematic, as for expression experiments from solid tumors or tissues, DNA content may be used as a surrogate if ploidy and DNA replication profiles are also characterized to prevent the introduction of a DNA content-based artifact.
-
以前的全基因组表达量数目中有多少我们已经解读错了？ How prevalent is misinterpretation of genome-wide expression data due to the assumption that cells produce similar levels of total RNA? The answer is likely related to the prevalence of regulatory mechanisms that globally amplify or suppress transcription.
最后对于RPKM的标准化可以用R包affy中的loess.normalize2来实现。说白了就是做回归平滑的时后可以选取所有的样本点（老方法），或者只选取一部分样本点（spiked-in的RNA，新方法）。
备注：
