一篇五分生信临床模型预测文章代码复现——Figure1 差异表达基因及预后基因筛选——下载数据（一）

楷然教你学生信

已于 2022-09-20 10:34:15 修改

阅读量1.7k

点赞数 2

分类专栏：一篇五分临床模型预测文章代码复现文章标签： python 开发语言

于 2022-08-14 13:35:08 首次发布

本文链接：https://blog.csdn.net/weixin_46500027/article/details/126329244

版权

一篇五分临床模型预测文章代码复现专栏收录该内容

36 篇文章 7 订阅 ¥299.90 ¥99.00

订阅专栏

超级会员免费看

之前讲过临床模型预测的专栏，但那只是基础版本，下面我们以自噬相关基因为例子，模仿一篇五分文章，将图和代码复现出来，学会本专栏课程，可以具备发一篇五分左右文章的水平：

本专栏目录如下：

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下载频道.https://download.csdn.net/download/weixin_46500027/85440520?spm=1001.2014.3001.5503

该数据里面包含了表达数据和临床数据，如果大家想用TPM格式的数据，可以到官网下载后自己转化格式，转化格式的代码在我的另一个专栏有相关介绍：

https://blog.csdn.net/weixin_46500027/category_11845300.html?spm=1001.2014.3001.5482https://blog.csdn.net/weixin_46500027/category_11845300.html?spm=1001.2014.3001.5482下面我们开始下载自噬相关基因，这是我们研究的某个方向，大家可以自己去找别的研究方向，模仿本专栏提供的课程和代码。

自噬相关基因下载如下：

自噬相关基因有一个数据库，在我们之前的文章中有讲过：

注意：由于自噬数据库最近打不开了，所以我将自噬相关基因上传到百度网盘上，供大家练习：

链接：https://pan.baidu.com/s/1dgUVqcoVoQvmLvQtKukSQA
提取码：z9m2
--来自百度网盘超级会员V5的分享

简单来说呢，就是把这些基因复制，粘贴到Excel里面。

全选，复制：

粘贴完了以后，我们需要对空格进行一下处理，因为自噬数据库呢，没有下载的地方我们只能这样手动粘贴，然后手动删减。

总之，处理好了以后，就剩下三列，我们最主要的就是要Symbol这一列。

注意：由于自噬数据库最近打不开了，所以我将自噬相关基因上传到百度网盘上，供大家练习：

链接：https://pan.baidu.com/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下载频道.https://download.csdn.net/download/weixin_46500027/84997590?spm=1001.2014.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)
grep
tumor <- 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)
grep1
normal <- 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))
design
colnames(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,]

#          X                              Name 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] 20109  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")

这样我们的文件夹里面就多出来两个数据：