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')
## 输出结果为
> exons
group_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 1