2021.11.12 初版
2021.11.16 更新单因素COX输出结果，多了半边括号...
本系列主要用于通过使用训练集（train）、内部验证集（test1组）、外部验证集（test2）联合对感兴趣的基因进行Cox回归模型的构建，并用于预测患者的临床预后情况，通过生存曲线、ROC曲线、生存状态图对模型效能进行可视化
胡乱调侃，与正题无关，赶时间的友友直接刷掉哈~
误入BioInfor的大黄鸭 --一个喜欢把教程写着写着写成科普的本科临床医学生
大黄鸭自从误入了生信的大坑之后，也走了很多弯路（毕竟自己也只是个临床医学的苦逼大三生），花了很多钱不说，电脑游戏软件图标逐渐被Rstudio等软件给取代。。。甚至被舍友给逼问为啥不直接报计算机专业？？？诚然，，自学生信花的时间也太多了吧。。从R说起，卑微的我随后在淘宝上开了个小店（杂货铺？），以别人R语言报错改代码为生，就以这点微薄收入度日，主要不是为了挣钱，还是为了认识更多师兄师姐，而且帮别人改代码真的学到了很多很多，R语言水平也提升了很多很多（离谱的是还有人为了8块钱改个代码跟我讲过价。。。）之后就在想为啥不把这些经验分享给友友们，好让友友们学生信少走点弯路（嘻嘻有一说一，其实还是想让大家关注一下我，以后我修改代码的话也会实时在里面更新教程，有学员亲测看了我改代码过程一个月之后，他R语言啥报错再也没怕过，还说想把我微信给删了）
（本文第一次试水，多有不足，欢迎大家批评指正）
首先我们先来讲包和输入文件的准备
目录
- 0.脑子的准备
- 1.包的准备
- 2.输入文件的准备
- 3.分组
- 4.单因素Cox分析
- 5.LASSO回归
- 6.多因素Cox回归（构建Cox模型）
- 7.模型验证
- 8.输出结果文件
- 9.可视化
Cox回归模型介绍：参见百度百科
算了还是简单讲讲吧。。又是废话一通，请各位知道Cox回归模型是啥东西的大神自行略过。
在我们日常观看的电视电影里，经常看到此景：“医生，我父亲还能活几年啊？”，“10-15年”
为啥是10-15年呢？是医生随便说的，没有任何依据吗？我乱报个数，说20年行不行？
不是的。
这时Cox回归模型就开始刷它的存在感了
Cox模型是由英国统计学家D.R.Cox(1972)年提出的一种半参数回归模型，它以生存结局和生存时间为因变量，可同时分析众多因素对生存期的影响。在这里因素可以是临床因素，比如说肿瘤各种分期，年龄，患者体征，甚至可以是患者的生活习惯等，也可以是各种基因的表达量，可以纳入你所感兴趣的基因进行模型的构建，它是一个系数×表达量的公式，通过计算得出的风险值来预测患者的预后情况。本例以基因表达作为例子。
本例基本思路：
机器学习一般可将数据分为训练集（train）、内部验证集（test1组）、外部验证集（test2），在本教程中我们在TCGA获取的数据进行随机分组，随机分为训练集和内部验证集，在ICGC数据库中我们获取数据对模型进行外部验证。训练集负责训练模型，而内部验证集我们用于对Cox模型的效能进行检验，随机分组则体现出数据的随机性及普遍性，而为了进一步说明模型的可用性，我们使用其他数据库中的数据，使用模型检验效果，数据的来源可以是其他数据库中的基因表达和临床数据，也可以使用我们自己实验中获得的患者基因表达和临床数据。
单因素Cox分析主要是分析其因素单独有没有对生存有影响，跟Kaplan-Meier生存曲线差不多，算法不一样
1.包的准备
install.packages("glmnet")
install.packages("caret")
install.packages("survival")
install.packages("survivalROC")
install.packages("survminer")
直接跑代码就完事了，后续我更新一个方便批量安装包的教程哦~（疯狂暗示-关注一波大黄鸭嘻嘻）
2.输入文件的准备
我们需要的数据文件是TCGA数据库中FPKM的表达数据文件，以及临床数据文件（包含患者生存状态及生存时间的），将其整理到一个文件中
表格样式为：
| ID | time | status | Gene1 | Gene2 | ... | GeneN |
|---|---|---|---|---|---|---|
| TCGA样本名 | 生存时间（天） | 生存状态（0/1 活着/死亡） | Gene1exp | Gene2exp | ... | GeneNexp |
样本名用TCGA-XX-XXXX的格式，可以用Excel的（mid函数提取前12位）免得到时候再代码里删，因为后面的部分为样本的信息，前面的才是患者信息，如果遇到重复行名（rowname），需要在Excel上寻找两个一样的，删掉其中一个，生存时间一定要用年啊！以天为单位的要除以365，这些工作直接用Excel就完事了。可以用VLOOKUP函数，也可以手动整理，工作量一般般吧，或者等我后续更新的教程也行，关注大黄鸭我就完事了哈哈哈~
当然，train集和test1集都是在同一个表达矩阵中进行随机分组的，而test2集需要通过外部数据库（ICGC或GEO等）或者自己实验的数据（由于本例纳入计算的数据为TCGA和ICGC数据库的Fpkm值，其他的标准可能与模型不适用，建议还是以同一标准处理的RNA-Seq量作为外部验证组的验证（比如说用edgeR标准化的TCGA Counts矩阵就不能用Fpkm标准化矩阵进行外部验证））整理出另外一个矩阵
接下来就是读取数据啦，要记得更改工作路径哦~（setwd()走一波）
我们这文件保存为exp_time1.txt（train集和test1集），exp_time2.txt（test2集）（timeexp不好看。。两个e给整重叠了），使用代码读取：
exptime=read.table("expTime1.txt" , header=T , row.names=1 , sep="\t" , check.names=F)
test2=read.table("expTime2.txt" , header=T , row.names=1 , sep="\t" , check.names=F)
header，row.names表示设置表头和行名，sep的这个\t表示识别空格为分隔符。。毕竟俺们保存为txt的时候，每个词之间是变成空格的嘛，check.names这个如果不设置成F，到时候列名前面多个“X.”的话别来找我
3.分组
咱前面说了，这个TCGA得随机分组，也就是分成训练集和内部验证集，这一块有人55分的，也有人64分的，73分的都有，咱这用的是73分，就是训练集7份，内部验证集3份，当然怎么分还是得看您的心情哈哈哈
随机分组，结果不一定好啊，所以这里咱得用循环，结果不好的话就重新再随机分组，一直到出现好的结果就输出。所以，这里得弄个循环，凭咱电脑配置，咱这设置为最多循环800次，也够跑半天了，当然您电脑差点的话，稍微调低一点也行，咱也不敢反对，否则怕您到时候电脑出问题啥的，或者跑了半天也没出结果，来找我麻烦也不太好。如果您电脑是几颗锐龙3990x+128G多少频率的内存，调太高也不好，毕竟分组几万次才出个好结果，这也太侥幸了吧，太过违背随机原则了。
for(i in 1:800){
...
}
构建一个这样的循环，本文下面的所有代码都放到里面去就好了。
咱这用了createDataPartition这个函数，当然您用随机数生成的seed也行，只是这个函数专门为机器学习分组的，真的很香啊，他可以从我们样本中随机抽出train和test，死亡和存活的样本都有，而且比例不是很过分，而用随机数的话，假如倒霉，train组大部分抽到死亡的，test组大部分抽到生存的，就不好了。所以建议你们还是用这个函数，真的很方便机器学习这块。。当然随机数生成seed也行（毕竟那么倒霉的几率很小），这块还得友友们亲自去百度，后续如果我觉得价值大可能我也会更新。
grouplist<-createDataPartition(y=exptime[,2] , p=0.7 , list=F)
train<-exptime[grouplist,]
test1<-exptime[-grouplist,]
代码解析：y=exptime[,2]是为了以第二排作为判断分类的标准，生存0，死亡1，以免出现train组大部分抽到死亡的，test组大部分抽到生存的
p=0.7当然就是抽70%了
list=F是不输出为list这个数据类型（有点难用），直接输出为矩阵
下面就是生成train组，test1组矩阵了
本教程就先讲到这啦，后续随后更新，欢迎大家关注支持～大家关注一下我公众号：误入BioInfor的大黄鸭，回复“TCGACox”获取完整版的代码
网友评论