之前讲过临床模型预测的专栏,但那只是基础版本,下面我们以自噬相关基因为例子,模仿一篇五分文章,将图和代码复现出来,学会本专栏课程,可以具备发一篇五分左右文章的水平:
本专栏目录如下:
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格式,并放在我的资源中,可以直接下载使用:
该数据里面包含了表达数据和临床数据,如果大家想用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")
这样我们的文件夹里面就多出来两个数据:
下一节将进行火山图,Venn图和热图的绘制。