把基因画在染色体上

参考简书：链接：https://www.jianshu.com/p/e924e33c60c9

需要六个R包：RIdeogram,RColorBrewer,ggbio,gggenes,Gviz,ggplot2

library(RIdeogram)
library(RColorBrewer)
library(ggbio)
library(gggenes)
library(Gviz)
library(ggplot2)
#1、准备染色体长度文件
DmChromosome <- data.table::fread("grape38_chr_len",data.table=FALSE) 
#2、gene密度信息
Dmdensity <- GFFex(input="grape38.gff3",
                   karyotype="grape38_chr_len",
                   feature="gene",window=10000)
#input为gff文件类型，karyotype染色体文件。
#Dmdensity文件最后为一个四列文件：染色体，染色体起始位置，染色体终止位置，以及这段windows中的gene个数。

#marker型绘图，有marker/line/polygon/heatmap四种绘图方法。
#3、marker型需要Type/ Shape/ Chr/ Start/ End/ color。
#它仅对type一栏绘制图例，marker型只适合展示少量基因,如下展示7个基因。
gene_list<- data.table::fread('geneinfo')
head(gene_list)
ideogram(karyotype=DmChromosome, #染色体信息
               overlaid=Dmdensity, #基因密度
               label=gene_list, #目标基因标签
               label_type="marker",width=135)

convertSVG("chromosome.svg",device="png")

#4、line/polygon/heatmap绘图
#需要制作一个类似于基因密度信息的表，包含Chr/ Start/ End/ Value_1/ Color_1/ Value_2/ Color_2。
gene_list_v<- data.table::fread('SNPdensity')
ideogram(karyotype=DmChromosome, #染色体信息
                   overlaid=Dmdensity, #基因密度
                   label=gene_list_v, #目标基因标签
                   label_type="line",width=135)
convertSVG("chromosome.svg",device="png")

#heatmap型绘图
ideogram(karyotype=DmChromosome,
         overlaid=Dmdensity,
         label=gene_list_v,
         label_type="heatmap",width=135,
         colorset1 = c("#f7f7f7", "#2c7fb8"),
         colorset2 = c("#f7f7f7", "#e34a33"))
convertSVG("chromosome.svg",device="png")

heatmap图没有line图清楚。

二、gggenes,把具体基因画在染色体上

head(example_genes)
dim(example_genes)
p1 <- ggplot(example_genes,aes(xmin=start,xmax=end,y=molecule,fill=gene)) +
     geom_gene_arrow() +
     facet_wrap(~ molecule,scales="free",ncol=1) +
     scale_fill_brewer(palette="Set3") +
     theme_genes()
p1

#对齐某一基因
dummies <- make_alignment_dummies(example_genes,aes(xmin=start,xmax=end,y=molecule,id=gene),on="genE")
p1 + geom_blank(data=dummies)

#写上基因名
p2 <- ggplot(example_genes,aes(xmin=start,xmax=end,y=molecule,fill=gene)) +
     geom_gene_arrow(arrowhead_height=unit(3,"mm"),arrowhead_width=unit(1,"mm")) +
     facet_wrap(~ molecule,scales="free",ncol=1) +
     scale_fill_brewer(palette="Set3") +
     theme_genes()
p2
p2 + geom_gene_label(aes(label=gene),align="left")

#区分正负链
p3 <- ggplot(example_genes,aes(xmin=start,xmax=end,y=molecule, fill=gene,forward=direction)) +
     geom_gene_arrow() +
     facet_wrap(~ molecule,scales="free",ncol=1) +
     scale_fill_brewer(palette="Set3") +
     theme_genes()
p3

#基因亚结构
head(example_subgenes);dim(example_subgenes)
p4 <- ggplot(example_genes,aes(xmin=start,xmax=end,y=molecule)) +
     facet_wrap(~ molecule, scales="free",ncol = 1) +
     geom_gene_arrow(fill="white") +
     geom_subgene_arrow(data=example_subgenes, aes(xmin=start,xmax=end,y=molecule,fill=gene,xsubmin=from,xsubmax=to),color="black",alpha=.7) +
     theme_genes()
p4