RNA 8. SCI文章中差异基因表达--热图 (heatmap)

2023-11-06 06:59
文章标签 文章 sci 表达 heatmap rna 差异基因 热图

本文主要是介绍RNA 8. SCI文章中差异基因表达--热图 (heatmap)，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

前言

大多我们在做完差异表达之后都会看下我们的差异基因筛选的是否能将分组结果展现出来，都会选择热图，主要是热图技能聚类，又可以展现表达量的大小，非常直观，所以这期我们就说下热图的绘制方法。

实例解析

1. 数据读取

数据的读取我们仍然使用的是 TCGA-COAD 的数据集，表达数据的读取以及临床信息分组的获得我们上期已经提过，我们使用的是edgeR 软件包计算出来的差异表达结果，合并了原始的 Count 数据，如下：

DEG=read.table("DEG-resdata.xls",sep="\t",check.names=F,header = T)
PCAmat<-DEG[,8:ncol(DEG)]
rownames(PCAmat)=DEG[,1]
PCAmat[1:5,1:3]##                 TCGA-3L-AA1B-01A-11R-A37K-07 TCGA-4N-A93T-01A-11R-A37K-07
## ENSG00000142959                           20                           15
## ENSG00000163815                          175                          108
## ENSG00000107611                           49                           13
## ENSG00000162461                           55                           89
## ENSG00000163959                          153                          259
##                 TCGA-4T-AA8H-01A-11R-A41B-07
## ENSG00000142959                           49
## ENSG00000163815                           59
## ENSG00000107611                            6
## ENSG00000162461                           48
## ENSG00000163959                          102

读取分组信息癌组织样本 478 和癌旁组织样本 41，如下：

group<-read.table("DEG-group.xls",sep="\t",check.names=F,header = T)
head(group)##                         Sample Group
## 1 TCGA-D5-6530-01A-11R-1723-07    TP
## 2 TCGA-G4-6320-01A-11R-1723-07    TP
## 3 TCGA-AD-6888-01A-11R-1928-07    TP
## 4 TCGA-CK-6747-01A-11R-1839-07    TP
## 5 TCGA-AA-3975-01A-01R-1022-07    TP
## 6 TCGA-A6-6780-01A-11R-1839-07    TP

table(group$Group)## 
##  NT  TP 
##  41 478

2. Count 转换为 FPKM

我们之前下载的数据是 Count 用来做差异表达分析，而在绘制热图时需要我们使用RPKM或者时FPKM，那么我们利用差异基因的结果，提取对应的基因，得到Count 矩阵，下载基因长度，来计算FPKM即可，还记得我们之前讲过的FPKM的计算公式，

FPKM (Fragments Per Kilobase of exon model per Million mapped fragments)

  1. FPKM：每千个碱基的转录每百万映射读取的fragments，主要是针对pair-end测序表达量进行计算。

  2. FPKM和RPKM的区别就是一个是fragment，一个是read。

那么我们首先下载外显子的注释信息，GDC Reference Files | NCI Genomic Data Commons (cancer.gov) 找到如下信息，下载即可。

外显子注释文件的处理，我们希望得到第一列时基因，第二列时基因长度。

安装并加载软件包 GenomicFeatures，这个软件包在处理注释文件时非常方便，如下：

if(!require(GenomicFeatures)){BiocManager::install("GenomicFeatures")
}library(GenomicFeatures)

读取下载的注释文件，并获取外显子的总长度，如下：

txdb <- makeTxDbFromGFF("gencode.v22.annotation.gtf.gz",format = "gtf")
exons <- exonsBy(txdb,by="gene")exons_length<-lapply(exons,function(x){sum(width(reduce(x)))})
class(exons_length)## [1] "list"length(exons_length)## [1] 60483exons_length<-as.data.frame(exons_length)
dim(exons_length)## [1]     1 60483exons_length1<-t(exons_length)
exons_length1<-as.data.frame(exons_length1)
dim(exons_length1)## [1] 60483     1head(exons_length1)##                      V1
## ENSG00000000003.13 4535
## ENSG00000000005.5  1610
## ENSG00000000419.11 1207
## ENSG00000000457.12 6883
## ENSG00000000460.15 5967
## ENSG00000000938.11 3474

我们发现差异基因后面时不带 “.” 和数字的，所以需要我们去掉基因后面的数字，如下：

colnames(exons_length1)="Length"
Gene<-gsub("\\.(\\.?\\d*)","",rownames(exons_length1))
exons_length1$Gene=as.data.frame(Gene)[,1]
head(exons_length1)##                    Length            Gene
## ENSG00000000003.13   4535 ENSG00000000003
## ENSG00000000005.5    1610 ENSG00000000005
## ENSG00000000419.11   1207 ENSG00000000419
## ENSG00000000457.12   6883 ENSG00000000457
## ENSG00000000460.15   5967 ENSG00000000460
## ENSG00000000938.11   3474 ENSG00000000938

根据提取差异基因并与 Count 矩阵合并，矩阵最后一列时基因长度，如下：

PCAmat$Gene<-rownames(PCAmat)
####count 矩阵，长度合并
count_length<-merge(PCAmat,exons_length1,by="Gene")

计算 FPKM，利用合并后的结果，样本有519个，那么就是从2列到520列，结合公式开始计算 FPKM，如下：

kb<-count_length$Length/1000
countdata<-count_length[,2:520]
rpk <- countdata/kb
fpkm <-t(t(rpk)/colSums(countdata) * 10^6)
rownames(fpkm)=rownames(PCAmat)
fpkm[1:5,1:3]##                 TCGA-3L-AA1B-01A-11R-A37K-07 TCGA-4N-A93T-01A-11R-A37K-07
## ENSG00000142959                    1.5336078                     1.826404
## ENSG00000163815                   51.1904250                   388.921103
## ENSG00000107611                  162.7697536                    52.262338
## ENSG00000162461                   69.9931252                    93.087394
## ENSG00000163959                    0.2628356                     0.000000
##                 TCGA-4T-AA8H-01A-11R-A41B-07
## ENSG00000142959                    0.7546987
## ENSG00000163815                  178.9371541
## ENSG00000107611                   13.3500151
## ENSG00000162461                   90.5084641
## ENSG00000163959                    0.0000000

FPKM 数据进行标准，如下：

data=log2(fpkm+1)
dat=t(scale(t(data))) # 'scale'可以对log(fpkm+1)数值进行归一化
#处理数据
dat[dat>2]=2
dat[dat<(-2)]= -2
dim(dat)## [1] 4128  519dat[1:5,1:3]##                 TCGA-3L-AA1B-01A-11R-A37K-07 TCGA-4N-A93T-01A-11R-A37K-07
## ENSG00000142959                  0.312929749                   0.51326990
## ENSG00000163815                 -0.451067532                   1.30683599
## ENSG00000107611                  0.811371990                   0.01667015
## ENSG00000162461                 -0.076316888                   0.23114035
## ENSG00000163959                 -0.004643292                  -0.62802639
##                 TCGA-4T-AA8H-01A-11R-A41B-07
## ENSG00000142959                   -0.3600191
## ENSG00000163815                    0.6308436
## ENSG00000107611                   -0.9112186
## ENSG00000162461                    0.2008002
## ENSG00000163959                   -0.6280264

3. 绘制热图

热图绘制使用软件包 pheatmap，这个软件包使用起来非常方便，参数多能够满足各种需求，更重要的是提取聚类数据也可以自己二次处理聚类结果，注释方便。

安装并加载软件包，如下：

if(!require(pheatmap)){install.packages("pheatmap")
}library(pheatmap)

最简单的一种方式绘制热图，如下：

###########pheatmap 
#调整参数，换颜色
pheatmap(dat,cluster_rows = TRUE,show_rownames=FALSE,show_colnames = FALSE,scale="none",cluster_cols = F,fontsize_row = 10,fontsize_col = 10,#color = colorRampPalette(c("navy", "white", "firebrick3"))(100),color = colorRampPalette(c("green3", "white", "blue4"))(100),#换颜色angle_col = 45, #修改横轴坐标名倾斜度filename = 'cor.fpkm1.png',
)

图片

添加加注释信息，需要对分组信息做一下出来，

首先是样本的分组注释，根据癌症组和正常组整理，也就是列注释数据处理，如下：

#######加注释信息
rownames(group)=group$Sample
head(group)##                                                    Sample Group
## TCGA-D5-6530-01A-11R-1723-07 TCGA-D5-6530-01A-11R-1723-07    TP
## TCGA-G4-6320-01A-11R-1723-07 TCGA-G4-6320-01A-11R-1723-07    TP
## TCGA-AD-6888-01A-11R-1928-07 TCGA-AD-6888-01A-11R-1928-07    TP
## TCGA-CK-6747-01A-11R-1839-07 TCGA-CK-6747-01A-11R-1839-07    TP
## TCGA-AA-3975-01A-01R-1022-07 TCGA-AA-3975-01A-01R-1022-07    TP
## TCGA-A6-6780-01A-11R-1839-07 TCGA-A6-6780-01A-11R-1839-07    TPgroup1=group[colnames(dat),-1,drop=FALSE]
head(group1)##                              Group
## TCGA-3L-AA1B-01A-11R-A37K-07    TP
## TCGA-4N-A93T-01A-11R-A37K-07    TP
## TCGA-4T-AA8H-01A-11R-A41B-07    TP
## TCGA-5M-AAT4-01A-11R-A41B-07    TP
## TCGA-5M-AAT5-01A-21R-A41B-07    TP
## TCGA-5M-AAT6-01A-11R-A41B-07    TP

加入列注释，再进行绘制热图，如下：

pheatmap(dat,cluster_rows = FALSE,show_rownames=FALSE,show_colnames = FALSE,scale="none",cluster_cols = FALSE,fontsize_row = 10,fontsize_col = 10,annotation_col = group1,color = colorRampPalette(c("navy", "white", "firebrick3"))(100),#color = colorRampPalette(c("green3", "white", "blue4"))(100),#换颜色angle_col = 45, #修改横轴坐标名倾斜度filename = 'cor.fpkm2.png',
)

图片

根据差异结果中 Up, Down 分组来调整基因的顺序，在对行进行处理，这其实就是对基因的一个分组，分为Up 和 Down，如下：

ann_row=as.data.frame(DEG[,c(1,7)])
colnames(ann_row)=c("Sample","Sig")
rownames(ann_row)=DEG$Row.names
ann_row=ann_row[order(ann_row$Sig),]
rownames(ann_row)=ann_row$Sample
ann_row=as.data.frame(ann_row)
ann_row1=as.data.frame(ann_row[,2])
rownames(ann_row1)=ann_row$Sample
colnames(ann_row1)="Sig"
head(ann_row1)##                  Sig
## ENSG00000142959 Down
## ENSG00000163815 Down
## ENSG00000107611 Down
## ENSG00000162461 Down
## ENSG00000163959 Down
## ENSG00000144410 Down

加入行列注释，再进行绘制热图，如下：

#######根据差异结果中 Up, Down 分组来调整基因的顺序
pheatmap(dat[rownames(ann_row1),],cluster_rows = FALSE,show_rownames=FALSE,show_colnames = FALSE,scale="none",cluster_cols = FALSE,fontsize_row = 10,fontsize_col = 10,annotation_col = group1,annotation_row = ann_row1,color = colorRampPalette(c("navy", "white", "firebrick3"))(100),#color = colorRampPalette(c("green3", "white", "blue4"))(100),#换颜色angle_col = 45, #修改横轴坐标名倾斜度filename = 'cor.fpkm3.png',
)

图片

拆分多个模块，我们先看下基因上下调个数，下调基因数为1296，如下

table(ann_row1$Sig)## 
## Down   Up 
## 1296 2832

根据下调基因的个数，设置 gaps_row=1296，如下：

#########拆分多个模块
pheatmap(dat[rownames(ann_row1),],cluster_rows = FALSE,show_rownames=FALSE,show_colnames = FALSE,scale="none",cluster_cols = FALSE,fontsize_row = 10,fontsize_col = 10,annotation_col = group1,annotation_row = ann_row1,color = colorRampPalette(c("navy", "white", "firebrick3"))(100),#color = colorRampPalette(c("green3", "white", "blue4"))(100),#换颜色angle_col = 45, #修改横轴坐标名倾斜度filename = 'cor.fpkm4.png',gaps_row=1296#拆分位置
)

图片

完成，非常简单，但是也需要大家仔细认真地把数据准备好，关注公众号，扫码进群，精彩内容不断更新！

在这里插入图片描述

桓峰基因

生物信息分析，SCI文章撰写及生物信息基础知识学习：R语言学习，perl基础编程，linux系统命令，Python遇见更好的你

45篇原创内容

公众号

这篇关于RNA 8. SCI文章中差异基因表达--热图 (heatmap)的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！


http://www.chinasem.cn/article/355115。

相关文章

java计算机毕设课设—停车管理信息系统(附源码、文章、相关截图、部署视频)

java计算机毕设课设—停车管理信息系统(附源码、文章、相关截图、部署视频)

这是什么系统？ 资源获取方式在最下方 java计算机毕设课设—停车管理信息系统(附源码、文章、相关截图、部署视频) 停车管理信息系统是为了提升停车场的运营效率和管理水平而设计的综合性平台。系统涵盖用户信息管理、车位管理、收费管理、违规车辆处理等多个功能模块，旨在实现对停车场资源的高效配置和实时监控。此外，系统还提供了资讯管理和统计查询功能，帮助管理者及时发布信息并进行数据分析，为停车场的科学
阅读更多...
文章解读与仿真程序复现思路——电力自动化设备EI\CSCD\北大核心《考虑燃料电池和电解槽虚拟惯量支撑的电力系统优化调度方法》

文章解读与仿真程序复现思路——电力自动化设备EI\CSCD\北大核心《考虑燃料电池和电解槽虚拟惯量支撑的电力系统优化调度方法》

本专栏栏目提供文章与程序复现思路，具体已有的论文与论文源程序可翻阅本博主免费的专栏栏目《论文与完整程序》 论文与完整源程序_电网论文源程序的博客-CSDN博客https://blog.csdn.net/liang674027206/category_12531414.html 电网论文源程序-CSDN博客电网论文源程序擅长文章解读,论文与完整源程序,等方面的知识,电网论文源程序关注python
阅读更多...
【Linux】萌新看过来！一篇文章带你走进Linux世界

【Linux】萌新看过来！一篇文章带你走进Linux世界

🚀个人主页：奋斗的小羊 🚀所属专栏：Linux 很荣幸您能阅读我的文章，诚请评论指点，欢迎欢迎 ~ 目录 前言💥1、初识Linux💥1.1 什么是操作系统？💥1.2 各种操作系统对比💥1.3 现代Linux应用💥1.4 Linux常用版本 💥2、Linux 和 Windows 目录结构对比💥2.1 文件系统组织方式💥2.2
阅读更多...
多线程的系列文章

多线程的系列文章

Java多线程学习（一）Java多线程入门 Java多线程学习（二）synchronized关键字（1）   Java多线程学习（二）synchronized关键字（2） Java多线程学习（三）volatile关键字 Java多线程学习（四）等待/通知（wait/notify）机制 Java多线程学习（五）线程间通信知识点补充 Java多线程学习（六）Lock锁的使用 Java多
阅读更多...
缓存的常见问题 以及解决博客文章

缓存的常见问题 以及解决博客文章

1.jedispool 连 redis 高并发卡死  （子非鱼yy） https://blog.csdn.net/ztx114/article/details/78291734 2. Redis安装及主从配置 https://blog.csdn.net/ztx114/article/details/78320193 3.Spring中使用RedisTemplate操作Redis（sprin
阅读更多...
java计算机毕设课设—企业员工信息管理系统(附源码、文章、相关截图、部署视频)

java计算机毕设课设—企业员工信息管理系统(附源码、文章、相关截图、部署视频)

这是什么系统？ 获取资料方式在最下方 java计算机毕设课设—企业员工信息管理系统(附源码、文章、相关截图、部署视频) 企业员工信息管理系统旨在为公司提供高效的员工信息管理解决方案。该系统的核心功能涵盖密码修改、员工管理、部门管理、出勤管理、工资管理、请假审核等方面，帮助企业优化人力资源管理流程。系统结构如下： （1）前端（员工端）： 1.密码修改：员工可以修改自己的密码，提升账户的安全
阅读更多...
android的工程和代码的命名规范（第一篇文章，勿喷）

android的工程和代码的命名规范（第一篇文章，勿喷）

1。首先我们从编译代码的工具说起吧：工程中的注释一般都是中文写的（毕竟大家都是中国人，还是习惯于中文）这样就设计到乱码的问题了；对于这类问题，我们一般最好的处理方法就是将工程设置成 UTF-8 的格式；下面就说说怎么将工作空间或者是工程设置成UTF-8 的格式吧(当然我这里面说的是eclips
阅读更多...
C#/.NET/.NET Core推荐学习路线文档文章

C#/.NET/.NET Core推荐学习路线文档文章

前言 专门为C#/.NET/.NET Core推荐学习路线&文档&文章提供的一个Issues，各位小伙伴可以把自己觉得不错的学习路线、文档、文章相关地址分享出来🤞。 https://github.com/YSGStudyHards/DotNetGuide/issues/10 🏷️C#/.NET/.NET Core优质学习资料 📚.NET 入门教程 📚
阅读更多...
【java 走进NLP】simhash 算法计算两篇文章相似度

【java 走进NLP】simhash 算法计算两篇文章相似度

python 计算两篇文章的相似度算法simhash见： https://blog.csdn.net/u013421629/article/details/85052915 对长文本 是比较合适的（超过500字以上） 下面贴上java 版本实现： pom.xml 加入依赖 <dependency><groupId>org.jsoup</groupId><artifactId>jsoup</a
阅读更多...
【python 走进NLP】simhash 算法计算两篇文章相似度

【python 走进NLP】simhash 算法计算两篇文章相似度

互联网网页存在大量的重复内容网页，无论对于搜索引擎的网页去重和过滤、新闻小说等内容网站的内容反盗版和追踪，还是社交媒体等文本去重和聚类，都需要对网页或者文本进行去重和过滤。最简单的文本相似性计算方法可以利用空间向量模型，计算分词后的文本的特征向量的相似性，这种方法存在效率的严重弊端，无法针对海量的文本进行两两的相似性判断。模仿生物学指纹的特点，对每个文本构造一个指纹，来作为该文本的标识，从形式上来
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2