2022.3.20 初版
2022.3.28 加入代码
众所周知，差异基因的分析方法有很多种：①普通的差异基因计算方法（limma、edgeR、DESeq2等R包）（缺点：针对不同数据集和不同样本来说，批次效应较大，有时候做出来的预测模型不能用于其他的数据，鲁棒性较低），②Robust rank aggregation(RRA)排秩法，通过对几个数据集差异基因进行排序，筛选出普遍适用的差异基因（优点：能高效去除批次效应，鲁棒性较高，结果、模型等普遍适用于其他的数据集；缺点：同时分析多个数据集，单个数据集无法使用）（这个方法网上很多教程也有，只是不是很详细，你们可以去读读然后思考思考，我看情况今后也会讲讲）③基于基因秩序表达秩序关系识别稳定逆转基因对（就是接下来我要讲的这个）（优缺点待会你们就知道啦）
误入BioInfor的大黄鸭 --一个喜欢把教程写着写着写成科普的本科临床医学生
在基因测序的统计学分析中，我们难免会遇到实验批次、时间环境因素、芯片平台等使测序结果带来的差异性，有时候我们得出的结果、预后模型、诊断模型，也仅仅只在我们的这个数据集内可用，利用结果去测试其他数据就会发现模型无法使用或者效果不佳。此时就有很多大神创造出许多去除批次效应的方法。
今天我们讲的是基于基因秩序表达秩序关系识别稳定逆转基因对的方法，这是我偶然看到的一篇电子科技大学的硕士论文上的一个方法[1]。当时也是第一次听说这个方法，其实我脑子里一直在想一种去除批次效应的方法，一直有这种类似于基因编码的形式，但是形不成系统，最近看到这篇文章后可以说是“不识庐山真面目，只缘身在此山中”了。然后我在知网和PubMed上就查了这个方法，感觉之前还挺多人用的，很感谢这个硕士论文的作者，发现了新大陆。找了一天的百度，一直没找到这个方法的教程，按照以前在文章中看到挺好用的方法我都会马上百度搜索怎么使用，现在这种情况太难了，因此打算自己按图索骥，按照这个方法自己花三天时间写出了个R代码。
稳定逆转基因对（relative expression orderings，REOs），即一种基因的秩序存在的模式，假设基因i（Gi）和基因j（Gj），若在95%（或90%[a]）的正常样本中Gi>Gj（或Gi<Gj），我们称之为稳定基因对。而在95%（或90%）肿瘤样本中情况相反，则我们定义为稳定逆转基因对（REOs），稳定逆转基因对直接用“是”或“否”来定义是否具有差异，故该方法计算结果无差异倍数，但是对每个样本进行了量化，按照每个样本的情况来定义的，不像普通差异分析方法是按照平均表达量来计算，因此批次效应会低很多。因此此差异分析方法仅适用于各模型的构建及特征基因的筛选。
[a]其实这个数是自己随便主观定义的数，只要筛选出的稳定逆转基因对的最后的模型结果够好就行了，80%也可以，我推荐使用90%或95%，因为毕竟有参考文献支持，审稿人问你这个数怎么来的，你也可以有底气的回答。
例：
| 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定义为稳定逆转基因对。
应用：
稳定逆转基因对的应用主要在于特征基因的筛选，如输入最小冗余最大关联算法（mRMR）、随机森林算法（RF）、LASSO、支持向量机（SVM）等，还有人工神经网络模型的构建。可作为分类器的特征筛选使用，对亚型分型、个性化治疗、早期诊断等提供帮助。
Rstudio操作：
我们先准备好一个正常样本的矩阵，一个肿瘤样本的矩阵，可以是基因，外泌体，lncRNA，miRNA，circRNA，可以是转录组，蛋白组，影像组，代谢组，糖组，免疫组等（以circRNA为例）（Tumor.txt和Normal.txt）
QQ截图
文件格式是这样的：
QQ截图
代码如下：
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) #获取稳定逆转基因对
GenePair=CalStabRvsPair(GetPair,Normal,Tumor) #获取稳定逆转基因对矩阵
StabRvsPair包是我自己写的R包，不归CRAN管理，也不归BioConductor管理，我还没放上GitHub。这是我人生中辛辛苦苦写的第一个R包，真的挺不简单的，唉，做的不好，还是得继续努力！如果想要的朋友可以联系我获取哦！
注意事项：
一个10000多的基因的表达矩阵，至少能筛选出上亿基因对，并进行基因对的对比，需要花很多时间，很多算力（浅算一下可能得花一个月的时间）。因此建议内存128G的服务器参与计算，或者进行多线程并行计算。如果是个人电脑，建议抛弃R使用Python语言，或者选一个基因相对较少（建议2000-5000个基因的表达矩阵进行计算，可能需要在48G内存的服务器上挂一个天，淘宝价格大约是100块钱左右）（实测1000个基因在48G运行内存服务器上需要跑3小时出结果，当然是我算法效率的问题，优化算法或使用多线程并行会快几倍）（1000个基因的表达矩阵才能跑出3个稳定逆转基因对，还是用80%的条件）（后来发现这个速度也不完全靠内存，CPU的单核主频非常重要，建议还是用好点的电脑跑）
参考文献：
[1]张子梅. 利用机器学习方法识别肝癌早期诊断标志[D].电子科技大学,2021.DOI:10.27005/d.cnki.gdzku.2021.002251.
由于StabRvsPair包是我辛辛苦苦自己写的，来之不易，想用的朋友在微公号发送“基因对”。
本教程就先讲到这啦，后续随后更新，欢迎大家关注支持～大家关注一下我公号：误入BioInfor的大黄鸭
网友评论