使用RIdeogram绘制染色体图谱

参考文章：https://www.jianshu.com/p/07ae1fe18071
需求描述：我想绘制差异表达基因在染色体上的分布，或开放染色质在染色体上的分布或
CTCF在染色体上的结合位点或突变位点在染色体上的分布或DNA甲基化在染色体上的分布

library("RIdeogram")
library(stringr)
data(human_karyotype, package="RIdeogram")
#包含Chr Start       End  CE_start    CE_end 五列
#karyotype文件
# 第一列：染色体号
# 第二列：起始
# 第三列：终止
# 第四列：中心粒起始位置
# 第五列：中心粒终止位置
#用于染色体的基本展示，有多少染色体就有多少行
data(gene_density, package="RIdeogram")
#包含Chr   Start     End Value 四列
#基因密度文件
# 第一列：染色体号
# 第二列：起始
# 第三列：终止
# 第四列：基因密度值
data(Random_RNAs_500, package="RIdeogram")
#包含Type    Shape Chr    Start      End  color 六列
#染色体旁的标记文件
# 第一列：标记类型
# 第二列：标记形状
# 第三列：染色体号
# 第四列：起始
# 第五列：终止
# 第六列：颜色
#该R包中还提供了一个GFFex函数用来从GFF文件中提取绘制染色体上热图的信息（例如基因密度）
#首先，需要准备物种的karyotype文件，格式与上述的相同，且保证第一列染色体号与GFF文件中的相同。然后，使用GFFex提取基因密度信息
#gene_density <- GFFex(input = "gencode.v32.annotation.gff3", karyotype = "human_karyotype.txt", feature = "gene", window = 1000000)
#其中，feature选项可以改为要绘制的特征类型，window选项可以更改统计的窗口大小
#注意：实测这个函数只是用来计算在限定窗口下染色体区域的gene的数目
#-------------------------------------------------------
#Function:绘图功能
#-------------------------------------------------------
#基本染色体绘制
ideogram(karyotype = human_karyotype)
convertSVG("chromosome.svg", device = "png")
#基因密度热图绘制
ideogram(karyotype = human_karyotype, overlaid = gene_density,output = "chromosome_sample.svg")
convertSVG("chromosome_sample.svg", device = "png",file = "chromosome_sample")
#画前五条染色体
# human_karyotype<-human_karyotype[1:5,]
# ideogram(karyotype = human_karyotype)
# convertSVG("chromosome.svg", device = "png")
# ideogram(karyotype = human_karyotype, overlaid = gene_density)
# convertSVG("chromosome.svg", device = "png")
#标记类型绘制
ideogram(karyotype = human_karyotype, label = Random_RNAs_500, label_type = "marker")
convertSVG("chromosome.svg", device = "png")
#染色体，基因密度和标记同时绘制
ideogram(karyotype = human_karyotype, overlaid = gene_density, label = Random_RNAs_500, label_type = "marker")
convertSVG("chromosome.svg", device = "png")
#以下提到的是RIdeogram的核心函数----ideogram
#ideogram中colorset1可以调整热图色阶，colorset2可以调节label的色阶，width调整染色体宽度，lx和ly用来调试图例的位置
#ideogram图总体由两部分组成，一部分是染色体，另一部分是label，染色体上的数据主要由overlaid参数传递，label的数据主要由label参数传递，由label_type负责展示方式
#比较难理解的是其中的overlaid,label,label_type三个参数
#其中overlaid是在染色体轮廓上进行区域性上色（染色体上的热图）
#labe是在染色体旁边的位置绘上注释
#label_type是决定label的绘图形式，分为marker,heatmap,line,polygon
#注意，对于离散型数据，使用marker，label接受的数据格式可以参考包中自带的数据集data(Random_RNAs_500, package="RIdeogram")
#对于连续性的value数据，可以使用其它形式的label，其中heatmap即在原染色体旁另绘制一根染色体，然后再该染色体上上色
data(human_karyotype, package="RIdeogram") #reload the karyotype data
ideogram(karyotype = human_karyotype, overlaid = gene_density, label = LTR_density, label_type = "heatmap", colorset1 = c("#f7f7f7", "#e34a33"), colorset2 = c("#f7f7f7", "#2c7fb8")) #use the arguments 'colorset1' and 'colorset2' to set the colors for gene and LTR heatmaps, separately.
convertSVG("chromosome.svg", device = "png")
#line是在染色体旁根据value的数值大小以线的波动予以展示，分为单线型和双线型，具体的数据集参考包中自带数据集data(Pi_for_CE, package="RIdeogram") 单线型，data(Pi_for_CE_and_CW, package="RIdeogram")双线型
#polygon是以线下面积的形式展示value值，与line的数据格式一样，只是展示方式不同，具体参见教程
#-------------------------------------------------------
#Function:绘制基因组间共线性
#-------------------------------------------------------
data(synteny_dual_comparison,package = "RIdeogram")
ideogram(karyotype = karyotype_dual_comparison, synteny = synteny_dual_comparison)
convertSVG("chromosome.svg", device = "png")
#其它的详细见简书教程

如果我们需要DIY自己的染色体map，那需要自己准备以下文件1、核型文件，2、gene的位置文件，3、其它与基因相关的注释文件（可选）

#-------------------------------------------------------
#Function:制作用于构建ideogram的基本文件
#-------------------------------------------------------
#1、染色体文件,hg38文件
human_karyotype_hg38<-read.table("E:\\xxx/Homo_sapiens.GRCh38.100.gff3/Homo_sapiens.GRCh38.100.gff3",stringsAsFactors = F, header = F, 
              comment.char = "#", sep = "\t", quote = "")
#gff中第三列中存在chromosome的标签用来标注染色体信息
human_karyotype_hg38<-human_karyotype_hg38[human_karyotype_hg38$V3=="chromosome",]
human_karyotype_hg38<-human_karyotype_hg38[,c(1,4,5)]
colnames(human_karyotype_hg38)<-c("Chr","Start","End")
human_karyotype_hg38<-human_karyotype_hg38[c(order(as.numeric(human_karyotype_hg38$Chr[1:22])),24:25),]
ideogram(karyotype = human_karyotype_hg38)
convertSVG("chromosome.svg", device = "png")
write.table(human_karyotype_hg38,"human_karyotype_hg38.txt",sep = "\t",row.names = F)
#注意：暂时无法获取着丝粒数据，因此染色体上没有凹陷
#发现本包中保存的核型文件是来自genecode库中的hg38版本，与我在ensembl中下载的hg38得到的核型文件的染色体起始位点相同(除了y染色体)
human_karyotype_hg38[,4:5]<-human_karyotype[,4:5]
human_karyotype_hg38[24,4:5]<-0
#2、准备基因信息文件
#gene_density_hg38<-GFFex(input = "E:\\xxx/Homo_sapiens.GRCh38.100.gff3/Homo_sapiens.GRCh38.100.gff3",karyotype="human_karyotype_hg38.txt",feature="gene",window=10000)
#发现这个函数不好用啊，这个函数使用来计算在特定长度的窗口中的基因的个数，，，，，
gene_identity_hg38<-x[x$V3=="gene",]
gene_identity_hg38<-gene_identity_hg38[,c(1,4,5,9)]
gene_identity_hg38$gene<-str_extract(gene_identity_hg38$V9,"Name=[0-9a-zA-Z]{1,};")
gene_identity_hg38$gene<-str_extract(gene_identity_hg38$gene,"=[0-9a-zA-Z]{1,}")
gene_identity_hg38$gene<-str_extract(gene_identity_hg38$gene,"[0-9a-zA-Z]{1,}")
gene_identity_hg38<-gene_identity_hg38[,c(1,2,3,5)]
colnames(gene_identity_hg38)<-c("Chr","Start","End","Gene")
gene_identity_hg38$Value<-sample(1:10,length(gene_identity_hg38$Chr),replace = T)
ideogram(karyotype = human_karyotype_hg38,overlaid = gene_identity_hg38)
convertSVG("chromosome.svg", device = "png")

最后的成图