R语言画基因突变结构图
2024-05-12 16:34:10
R语言画基因突变结构图
做遗传病的同志们经常头痛的一个事应该是怎么画突变示意图，以前都是PPT直接画，但是最近碰到一个问题，综述里涉及到数个基因的数百个突变位点，PPT画的画得累死，于是开始搜索怎么用R来画基因突变结构图
1.gggenes包实现
gggenes包是一个用以可视化基因组或染色体基因结构的ggplot2扩展包
首先安装gggenes
#建议安装github版
devtools::install_github("wilkox/gggenes")
#以下代码参考自官方教程
library(gggenes)
library(tidyverse)
library(rtracklayer)
library(ggrepel)
rm(list=ls())## 基本原理就是分别提取基因及其对应外显子的位置坐标
## 分别使用gggenes的geom_gene_arrow和geom_subgene_arrow画出对应的图## 利用TxDb包获取基因组坐标信息
library(TxDb.Hsapiens.UCSC.hg38.knownGene)txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene## 提取外显子序列并去除重叠部分，提取基因序列
ex_db <- reduce(exonsBy(txdb, "gene"))
gene_db <- genes(txdb, single.strand.genes.only=F)## 输入你要提取的基因名
keys =“AMELX"
## 利用y叔的clusterProfiler包里的bitr函数进行基因ID转换
library(clusterProfiler)
library(org.Hs.eg.db)
ENTREZID <- bitr(keys, fromType = 'SYMBOL', toType = "ENTREZID", OrgDb = org.Hs.eg.db)## 提取感兴趣的基因的外显子信息，这里我只需要primary assembly版本
exons <- subset(as.data.frame(ex_db) , group_name %in% ENTREZID$ENTREZID & !str_detect(seqnames, 'alt'))
## 提取感兴趣的基因的信息，这里我只需要primary assembly版本genes <- subset(as.data.frame(gene_db), group_name %in% ENTREZID$ENTREZID & !str_detect(seqnames, 'alt'))exons <- merge(exons, ENTREZID, by.x = 'group_name', by.y = 'ENTREZID')
## 输出结果为
> exonsgroup_name group seqnames start end width strand SYMBOL
1 265 15832 chrX 11293413 11293468 56 + AMELX
2 265 15832 chrX 11294777 11294842 66 + AMELX
3 265 15832 chrX 11296779 11296826 48 + AMELX
4 265 15832 chrX 11298103 11298144 42 + AMELX
5 265 15832 chrX 11298236 11298277 42 + AMELX
6 265 15832 chrX 11298548 11298973 426 + AMELX
7 265 15832 chrX 11300607 11300761 155 + AMELXgenes <- merge(genes, ENTREZID, by.x = 'group_name', by.y = 'ENTREZID')## 输出结果为> genesgroup_name group seqnames start end width strand SYMBOL
1 265 15832 chrX 11293413 11300761 7349 + AMELX#开始作图library(ggsci)library(ggrepel)p <- ggplot(genes, aes(xmin = start, xmax = end, y = seqnames)) +geom_gene_arrow(aes(fill = SYMBOL)) +geom_gene_label(aes(label = SYMBOL)) +geom_subgene_arrow(data = exons,aes(xsubmin = start, xsubmax = end, fill = "royalblue")) + theme_genes() + scale_fill_brewer("Set3")+theme(legend.position = "none", axis.title = element_blank())
print(p)
## 然后是导入突变信息
## 一般我们从文献里整理或者从clinvar下载的突变信息是cDNA坐标的
## 我们可以利用transvar等工具将其转换为基因组坐标
## 下载整理后的突变位点信息如下
> AMELX_l
# A tibble: 27 × 4gene cDNA gDNA location <chr> <chr> <chr> <chr> 1 AMELX c.2T>C chrX:g.11294790T>C inside_[cds_in_exon_2]2 AMELX c.11G>A chrX:g.11294799G>A inside_[cds_in_exon_2]3 AMELX c.11G>C chrX:g.11294799G>C inside_[cds_in_exon_2]4 AMELX c.14_22delTTTTATTTG chrX:g.11294802_11294810delTTTTATTTG inside_[cds_in_exon_2]5 AMELX c.120T>C chrX:g.11298120T>C inside_[cds_in_exon_4]6 AMELX c.152C>T chrX:g.11298243C>T inside_[cds_in_exon_5]7 AMELX c.155C>G chrX:g.11298246C>G inside_[cds_in_exon_5]8 AMELX c.155delC chrX:g.11298246delC inside_[cds_in_exon_5]9 AMELX c.185delC chrX:g.11298276delC inside_[cds_in_exon_5]
10 AMELX c.152C>T chrX:g.11298555C>T inside_[cds_in_exon_5]
# … with 17 more rows
## 利用正则式提取一下突变坐标
## feature plot
feature_location <- str_replace_all( AMELX_l$gDNA, "(.+\\.)(\\d+)", "\\2") %>% str_remove_all(pattern = "[:alpha:]|>")
AMELX_l<- AMELX_l%>% dplyr::mutate(location = feature_location, .before =1) AMELX_l<- AMELX_l %>% dplyr::mutate(end = case_when(str_detect(location, "_") ~ str_replace (location, "(\\d+)_(\\d+)", "\\2"),.default = location), start = str_replace(location, "(\\d+)_(\\d+)", "\\1")) %>% modify_at(c("start", "end"), as.integer)
## 输出如下
>AMELX_l[,-4]
# A tibble: 27 × 5gene cDNA gDNA end start<chr> <chr> <chr> <int> <int>1 AMELX c.2T>C chrX:g.11294790T>C 11294790 112947902 AMELX c.11G>A chrX:g.11294799G>A 11294799 112947993 AMELX c.11G>C chrX:g.11294799G>C 11294799 112947994 AMELX c.14_22delTTTTATTTG chrX:g.11294802_11294810delTTTTATTTG 11294810 112948025 AMELX c.120T>C chrX:g.11298120T>C 11298120 112981206 AMELX c.152C>T chrX:g.11298243C>T 11298243 112982437 AMELX c.155C>G chrX:g.11298246C>G 11298246 112982468 AMELX c.155delC chrX:g.11298246delC 11298246 112982469 AMELX c.185delC chrX:g.11298276delC 11298276 11298276
10 AMELX c.152C>T chrX:g.11298555C>T 11298555 11298555
# … with 17 more rows
## 加入染色体信息
AMELX_l$seqnames <- "chrX"
AMELX_l <- AMELX_l%>% arrange(start)
## 抽取前10个突变作为示例
##利用ggrepel的函数来给突变加上label
## 主要调整参数包括force，nudge_x, nudge_y等等，
突变越多，调整难度越大，100个以上建议还是R生成SVG，然后用AI手调
#具体参见官方文档https://ggrepel.slowkow.com/articles/examples.html#examples
p + geom_text_repel(data = AMELX_l[1:10,], aes(x = start, label = cDNA, ), size =2,force = 20, # 强制让各个label不重叠direction = "both",hjust = 0.2,segment.size = 0.2,max.iter = 1e4, max.time = 1, max.overlaps = Inf, nudge_x =1000, nudge_y = 0.5)
transPlotR
除此之外还可以使用junjunlab开发的transPlotR和BioSeqUtils来达到类似的效果，具体参见https://github.com/junjunlab/BioSeqUtils
mytest <- loadGenomeGTF(gtfPath = "hg38.ncbiRefSeq.gtf.gz",genomePath = "hg38.fa.gz")
library(transPlotR)keys = c("AMELX","ENAM","FAM20A","FAM83H")
# get intron
intron <- getIntronInfo(mytest,geneName = keys)# get gtf
gtfFile = data.frame(mytest@gtf)
gtfFile <- rbind(gtfFile,intron)# plot
p1 <- trancriptVis(gtfFile = gtfFile,gene = keys[2], selecType = "lt", exonColorBy = "type", #facetByGene = T, exonWidth = 0.1, textLabel = "gene_name",arrowCol = NA,addNormalArrow = F,newStyleArrow = T,show.legend = T) +ggsci::scale_fill_igv() + scale_x_continuous(expand = expansion(mult= c(0.08,0.05))) + guides(fill = guide_legend("Type"))
同样用ggrepel添加一下突变信息
## feature plot
##导入整理好的突变信息
ENAM_l <- read_csv("ENAM_mutation_location.csv")
feature_location <- str_replace_all(ENAM_l$gDNA, "(.+\\.)(\\d+)", "\\2") %>% str_remove_all(pattern = "[:alpha:]|>")
ENAM_l <- ENAM_l %>% dplyr::mutate(location = feature_location, .before =1) ENAM_l <- ENAM_l %>% dplyr::mutate(end = case_when(str_detect(location, "_") ~ str_replace(location, "(\\d+)_(\\d+)", "\\2"),.default = location
), start = str_replace(location, "(\\d+)_(\\d+)", "\\1")) %>% modify_at(c("start", "end"), as.integer)
ENAM_l$seqnames <- "chr4"
ENAM_l <- ENAM_l %>% arrange(start)
library(ggrepel)
## 这里y坐标根据transPlot的内含子的y坐标，相对比较好看
p1 + geom_text_repel(data = ENAM_l[1:10,], aes(x = start, y= 1.016667, label = cDNA, ), size =2,force = 40, # do not pull toward data pointsdirection = "x",hjust = 0.3,segment.size = 0.2,max.iter = 1e4, max.time = 1, max.overlaps = Inf, nudge_x =1000, nudge_y = 0.2)
ggsave("ENAM_mut.jpg", units = "cm", width = 30, height = 15, device = "jpg")
完成，这样至少能在突变少于等于30个的时候比较自动化了，当然还是要手调…
R语言画基因突变结构图相关推荐
- 利用R语言画简单时间序列图
R 语言无法自动将读取的数据转化为时间序列格式, 所以利用R语言画时间序列图的一个关键步骤就是将读取的数据转变为时间序列格式, 下面是一个简单的程序: # 读取数据, 首先将excel 格式的转化为 ...
- R语言画个中国地图使用shp文件
R语言画个中国地图,使用shp文件 前几天帮人用R语言画个一个河北省的地图,河北省各个市被填充上颜色,颜色的深浅和这个市的GDP有关系,效果如下: 然后加上各个城市的名字并加上指北针.再把背景的网格去 ...
- R语言画森林图方法4
获取更多R语言知识,请关注公众号:医学和生信笔记 医学和生信笔记,专注R语言在临床医学中的使用,R语言数据分析和可视化.主要分享R语言做医学统计学.meta分析.网络药理学.临床预测模型.机器学习.生 ...
- 【R语言画生日蛋糕】
R语言画生日蛋糕 直接上代码! rm(list = ls())a100 <- seq(0, 100, by = 1) a50 <- seq(12, 88, length.out = 101 ...
- 用R语言画一朵玫瑰花
首先,需要安装并加载绘图包 "ggplot2". install.packages("ggplot2") library(ggplot2) 然后,使用函数 gg ...
- 用r语言画出y = ax^2 + bx + c,R语言中如何使用最小二乘法
这里只是介绍下R语言中如何使用最小二乘法解决一次函数的线性回归问题. 代码如下: > x > y > lsfit(x,y) 结果如下: $coefficients Intercept ...
- R语言画词云图——建模常用软件
在数学建模--软件篇介绍了我常用的软件,借着今年的华为杯,想做个LOGO.然后就使用R语言做了词云图,不过最后由于队友做得比较好,就没采用.本文吧效果和代码贴出来记录一下. library(wordc ...
- 包包的结构制图_原来可以用R这么画基因结构图
gggenes是一款基于ggplot2开发的R包,可以很方便的画出下图所示的基因结构图. 1. 安装R包 直接从CRAN官方源来安装: install.packages("gggenes&q ...
- R语言对基因进行GO注释（附代码）
通常在挑选出一些基因之后,需要对这些基因进行GO注释,分析该基因的功能,现在把R分析GO的代码放在下面 所需要的文件,有一列geneid即可 rm(list=ls()) genelist <- ...
最新文章
- 计算机的来源知识,如何理解计算机知识及计算机发展史
- Hadoop大数据--Mapreduce程序运行并发度
- Mybatis注解实现一对多关联映射(@Many)
- python之socket
- python第五次作业——陈灵院
- 开源数据库中间件-MyCa初探与分片实践
- 【论文翻译】学习新闻事件预测的因果关系
- 资管新规这样规定，我的货基该怎么办？
- everything搭配什么软件_如果你在用Everything，那这个插件你一定会毫不犹豫就装上！...
- Pr 入门教程：如何向影片中的剪辑添加过渡效果？
- html 菜单栏横向排列,响应式可弹出横向导航栏
- python 破解验证码之二：OCR tesseract识别验证码
- 使用Houdini快速将图片转换成文字模型
- 使用Windows服务启动C#桌面应用程序问题解决
- Promise晋级—完全吃透
- Shiro 通过配置Cookie 解决多个二级域名的单点登录问题。
- oracle rfs进程过多,【DB笔试面试755】在Oracle的DG中，RFS、LNSn、MRP、LSP进程的作用分别是什么？...
- [图]WPS Office 2019上架微软商城：引入全新用户界面
- 【规范化】项目规范化
- error: unrecognized arguments:问题