Fission：筛选多时间点测序数据的差异表达基因

  研究动态生物学过程的时候，例如发育过程或药物反应，往往需要设计多个时间点的实验（Time Course experiment，即TC）。一般的RNA-seq转录组分析是两个样本与两个样本比较，或两组样本与两组样本比较，显然这样的分析方法，对TC的RNA-seq数据是不恰当的，最大的缺陷就是没有考虑某个时间点的基因表达与前或后一个时间点之间的关系，当然对数据标准化、噪音校正、计算差异的时候也没考虑TC实验的随机变动、整体时间的变化趋势、时间位移的影响。比如，某种药物处理之后，细胞群体的代谢会变慢，导致基因表达模式的改变，对于这种延迟效应不应该采取两两比较的方式，而应该通盘考虑各个时间点的变化。

  R包fission就是专门针对TC实验的分析工具，可以在bioconductor网站下载。本文重点介绍它的使用方法。

  输入文件有两个，一个是“sample.txt”，记录了样本信息，如下图：

  另一个文件是“count.txt”，记录了每个基因在每个样本中的表达值，即原始的reads count，一行代表一个基因，一列代表一个样本。如下图：

library("fission")

读入表达值：

countData

读入样本信息：

colData

数据格式转化为因子：

colData$minute

表明计算差异的依据（这里是strain和minute，即不同菌株在不同时间点对药物的反应）：

dds

计算差异表达基因（方法为likelihood ratio test，即LRT）：

#Genes at one or more time points aftertime 0 showed a strain-specific effect.

ddsTC

resTC

添加基因的名字：

resTC$symbol

head(resTC[order(resTC$padj),],4)

画出adjusted p value最小的基因的表达模式（不同菌株在不同时间的表达值）

data

ggplot(data, aes(x=minute, y=count, color=strain, group=strain)) + geom_point() + stat_smooth(se = FALSE, method = "loess") + scale_y_log10()

对差异基因聚类

betas

差异最显著（按照adjusted p value排序）的前20个基因

topGenes

对log2（Foldchange）画图：

mat thr]

pheatmap(mat, breaks=seq(from=-thr, to=thr,length=101),cluster_col=FALSE)

##