2022.3.20 初版
2022.3.28 加入代碼
眾所周知，差異基因的分析方法有很多種：①普通的差異基因計算方法（limma、edgeR、DESeq2等R包）（缺點：針對不同數(shù)據(jù)集和不同樣本來說，批次效應較大，有時候做出來的預測模型不能用于其他的數(shù)據(jù)，魯棒性較低），②Robust rank aggregation(RRA)排秩法，通過對幾個數(shù)據(jù)集差異基因進行排序，篩選出普遍適用的差異基因（優(yōu)點：能高效去除批次效應，魯棒性較高，結果、模型等普遍適用于其他的數(shù)據(jù)集；缺點：同時分析多個數(shù)據(jù)集，單個數(shù)據(jù)集無法使用）（這個方法網(wǎng)上很多教程也有，只是不是很詳細，你們可以去讀讀然后思考思考，我看情況今后也會講講）③基于基因秩序表達秩序關系識別穩(wěn)定逆轉基因對（就是接下來我要講的這個）（優(yōu)缺點待會你們就知道啦）
誤入BioInfor的大黃鴨 --一個喜歡把教程寫著寫著寫成科普的本科臨床醫(yī)學生
在基因測序的統(tǒng)計學分析中，我們難免會遇到實驗批次、時間環(huán)境因素、芯片平臺等使測序結果帶來的差異性，有時候我們得出的結果、預后模型、診斷模型，也僅僅只在我們的這個數(shù)據(jù)集內(nèi)可用，利用結果去測試其他數(shù)據(jù)就會發(fā)現(xiàn)模型無法使用或者效果不佳。此時就有很多大神創(chuàng)造出許多去除批次效應的方法。
今天我們講的是基于基因秩序表達秩序關系識別穩(wěn)定逆轉基因對的方法，這是我偶然看到的一篇電子科技大學的碩士論文上的一個方法[1]。當時也是第一次聽說這個方法，其實我腦子里一直在想一種去除批次效應的方法，一直有這種類似于基因編碼的形式，但是形不成系統(tǒng)，最近看到這篇文章后可以說是“不識廬山真面目，只緣身在此山中”了。然后我在知網(wǎng)和PubMed上就查了這個方法，感覺之前還挺多人用的，很感謝這個碩士論文的作者，發(fā)現(xiàn)了新大陸。找了一天的百度，一直沒找到這個方法的教程，按照以前在文章中看到挺好用的方法我都會馬上百度搜索怎么使用，現(xiàn)在這種情況太難了，因此打算自己按圖索驥，按照這個方法自己花三天時間寫出了個R代碼。
穩(wěn)定逆轉基因對（relative expression orderings，REOs），即一種基因的秩序存在的模式，假設基因i（Gi）和基因j（Gj），若在95%（或90%[a]）的正常樣本中Gi>Gj（或Gi<Gj），我們稱之為穩(wěn)定基因對。而在95%（或90%）腫瘤樣本中情況相反，則我們定義為穩(wěn)定逆轉基因對（REOs），穩(wěn)定逆轉基因對直接用“是”或“否”來定義是否具有差異，故該方法計算結果無差異倍數(shù)，但是對每個樣本進行了量化，按照每個樣本的情況來定義的，不像普通差異分析方法是按照平均表達量來計算，因此批次效應會低很多。因此此差異分析方法僅適用于各模型的構建及特征基因的篩選。
[a]其實這個數(shù)是自己隨便主觀定義的數(shù)，只要篩選出的穩(wěn)定逆轉基因對的最后的模型結果夠好就行了，80%也可以，我推薦使用90%或95%，因為畢竟有參考文獻支持，審稿人問你這個數(shù)怎么來的，你也可以有底氣的回答。
例：
| ID | 樣本A_con | 樣本B_con | 樣本C_con | 樣本D_con | 樣本E_con | 樣本F_treat | 樣本G_treat | 樣本H_treat | 樣本I_treat | 樣本J_treat |
|---|---|---|---|---|---|---|---|---|---|---|
| Gene A | 0.85 | 0.82 | 0.73 | 0.89 | 0.83 | 0.24 | 0.28 | 0.46 | 0.25 | 0.33 |
| Gene B | 0.47 | 0.53 | 0.26 | 0.37 | 0.43 | 0.67 | 0.73 | 0.83 | 0.66 | 0.76 |
若Gene A>Gene B，我們用1表示，反之用0表示，則：
| ID | 樣本A_con | 樣本B_con | 樣本C_con | 樣本D_con | 樣本E_con | 樣本F_treat | 樣本G_treat | 樣本H_treat | 樣本I_treat | 樣本J_treat |
|---|---|---|---|---|---|---|---|---|---|---|
| Gene A / Gene B | 1 | 1 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 0 |
在95%（或90%）的樣本中，正常組Gene A都大于Gene B，腫瘤組Gene A都小于Gene B，因此我們把Gene A / Gene B定義為穩(wěn)定逆轉基因對。
應用：
穩(wěn)定逆轉基因對的應用主要在于特征基因的篩選，如輸入最小冗余最大關聯(lián)算法（mRMR）、隨機森林算法（RF）、LASSO、支持向量機（SVM）等，還有人工神經(jīng)網(wǎng)絡模型的構建。可作為分類器的特征篩選使用，對亞型分型、個性化治療、早期診斷等提供幫助。
Rstudio操作：
我們先準備好一個正常樣本的矩陣，一個腫瘤樣本的矩陣，可以是基因，外泌體，lncRNA，miRNA，circRNA，可以是轉錄組，蛋白組，影像組，代謝組，糖組，免疫組等（以circRNA為例）（Tumor.txt和Normal.txt）
文件格式是這樣的：
代碼如下：
setwd("D:\\Users\\123\\Desktop\\StabRvsPair") #設定工作目錄
Normal = read.table("Normal.txt",header=T,sep="\t",check.names=F,row.names=1) #讀取Normal矩陣
Tumor = read.table("Tumor.txt",header=T,sep="\t",check.names=F,row.names=1) #讀取Tumor矩陣
#加載包
library(StabRvsPair)
library(stringr)
PairTime=testRunTime(Normal,Tumor,0.9,0.1) #測試電腦運行速度
GetPair=FindStabRvsPair(Normal,Tumor,0.9,0.1) #獲取穩(wěn)定逆轉基因對
GenePair=CalStabRvsPair(GetPair,Normal,Tumor) #獲取穩(wěn)定逆轉基因對矩陣
StabRvsPair包是我自己寫的R包，不歸CRAN管理，也不歸BioConductor管理，我還沒放上GitHub。這是我人生中辛辛苦苦寫的第一個R包，真的挺不簡單的，唉，做的不好，還是得繼續(xù)努力！如果想要的朋友可以聯(lián)系我獲取哦！
注意事項：
一個10000多的基因的表達矩陣，至少能篩選出上億基因對，并進行基因對的對比，需要花很多時間，很多算力（淺算一下可能得花一個月的時間）。因此建議內(nèi)存128G的服務器參與計算，或者進行多線程并行計算。如果是個人電腦，建議拋棄R使用Python語言，或者選一個基因相對較少（建議2000-5000個基因的表達矩陣進行計算，可能需要在48G內(nèi)存的服務器上掛一個天，淘寶價格大約是100塊錢左右）（實測1000個基因在48G運行內(nèi)存服務器上需要跑3小時出結果，當然是我算法效率的問題，優(yōu)化算法或使用多線程并行會快幾倍）（1000個基因的表達矩陣才能跑出3個穩(wěn)定逆轉基因對，還是用80%的條件）（后來發(fā)現(xiàn)這個速度也不完全靠內(nèi)存，CPU的單核主頻非常重要，建議還是用好點的電腦跑）
參考文獻：
[1]張子梅. 利用機器學習方法識別肝癌早期診斷標志[D].電子科技大學,2021.DOI:10.27005/d.cnki.gdzku.2021.002251.
由于StabRvsPair包是我辛辛苦苦自己寫的，來之不易，想用的朋友在微公號發(fā)送“基因對”。
本教程就先講到這啦，后續(xù)隨后更新，歡迎大家關注支持～大家關注一下我公號：誤入BioInfor的大黃鴨
