一篇五分生信临床模型预测文章代码复现——Figure1 差异表达基因及预后基因筛选——下载数据（一）
之前讲过临床模型预测的专栏，但那只是基础版本，下面我们以自噬相关基因为例子，模仿一篇五分文章，将图和代码复现出来，学会本专栏课程，可以具备发一篇五分左右文章的水平：
本专栏目录如下：
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")这样我们的文件夹里面就多出来两个数据：
下一节将进行火山图，Venn图和热图的绘制。
一篇五分生信临床模型预测文章代码复现——Figure1 差异表达基因及预后基因筛选——下载数据（一）相关推荐
- 一篇五分生信临床模型预测文章代码复现——Figure 7 外部数据集验证模型
之前讲过临床模型预测的专栏,但那只是基础版本,下面我们以自噬相关基因为例子,模仿一篇五分文章,将图和代码复现出来,学会本专栏课程,可以具备发一篇五分左右文章的水平: 本专栏目录如下: Figure 1 ...
- 五分钟学会用Simulink模型生成HDL代码
五分钟学会用Simulink模型生成HDL代码 1 核心步骤 2 视频展示 3 生成HDL代码的注意事项 3.1 HDL支持的库和模块 3.2 设置simulink模型为可生成 hdl 的模式 3.3 ...
- R语言与临床模型预测——LASSO回归，单因素多因素cox，差异表达分析，Venn图，森林图，列线图，矫正曲线，ROC全套代码及解析——第十三部分 校准曲线 本专栏可免费答疑
1.下载数据 2. 匹配基因 3. 基因去重复 4.匹配临床数据 5.批量cox回归分析 6.差异表达基因筛选 7.取交集,选出预后相关的差异表达基因 8.森林图绘制 9.lasso回归进一步排除具有 ...
- 生信技能-高通量测序工具bam、samtools、bedtools及conda的下载和安装
一.BWA 1.介绍 简介:用于建立 index:基于 BWT 算法,将 reads 比对到参考基因组:最新版本 bwa-mem2,Intel实验室对计算效率进行了优化. 详情:baw是一款将序列比对 ...
- 重磅！这个生信神器助你文章秒出图——miRNA与基因互作数据库
我们熟知,在特定情况下,microRNA(miRNA)可以直接或间接激活和抑制基因表达.但是,尚没有基于多组学的数据库能够证明对激活与抑制以及正常与癌症状况之间相互作用模式转换的系统数据.今天我们为大 ...
- CV领域Transformer这一篇就够了（原理详解+pytorch代码复现）
文章目录 前言 一.注意力机制 1.1注意力机制通俗理解 1.2注意力机制计算公式 1.3注意力机制计算过程 1.4注意力机制代码 二.自注意力机制 2.1 注意力机制和自注意力机制的区别 2.2 编 ...
- 生信宝典文章集锦，一站式学习生信！众多干货，有趣有料
生信的作用越来越大,想学的人越来越多,不管是为了以后发展,还是为了解决眼下的问题.但生信学习不是一朝一夕就可以完成的事情,也许你可以很短时间学会一个交互式软件的操作,却不能看完程序教学视频后就直接写程 ...
- 生信宝典教程大放送，一站式学习生信技术
生物信息学包含生物数据分析.数据可视化.重复工作程序化,是生物.医学科研必备的技能之一.生信宝典精心组织生信学习系列教程.生信工具精品教程,通过大量的生信例子.关键的注释.浓缩的语句和录制的视频帮助快 ...
- 送书 | 知乎阅读300w+的生信学习指南（更新版）
先送书 在上周的留言送书活动中,恭喜下面这位读者获得书籍"Oracle高性能系统架构实战大全",请及时与生信宝典编辑(shengxinbaodian)联系. 2020过去三分之一了 ...
最新文章
- 电子秤专用模拟/数字（A/D）转换器芯片 HX711
- 第K短路模板【POJ2449 / 洛谷2483 / BZOJ1975 / HDU6181】
- Java 8 forEach 示例
- OpenCASCADE：扩展数据交换(XDE)的简介
- 【枭·音频】注入灵魂—《暗影火炬城》角色语音后期处理
- datagrid页面获取表单一条数据的例子
- Linux安装mysql（解决E: Package ‘mysql-server‘ has no installation candidate与ERROR 1698 (28000)）
- 正则表达式五分钟快速复习
- jar反编译工具 比jd-gui 功能更强大的 Luyten 查看jar源码， 解决jd反编译代码中break labelxxx 、 static初始块中出现return 等问题
- 中科大一所学校撑起中国人工智能半壁江山
- DASCTF2022 7月赋能赛 crypto wp（DASCTF2022.07赋能赛Pwn easyheap）
- 如何 自定义starter?
- appcan使用心得体会
- COVID-19 抗原自检试剂盒行业研究及十四五规划分析报告
- Android案例分享__HomePageA__仿'58到家/百度糯米/豆果美食/美团外卖/手机京东'首页
- 软件测试基础-今日②问-4
- Codeforces 1293 E. Xenon‘s Attack on the Gangs —— 树上记忆化搜索，单点加改成区间加，有丶东西
- 我的大三一年职业规划，预期毕业目标
- yolov3识别的类别_YOLO v3实战之钢筋数量AI识别（一）
- Springboot 项目导出word文档(文档内容包括数据以及服务器图片)