对于实验生物学研究者来说,经常遇到的一个场景是,对照组测量3个值,实验组测量3个值,然后用T检验得出均值是否有显著性差异,然后做出结论。如下图所示:

那么T检验靠谱吗?本文采用模拟抽样的方法来探讨这件事情,由于本人统计学知识有限,还请各位读者批评指正。

我们通过构建两个随机总体,总体1的均值固定为1(模拟control),总体2的均值我们从依次取0.5到2.0(模拟treatment相对与control的变化),标准差我们取0.01、0.05、0.10(模拟测量误差1%,5%,10%)。总体数量为1000,每次随机抽取3个样本,抽样1000次,考察通过T检验得出正确结论的频率。代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
import numpy as np
import pandas as pd
import seaborn as sns
from scipy import stats
import matplotlib.pyplot as plt
plt.xkcd()
sns.set(style="white", context="talk")
def norm_t(loc1, loc2, scale1, scale2, n):
normSpace1 = stats.norm.rvs(loc = loc1, scale = scale1, size = 1000)
normSpace2 = stats.norm.rvs(loc = loc2, scale = scale2, size = 1000)
count = 0
for i in range(1000):
sample1 = np.random.choice(normSpace1, n)
sample2 = np.random.choice(normSpace2, n)
t, p = stats.ttest_ind(sample1, sample2)
if (p < 0.05 and (sample1.mean() > sample2.mean()) == (loc1 > loc2)):
count = count + 1
return count / 1000
def plot(stds):
df = pd.DataFrame()
for std in stds:
x = np.arange(0.5, 2.1, 0.1)
y = [norm_t(loc1 = 1, loc2 = i, scale1 = std, scale2 = std, n = 3) for i in x]
z = [std for i in x]
df_std = pd.DataFrame([x, y, z]).T
df_std.columns = ["difference", "sensitivity", "std"]
if df.shape == (0, 0):
df = df_std
else:
df = df.append(df_std)
sns.pointplot(x = "difference", y = "sensitivity", hue = "std", data = df)
plot([0.01, 0.05, 0.10])

通过上图,我们发现,当标准差我们取0.01(模拟测量误差1%)时,增大到1.1倍或者减小到0.9倍,正确率可以达到100%。而当标准差取0.05和0.10(模拟测量误差5%和10%),效果就不理想了,本来是有差异的,还是有很大几率通过T检验得出错误的结论。显而易见的是,随着标准差的增加,T检验的效果变差了。

那么当标准差取0.10时,随着抽样样本数的增加,T检验判断正确的频率会是怎么样的呢?请看如下代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
def plot2(diffs):
df = pd.DataFrame()
for diff in diffs:
x = list(range(3, 60, 3))
y = [norm_t(loc1 = 1, loc2 =1.1, scale1 = 0.10, scale2 = 0.10, n = i) for i in x]
z = [diff for i in x]
df_diff = pd.DataFrame([x, y, z]).T
df_diff.columns = ["samle number", "sensitivity", "difference"]
if df.shape == (0, 0):
df = df_diff
else:
df = df.append(df_diff)
sns.pointplot(x = "samle number", y = "sensitivity", hue = "difference", data = df)
plot2([0.5, 0.8, 1.2, 2.0])

通过上图,我们发现,即使标准差达到0.10,当样本数量达到30以上时,T检验的效果就非常好了。

那么如果两个总体没有差异,T检验得出有差异的错误结论的情况是什么样子的呢?我们构建两个均值为1容量为1000的总体,每次随机抽取3个样本,在不同标准差及样本数目的情况下,考察T检验得出正确结论的频率。请看代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
def norm_t2(loc, scale1, scale2, n):
normSpace1 = stats.norm.rvs(loc = loc, scale = scale1, size = 1000)
normSpace2 = stats.norm.rvs(loc = loc, scale = scale2, size = 1000)
count = 0
for i in range(1000):
sample1 = np.random.choice(normSpace1, n)
sample2 = np.random.choice(normSpace2, n)
t, p = stats.ttest_ind(sample1, sample2)
if p > 0.05:
count = count + 1
return count / 1000
def plot3(stds):
stds = [0.1, 0.5, 1]
df = pd.DataFrame()
for std in stds:
x = list(range(3, 60, 3))
y = [norm_t2(loc = 1, scale1 = std , scale2 = std, n = i) for i in x]
z = [std for i in x]
df_std = pd.DataFrame([x, y, z]).T
df_std.columns = ["samle number", "specificity", "std"]
if df.shape == (0, 0):
df = df_std
else:
df = df.append(df_std)
sns.pointplot(x = "samle number", y = "specificity", hue = "std", data = df)
plt.ylim(0, 1.2)
plot3([0.1, 0.5, 1])

通过上图,我们惊讶地发现,当两个总体没有差异时,不管标准差是多少和取样数目是多少,T检验都有非常好的表现,即得到差异不显著的的结论。

这里构建的总体都是正太总体,其实对于其他分布的总体,结论是一样的,这里就不展示了。

通过上面的探索,我们可以得到如下结论:

  1. 对于两个未知总体,通过T检验考察均值是否有显著差异,如果得出结论是差异不显著,那么进一步分析这些数据的标准差是否太大了,考虑是否增加抽样样本数做进一步分析。
  2. 对于两个未知总体,通过T检验考察均值是否有显著差异,如果得出结论是差异显著,那么请相信它吧!

如果差异不显著,我们通过增大样本数量,使得差异显著;如果差异显著,那就是一个好结果嘛!
哈哈,T检验是实验生物学家的利器啊!

FASTA

FASTA文件是DNA/RNA序列文件。第一行是由大于符号>打头的任意文字说明,主要为标记序列用。从第二行开始是序列本身,标准核苷酸符号或氨基酸单字母符号。通常核苷酸符号大小写均可,而氨基酸一般用大写字母。文件中和每一行都不要超过80个字符(通常60个字符)。

1
2
>sequence info
GATGGGATTGGGGTTTTCCCCTCCCATGTGCTCAAGACTGGCGCTAAAAGTTTTGAGCTTCTCAAAAGTC

FASTQ

FASTQ格式的序列一般都包含有四行,第一行由’@’开始,后面跟着序列的描述信息,这点跟FASTA格式是一样的。第二行是序列。第三行由’+’开始,后面也可以跟着序列的描述信息。第四行是第二行序列的质量评价,字符数跟第二行的序列是相等的。

1
2
3
4
@SEQ_ID
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT
+
!”*((((***+))%%%++)(%%%%).1***-+*”))**55CCF>>>>>>CCCCCCC65

BED

BED文件是一种基因组注释文件,或者叫做基因组坐标文件。每行表示一段序列,只不过不是ATCG碱基,而是它们在基因组上的坐标。每一行有由tab分割的多列,BED行有三个必须的列和九个额外可选的列。 每行的数据格式要求一致。
3个必须列:

  1. chrom,染色体或scafflold 的名字
  2. chromStart,序列在染色体或scaffold的起始位置
  3. chromEnd,序列在染色体或scaffold的结束位置

9个可选列:
4. name 指定BED行的名字,这个名字标签会展示在基因组浏览器中的bed行的左侧。
5. score 0到1000的分值, 可以设置为’0’。
6. strand 定义链的方向,”+” 或者”-”
7. thickStart 起始位置(The starting position at which the feature is drawn thickly)(例如,基因起始编码位置)
8. thickEnd 终止位置(The ending position at which the feature is drawn thickly)(例如:基因终止编码位置)
9. itemRGB 是一个RGB值的形式, R, G, B (eg. 255, 0, 0), 如果itemRgb设置为’On”, 这个RBG值将决定数据的显示的颜色。可以设置为’0,0,0’
10. blockCount BED行中的block数目,也就是外显子数目。
11. blockSize 用逗号分割的外显子的大小, 这个item的数目对应于BlockCount的数目。
12. blockStarts 用逗号分割的列表, 所有外显子的起始位置,数目也与blockCount数目对应。

GFF/GTF

GFF/GTF文件和BED文件一样,也是一种基因组注释文件。每行表示一段序列,有9个tab分割的列共同描述这段序列的信息。
9个列定义如下:

  1. seq_id:序列的编号,一般为chr或者scanfold编号;
  2. source: 注释的来源,一般为数据库或者注释的机构,如果未知,则用点“.”代替;
  3. type: 注释信息的类型,比如Gene、cDNA、mRNA、CDS等
  4. start:该基因或转录本在参考序列上的起始位置;
  5. end: 该基因或转录本在参考序列上的终止位置;
  6. score: 得分,数字,是注释信息可能性的说明,可以是序列相似性比对时的E-values值或者基因预测是的P-values值,“.”表示为空;
  7. strand: 该基因或转录本位于参考序列的正链(+)或负链(-)上;
  8. phase: 仅对注释类型为“CDS”有效,表示起始编码的位置,有效值为0、1、2(对于编码蛋白质的CDS来说,本列指定下一个密码子开始的位置。每3个核苷酸翻译一个氨基酸,从0开始,CDS的起始位置,除以3,余数就是这个值,,表示到达下一个密码子需要跳过的碱基个数。该编码区第一个密码子的位置,取值0,1,2。0表示该编码框的第一个密码子第一个碱基位于其5’末端;1表示该编码框的第一个密码子的第一个碱基位于该编码区外;2表示该编码框的第一个密码子的第一、二个碱基位于该编码区外;如果Feature为CDS时,必须指明具体值。);
  9. attributes:一个包含众多属性的列表,格式为“标签=值”(tag=value),标签与值之间以空格分开,且每个特征之后都要有分号;(包括最后一个特征),其内容必须包括gene_id和transcript_id。以多个键值对组成的注释信息描述,键与值之间用”=”,不同的键值用”;”。

有空继续添加…

读本科那会儿,写博客、写日志、弄QQ空间的人很多。我的一位大学同学,坚持写了4年博客,我也看了4年。这期间,微博流行起来了,再后来又是微信,进入了移动互联网的时代。大家喜欢阅读短小的文字,喜欢在手机上阅读,博客似乎被人遗忘了。大学毕业后,这位同学的博客停止更新了。有时候,做一件事情是因为身边有人在做同样的事情,我写博客的源头或许是这位同学做了这样的事情。

2013年6月21日,注册了域名BioChen.com,用WordPress搭建了博客。我不是一个善于写作的人,博客以分享自己学习和工作的经验和感悟为主。跟预想的一样,最开始的一两年没什么访问量。随着文章越来越多,访问量慢慢地有所提升,到现在日访问量已经过千。一些网友的也会给我发邮件,寻求进一步的帮助;我也收到过几次网友的打赏。这些让我觉得写博客是有意义的,能帮助到一些人。

2020年3月17日,我决定放弃WordPress,改用Hexo;并把博客迁移到云主机,启用新的网址www.biochen.org/blog/cn。原博客网址blog.biochen.com还会保留一段时间。WordPress越来越臃肿,它的编辑器我也不喜欢;Hexo支持markdown,可以生成静态网站;这或许是我从WordPress转向Hexo的理由吧。博客中有些文章已经过时了,我也想借迁移博客的机会,重新整理博客,也是重新整理自己。

写作即表达。感谢这个互联网时代,让我这样的小人物,也可以自由表达自己。

0%