之前讲过临床模型预测的专栏，但那只是基础版本，下面我们以自噬相关基因为例子，模仿一篇五分文章，将图和代码复现出来，学会本专栏课程，可以具备发一篇五分左右文章的水平：
本专栏目录如下：
Figure 1：差异表达基因及预后基因筛选(图片仅供参考)
Figure 2. 生存分析，箱线图表达改变分析（图片仅供参考）
Figure 3. 基因富集分析（图片仅供参考）
Figure 4.构建临床预测模型（图片仅供参考）
Figure 5.训练集训练模型（图片仅供参考）
Figure 6.测试集测试模型（图片仅供参考）
Figure 7.外部数据集验证模型（图片仅供参考）
Figure 8.生存曲线鲁棒性分析（图片仅供参考）
FIgure 9.列线图构建，ROC分析，DCA分析（图片仅供参考）
Figure 10.机制及肿瘤免疫浸润（图片仅供参考）
########################### 下面是我们要做的图
Figure 1：差异表达基因及预后基因筛选(图片仅供参考)
################### 首先是下载数据 ##################
本专栏用到的数据是TCGA数据库里面的KIRC数据，我已经整理成了LCPM格式，并放在我的资源中，可以直接下载使用：
TCGA-KIRC-mRNA表达数据——肾透明细胞癌表达及临床数据集整理_TCGA-数据挖掘文档类资源-CSDN下载TCGA-KIRC数据集已经整理成LCPM格式，临床数据已经汇总整理。LCPM格式即log2(CPTCGA更多下载资源、学习资料请访问CSDN下载频道./download/weixin_46500027/85440520?spm=1001..3001.5503
该数据里面包含了表达数据和临床数据，如果大家想用TPM格式的数据，可以到官网下载后自己转化格式，转化格式的代码在我 的另一个专栏有相关介绍：
/weixin_46500027/category_11845300.html?spm=1001..3001.5482/weixin_46500027/category_11845300.html?spm=1001..3001.5482下面我们开始下载自噬相关基因，这是我们研究的某个方向，大家可以自己去找别的研究方向，模仿本专栏提供的课程和代码。
自噬相关基因下载如下：
自噬相关基因有一个数据库，在我们之前的文章中有讲过：
注意：由于自噬数据库最近打不开了，所以我将自噬相关基因上传到百度网盘上，供大家练习：
链接：/s/1dgUVqcoVoQvmLvQtKukSQA 提取码：z9m2 --来自百度网盘超级会员V5的分享
简单来说呢，就是把这些基因复制，粘贴到Excel里面。
全选，复制：
粘贴完了以后，我们需要对空格进行一下处理，因为自噬数据库呢，没有下载的地方我们只能这样手动粘贴，然后手动删减。
总之，处理好了以后，就剩下三列，我们最主要的就是要Symbol这一列。
注意：由于自噬数据库最近打不开了，所以我将自噬相关基因上传到百度网盘上，供大家练习：
链接：/s/1dgUVqcoVoQvmLvQtKukSQA 提取码：z9m2 --来自百度网盘超级会员V5的分享
下面我们得准备KIRC的测序数据和临床数据，大家可以取UCSC xena下载。不过下载的数据是FPKM格式的，需要传承TPM或者LCPM格式。
因为FPKM格式做的文章是容易受到质疑的。
在我的资源里面已经对TCGA的数据库进行了格式转换，id转换，大家也可以直接进入我的资源去下载：TCGA-KIRC-mRNA表达数据——肾透明细胞癌表达及临床数据集整理_黑色素瘤TCGA-数据挖掘文档类资源-CSDN下载TCGA-KIRC数据集已经整理成LCPM格式，临床数据已经汇总整理。LCPM格式即log2(CP黑色素瘤TCGA更多下载资源、学习资料请访问CSDN下载频道./download/weixin_46500027/84997590?spm=1001..3001.5503
下载完了以后解压缩，和前面的自噬相关基因放在一起。
下面我们挨个看看这些文件里面的内容：
第一个KIRC_clinicalMatrix:
这个临床文件是整理好了的，包含了所有KIRC的临床信息。
下面的KIRC_LCPM是RNA-sequence数据：
这个数据是id转换完了，也从FPKM转成了LCPM格式，包含肿瘤和正常组织的6万多个基因的测序数据，是可以直接拿来用的。
当然有些研究者也还在沿用TPM格式，后面我也会陆续上传到我的资源里面，或者大家自己转一转格式。
这里使用LCPM格式的原因，是因为我之前发的一篇六分的文章，审稿人说FPKM和TPM格式有点过时了，建议我转成LCPM格式，所以我才跟大家介绍这种格式。我也特地去查了查文献，确实LCPM格式的稳定性比其他两种要强一些。这些都是有文献支持的。
这是RNA-seq数据。
下面是自噬相关基因的数据：
我们将这些数据放在同一文件夹里面，然后开始进行差异表达分析，和基因匹配：
下面进行差异表达分析和基因匹配，代码如下：
############################# 差异表达基因筛选 #############################setwd("C:\\Users\\ASUS\\Desktop\\自噬相关基因临床模型预测\\1数据准备")dir()data <- read.csv("KIRC_lcpm.csv",sep=",",header=T) data[1:4,1:4]names=as.data.frame(data$X)head(names)names(names) <- "Name"grep <- grep("^TCGA[.]([a-zA-Z0-9]{2})[.]([a-zA-Z0-9]{4})[.]([0][0-9][A-Z])",colnames(data))length(grep)greptumor <- data[,grep]tumor[1:4,1:4]grep1 <- grep("^TCGA[.]([a-zA-Z0-9]{2})[.]([a-zA-Z0-9]{4})[.]([1][0-9][A-Z])",colnames(data))length(grep1)grep1normal <- data[,grep1]normal[1:4,1:4]data <- cbind(names,tumor,normal)data[1:5,1:5]library(limma)rt=as.matrix(data)rownames(rt) <- rt[,1]exp=rt[,2:ncol(rt)]dimnames=list(rownames(exp),colnames(exp))data33=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)library(limma)dim(data33)data33=avereps(data33)dim(data33)head(data33)dim(data33)dim(data)data33 <- normalizeBetweenArrays(data33)data <- as.data.frame(data33)class(data)length(grep)length(grep1)tumorNum <- length(grep)NormalNum <- length(grep1)modType=c(rep("tumor",tumorNum),rep("normal",NormalNum))design <- model.matrix(~0+factor(modType))designcolnames(design) <- c("con","treat")fit <- lmFit(data,design)cont.matrix<-makeContrasts(treat-con,levels=design)fit2 <- contrasts.fit(fit, cont.matrix)fit2 <- eBayes(fit2)allDiff=topTable(fit2,adjust='fdr',number=200000)write.csv(allDiff,file="TCGA_KIRC_diff_expression.csv")
然后匹配自噬相关基因：
################### 自噬相关差异表达基因匹配 #############gene <- read.csv("自噬相关基因临床.csv",header = T,sep = ",")match(gene$Symbol,rownames(allDiff))# > match(gene$Symbol,rownames(allDiff))[1] 8136 3854 10824 5309 9322 1715 2549 14382 340 4688 1416 6940 5219 11515 8966 11522 10366 7480 2789[20] 14321 13938 4424 11773 734 6949 12209 527 7182 12563 4102 746 907 9025 13908 1335 2889 NA NA[39] 8083 1714 13271 8276 11229 561 10819 359 3485 563 1725 12713 4563 2366 4311 8199 10 7433 4296[58] 10132 NA 13729 13055 NA 4770 28 9210 13614 13141 784 13230 10601 7180 2781 8790 11978 1738 2664[77] 11985 13234 12769 1637 9401 450 6744 485 11159 NA 1350 1635 2698 7035 3141 13796 8719 13704 243[96] 3625 1754 12560 2656 11668 9233 NA 9265 10290 7236 1551 7743 9645 1838 541 NA 9039 8789 NA[115] NA 9762 11765 11788 11025 3913 8719 13704 243 3625 1754 12560 2656 11668 9233 NA 7856 NA NA[134] 3620 11954 732 13199 7531 13158 NA 5801 12203 6523 5684 3120 9479 13974 11738 4073 130 1539 9292[153] 9172 3954 3806 3545 9488 NA 1464 4375 10457 NA 1804 1322 4136 9510 2948 10314 11587 2566 12104[172] 4729 3858 6207 4622 1818 6319 1376 12499 NA 5942 6835 1543 13483 2195 6486 8469 986 7790 11363[191] 4141 228 11684 3924 11882 10395 8163 6086 7128 14181 2699 4790 5948 5866 7562 8317 5603 12222 NA[210] NA 8661 3715 12295 NA NA 12853 8354 6108 7075 11324 9693 13685 12699 9675 NA 215 8353 9967[229] NA 13961 11118 4629
此处可以看到匹配的结果中有很多NA，说明没有匹配上，原因可能是两个数据库的同一个基因用了不同的基因名字，因此我们需要去genecard里面将这些NA的基因名字换成相同的基因名。
例如：我们发现我自噬基因数据的第115个基因是NA，我们查看一下第115个基因是什么：
gene[115,]#XName Symbol# 115 345611 immunity-related GTPase family, M IRGM
第115个基因是IRGM，然后我们去genecard里面查找一下它有没有别的名字：
我们将IRGM替换成IRGM1，然后重新匹配：
再次运行代码：
gene <- read.csv("自噬相关基因临床.csv",header = T,sep = ",")match(gene$Symbol,rownames(allDiff))gene[115,]# > match(gene$Symbol,rownames(allDiff))[1] 18795 7434 29495 10189 21064 3293 5170 36709 799 8714 2426 15466 11069 34093 23523 29276 30793 17671 5805[20] 52608 56472 5429 34466 1719 14982 40942 1222 14808 36598 8650 1467 1217 20852 35516 2690 5503 NA NA[39] 16473 3164 39131 18600 24767 1396 29198 1064 5710 1273 3162 39138 7651 4764 8276 16969 85 15273 8754[58] 22950 NA 44885 49259 NA 9276 198 21472 33642 39701 2095 46241 27259 14728 5254 18546 27863 3274 5123[77] 36294 41356 37320 3376 22490 1167 12977 1174 25311 NA 3089 3142 8078 15054 6293 48066 18701 47017 721[96] 7271 3091 35866 4752 35126 14812 8223 21938 25919 15047 3202 18640 24198 3978 1477 4774 19031 15577 3619[115] NA 22176 31806 31619 27856 7748 18701 47017 721 7271 3091 35866 4752 35126 14812 8223 15286 NA NA[134] 7343 35709 1772 47363 17612 47517 27884 10894 33738 12984 12213 6667 20243 57546 31847 8433 494 3128 26497[153] 9 8075 7762 7214 20955 20801 1561 8172 28321 45015 2787 2624 8905 22741 5542 24392 36559 6202 37921[172] 9843 8141 12418 9260 3769 13302 3547 36109 19999 12472 15032 3069 38172 4415 13918 17315 2185 16245 32085[191] 7622 608 33943 6702 30975 27510 17475 13318 14267 49774 4913 6385 6171 11366 16855 17025 10871 50009 NA[210] 5846 18813 6885 38144 15944 2051 38495 17349 12413 14191 33361 21448 55495 34836 21206 15946 669 18683 22127[229] NA 57431 30426 9450
可以看到没有匹配上：
如果实在匹配不上，那就说明数据库里面没有这个基因，那我们就可以将这个基因舍弃掉。
由于此处是教学，我们就不一一匹配了。
下面将匹配到的基因的表达数据读出，代码如下：
alldiff <- allDiff[match(gene$Symbol,rownames(allDiff)),]head(alldiff)data <- data[match(rownames(alldiff),rownames(data)),]data <- na.omit(data)dim(data)data[1:5,1:5]write.csv(data,"Autophagy.csv")
这样我们的文件夹里面就多出来两个数据：
下一节将进行火山图，Venn图和热图的绘制。
一篇五分生信临床模型预测文章代码复现——Figure1 差异表达基因及预后基因筛选——下载数据（一）
如果觉得《一篇五分生信临床模型预测文章代码复现——Figure1 差异表达基因及预后基因筛选——》对你有帮助，请点赞、收藏，并留下你的观点哦！
