針對轉(zhuǎn)錄組數(shù)據(jù)摄乒，平時(shí)分析中最常見兩組之間的比較悠反，比如不同處理或者不同突變體。面對這樣的數(shù)據(jù)用用 DESeq2 或者 edgeR 基本就差不多馍佑。如果樣本量或者不同條件很多的話可以還做WGCNA的分析斋否。但是生物體的生長發(fā)育和時(shí)間這個(gè)維度有著非常密切的關(guān)系。如果碰上了一組和時(shí)間有關(guān)系的數(shù)據(jù)可以怎么處理呢拭荤？
時(shí)序分析
所謂時(shí)序分析 (time series analysis) 在 data science 中是非常重要的一個(gè)方向茵臭。對大多數(shù)商業(yè)行為而言如果能夠通過已有不同時(shí)間數(shù)據(jù)來進(jìn)行預(yù)測就有可能大大提高自己的勝率。通常時(shí)間序列數(shù)據(jù)會包括趨勢部分和不規(guī)則部分舅世， 我們需要做的就是剔除不規(guī)則部分然后找到趨勢所在旦委，再進(jìn)行預(yù)測奇徒。在預(yù)測過程中通常可以采用移動平均法缨硝、局部加權(quán)回歸法摩钙、指數(shù)平滑法和自回歸整合移動平均等方法。
生物學(xué)時(shí)序分析
生物學(xué)的時(shí)間相關(guān)數(shù)據(jù)本身預(yù)測屬性和商業(yè)數(shù)據(jù)相比要弱很多查辩。一種是單一條件的純時(shí)間序列腺律，主要看不同基因的表達(dá)模式，根據(jù)相似的表達(dá)譜將基因歸為多個(gè)類有助于找到功能相似的基因宜肉。另一種情況是含有對照和處理的時(shí)間序列匀钧，需要再考察不同條件的差異基因。
可用的分析時(shí)間序列工具
關(guān)于時(shí)間序列轉(zhuǎn)錄組數(shù)據(jù)分析的工具谬返，近三年來有兩篇偏綜述和測評類的文章（一個(gè)人寫的）之斯。
- Dynamics in Transcriptomics: Advancements in RNA-seq Time Course and Downstream Analysis
- Comparative analysis of differential gene expression tools for RNA sequencing time course data
在這兩篇文章中還是提到了一些工具，但其中有一些用到matlab（這軟件貴扒猜痢）佑刷，有一些年久失修或者不維護(hù)或者和最新R版本不兼容，篩篩撿撿能用的且文章里認(rèn)為還不錯(cuò)的也就剩下三四個(gè)酿炸。
主要模型
- 負(fù)二項(xiàng)分布NB
來自于 DESeq 的方法瘫絮，下文中提到的 ImpluseDE2 和 MaSigPro 都使用了這種模型。
- 多項(xiàng)式回歸PR
來自于 maSigPro 方法填硕，所謂多項(xiàng)式回歸區(qū)別常見的線性回歸麦萤，會把一次特征轉(zhuǎn)換成高次特征的線性組合多項(xiàng)式，比使用直線擬合更加準(zhǔn)確扁眯。但是到底用幾次方需要具體分析壮莹，次數(shù)過高會出現(xiàn)過擬合。在能夠解釋自變量和因變量關(guān)系的前提下姻檀，次數(shù)應(yīng)該是越低越好命满，這也算是奧卡姆剃刀原則吧。
- 自回歸隱馬爾可夫模型AR-HMM
所謂自回歸是統(tǒng)計(jì)上一種處理時(shí)間序列的方法绣版，用 至
來預(yù)測本期
的表現(xiàn)并假設(shè)它們?yōu)榫€性關(guān)系胶台。簡單說就是用自己來預(yù)測自己，因?yàn)槭菑幕貧w分析中線性回歸發(fā)展而來只是用x預(yù)測x杂抽，所以叫自回歸诈唬。
差異基因檢測方法
- 對數(shù)似然比 log likelihood ratio
同樣是來自于 DESeq 的方法，下文中提到的 ImpluseDE2 和 MaSigPro 也都使用了這種方法默怨。
似然比檢驗(yàn) (likelihood ratio test,LRT) 用于比較兩個(gè)模型的擬合優(yōu)度進(jìn)而確定哪個(gè)模型與樣本數(shù)據(jù)擬合的更好讯榕。其中一個(gè)是具有一定數(shù)量項(xiàng)的完整模型，另一個(gè)是刪掉完整模型中一部分項(xiàng)的簡化模型。LRT 檢驗(yàn)中愚屁，自由度等于在簡化模型中減少的模型參數(shù)數(shù)目济竹，LR近似符合卡方分布。一個(gè)相對復(fù)雜的模型與一個(gè)簡單模型比較霎槐，如果可以顯著地適合一個(gè)特定數(shù)據(jù)集送浊，那么這個(gè)復(fù)雜模型的附加參數(shù)就能夠用在以后的數(shù)據(jù)分析中。
The LRT examines two models for the counts, a full model with a certain number of terms and a reduced model, in which some of the terms of the full model are removed. The test determines if the increased likelihood of the data using the extra terms in the full model is more than expected if those extra terms are truly zero.
為了測試多個(gè)時(shí)間點(diǎn)的任何差異丘跌，可以使用包含時(shí)間因子的設(shè)計(jì)和時(shí)間因子在簡化公式中被刪除的另一個(gè)設(shè)計(jì)袭景。對于包括對照和實(shí)驗(yàn)組的時(shí)間序列，可以使用包含條件因子闭树，時(shí)間因子和兩者相互作用的公式耸棒。在這種情況下，使用具有不包含相互作用項(xiàng)的簡化模型的似然比檢驗(yàn)將測試該條件是否在參考水平時(shí)間點(diǎn)（time 0）之后的任何時(shí)間點(diǎn)可以誘導(dǎo)基因表達(dá)的變化报辱。
- 經(jīng)驗(yàn)貝葉斯 empirical Bayesian
EBseq-HMM 采用的方法与殃，來自于 BEseq。
具體應(yīng)用
Mfuzz
這個(gè)軟件最早發(fā)表在2007年，相對老一些好在目前仍然在維護(hù)，其主要目的是給時(shí)序數(shù)據(jù)進(jìn)行基于模糊聚類算法的聚類呜魄。我們常見的聚類算法可以分為嚴(yán)格聚類(hard clustering)和模糊聚類(Fuzzy clustering )(也叫做寬松聚類 soft clustering)。嚴(yán)格聚類會將一個(gè)基因只聚到一類中爽篷，kmeans 就屬于嚴(yán)格聚類。而模糊聚類允許同一數(shù)據(jù)屬于多個(gè)不同的類慢睡，其聚類結(jié)果是一個(gè)數(shù)據(jù)對聚類中心的隸屬度逐工，0到1之間。對于分類很開的數(shù)據(jù)使用嚴(yán)格聚類是沒問題的一睁。但對于時(shí)序表達(dá)量數(shù)據(jù)來說钻弄，不同的類常常會有重疊，所以可以嘗試寬松聚類方法者吁。算法需要首先設(shè)定一些參數(shù)，若初始化參數(shù)不合適饲帅，可能影響聚類結(jié)果的正確性复凳。
在使用 Mfuzz 時(shí)首先應(yīng)該進(jìn)行數(shù)據(jù)標(biāo)準(zhǔn)化處理 ，可以使用類似于 FPKM 或者 TPM 的表達(dá)結(jié)果也可以使用 DESeq2 矯正后的結(jié)果進(jìn)行比較分析灶泵，另外不支持值為0的數(shù)據(jù)育八，所以需要加上 pseudocount 。除此之外赦邻，Mfuzz 接受的數(shù)據(jù)格式為 ExpressionSet髓棋，需要對矩陣進(jìn)行轉(zhuǎn)換。
這個(gè)包只能進(jìn)行聚類，是找不了有處理對照組的差異基因的按声。需要注意膳犹。
MaSigPro
運(yùn)行masigpro 主要有四步：
- 確定回歸模型
- 找到顯著基因
- 找到顯著差異
- 獲得顯著基因集
有兩點(diǎn)內(nèi)容需要注意：對于無對照的單一時(shí)序數(shù)據(jù)處理方法；以及處理轉(zhuǎn)錄數(shù)據(jù)時(shí)的特殊參數(shù)签则。因?yàn)檫@個(gè)包不會對數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化须床，所以應(yīng)該提前做好，使用 DESeq2 即可渐裂。
另外豺旬，在實(shí)際分析的時(shí)候可能會出現(xiàn) glm.fit: algorithm did not converge 的警告。這是由于進(jìn)行 logistic 回歸時(shí)柒凉，依照極大似然估計(jì)原則進(jìn)行迭代求解回歸系數(shù)族阅，glm 函數(shù)默認(rèn)最大迭代次數(shù)是 25，當(dāng)數(shù)據(jù)不太好時(shí) 25 次迭代可能還不收斂膝捞，一方面可以增大迭代次數(shù)坦刀。但當(dāng)增大迭代次數(shù)仍然不收斂就需要對數(shù)據(jù)進(jìn)行異常值檢驗(yàn)等進(jìn)一步處理。通常把一些表達(dá)量極低或者極高的基因刪除掉绑警，這個(gè)問題就可以解決求泰。
ImpulseDE2
ImpulseDE2 是最近才出來的一個(gè)R包，在前面提到的綜述評測文章中認(rèn)為這個(gè)包找時(shí)序數(shù)據(jù)中的差異基因效果最好计盒，它可以用來解決兩類問題渴频。
- Case-only differential expression analysis tests, whether the expression level of a gene changes over time.
- Case-control differential expression analysis tests, whether the expression trajectory of a gene over time differs between samples from a case and samples from a control condition.
這個(gè)包中，有一個(gè)plotHeatmap 函數(shù)北启，可以借助 ComplexHeatmap 對數(shù)據(jù)整體進(jìn)行熱圖的繪制同時(shí)提取不同類的基因卜朗，也可以使用plotGenes看某一個(gè)基因的表達(dá)情況。
在展示的熱圖中會出現(xiàn)四部分咕村，包括 transient and transition trajectorie场钉，其中每一種 tarjectorie 又包括 up 和 down 兩類。所謂的 transient 可以理解為時(shí)序數(shù)據(jù)在中間某一個(gè)時(shí)間點(diǎn)存在up 或者 down peak懈涛，即在某一個(gè)時(shí)間點(diǎn)存在表達(dá)的最大或者最小值逛万；而所謂的 transient 可以理解為一個(gè)持續(xù)的變化，比如持續(xù)的升高或者持續(xù)的降低批钠。
EBSeq-HMM
EBSeq-HMM 是基于 EBSeq 二次開發(fā)的工具宇植，主要用于分析時(shí)序數(shù)據(jù)。在計(jì)算的時(shí)候首先基于負(fù)二項(xiàng)分布對參數(shù)進(jìn)行估計(jì)埋心，然后利用自回歸隱馬模型將基因的表達(dá)進(jìn)行分類指郁。比較神奇的是，最終給到的結(jié)果會標(biāo)示為 Up-Up-Down-Down-Down 之類的若干 path拷呆，然后你可以選出你感興趣的 path 進(jìn)行后續(xù)分析闲坎。
后續(xù)感受
因?yàn)槟壳白龅臄?shù)據(jù)是沒有對照的單一時(shí)間序列數(shù)據(jù)疫粥，所以還不能體會哪一個(gè)找出的差異基因更準(zhǔn)確些。但是如果只是想把所有的基因根據(jù)不同的時(shí)間點(diǎn)分為若干表達(dá) pattern腰懂，似乎結(jié)合 Mfuzz 和 ImpulseDE2 就可以了梗逮。
當(dāng)然，涉及到聚類悯恍，尤其是非監(jiān)督聚類的時(shí)候通常主觀因素還是較強(qiáng)库糠，如果能對關(guān)鍵基因或者數(shù)據(jù)有一個(gè)大致的估計(jì)預(yù)判操作起來會相對輕松些，如果沒有涮毫，可能就需要結(jié)合不同類的生物學(xué)意義等角度來找合適的聚類數(shù)目了瞬欧。
參考資料
http://a-little-book-of-r-for-time-series.readthedocs.io/en/latest/src/timeseries.html
http://a-little-book-of-r-for-biomedical-statistics.readthedocs.io/en/latest/
https://laranikalranalytics.blogspot.com/2018/07/time-series-analysis-with-documentation.html
https://www.displayr.com/smoothing-time-series-data/?utm_medium=Feed&utm_source=Syndication
https://www.analyticsvidhya.com/blog/2015/12/complete-tutorial-time-series-modeling/
