一直不清楚MetaStat的脚本输入到底应该是Read数还是相对丰度,大概过了一下R语言版本的MetaStat代码的脚本,留此记录。
代码解析
R语言版本的MetaStat被封装到了一个R函数内部,针对不同的数据输入,MetaStat会采用不同的处理策略
两个比较组都是单样本无重复
此时MetatStat使用Fisher精确检验或卡方检验(每个样本内tag或reads数都大于20)
涉及到Fisher或卡方检验的输入都是列联表,肯定都是整数因此理想的输入的表格内都应该是整数
有重复
MetaStat对于有重复的样本(或者说是上述情况以外)的比较组的统计分为两步
- 对于输入的每一行,都会通过T检验计算其t统计量(这里使用的是每个样本的相对丰度),然后通过t统计量利用置换检验算出p值
- 重新遍历每一行,对于在两个比较组内平均tag数小于1的行,重新利用Fisher精确检验计算p值替换前一步得到的该物种的p值
感想
综上所述,应该使用count/reads值作为其函数的输入
吐槽
- 这个R代码内各种for循环~ 效率真是不高~
- 有空看一下置换检验,这里置换检验也是限制代码执行时间的坑~
好像mothur里面有C++的版本,有时间试一下