R包circlize绘制基因组重测序变异圈图示例
作为circos的R替代品，circlize在绘制基因组圈图中同样优秀。本篇随便找了个简单微生物基因组重测序分析中的数据（好吧其实都不算重测序，是瞎搞的，数据本身很有问题，但好在不影响作图......），通过circlize绘制基因组变异圈图展示。数据文件的百度盘链接：
https://pan.baidu.com/s/1QH3_YpNV_JKx2S43YV7lDA
测试数据来源于某重测序分析，这些文件都是在变异检测过程中，软件跑出来的标准文件。
（1）Bacillus_subtilis.str168.gff：参考基因组注释信息gff文件；
（2）Bacillus_subtilis.snp.vcf：SNP检测结果vcf文件；
（3）Bacillus_subtilis.indel.vcf：InDel检测结果vcf文件；
（4）Bacillus_subtilis.cnv：CNV检测结果文件（软件CNVnator）；
（5）Bacillus_subtilis.ctx：SV检测结果文件（软件Breakdancer）；
（6）Bacillus_subtilis.depth_base.txt：参考基因组碱基含量、测序reads覆盖深度统计结果文件，根据滑动窗口统计（滑窗大小使用2000bp）；该文件通过测序reads与参考基因组的比对bam文件统计所得。
接下来在R中处理，基于这些原始文件作统计，包括基因组测序覆盖度、平均GC含量、SNP/InDel密度分布、CNV/SV位点等，将统计后的数据通过circlize绘制基因组变异圈图，如下所示。
重测序变异检测圈图，由外到内依次为：
第一圈，参考物种基因组序列信息，包含序列名称（id）和位置刻度；
第二圈，参考基因组序列的GC含量曲线，这里以2000bp为一滑窗在基因组上滑动，统计的平均GC含量；虚线展示了参考基因组的平均GC含量；
第三圈，测序深度及覆盖度信息，这里以2000bp为一滑窗在基因组上滑动，统计的平均测序深度，反应出不同区域的reads覆盖情况；虚线展示了整体水平的平均reads覆盖度；
第四圈，展示了参考基因组中的基因编码区（CDS）以及非编码RNA区（rRNA、tRNA），以两圈表示，外圈代表正链，内圈代表负链；
第五圈，SNP（单核苷酸突变）密度分布，这里以2000bp为一滑窗在基因组上滑动，统计各区间内的SNP总数量，颜色越深代表SNP越多；
第六圈，InDel（基因组小片段的插入缺失）密度分布，这里以2000bp为一滑窗在基因组上滑动，统计各区间内的InDel总数量，颜色越深代表InDel越多；
第七圈，CNV（基因组大片段的拷贝数变异）信息，绿色表示大片段缺失（Deletion），红色表示大片段重复（Duplication，本示例图中无Duplication结构）；
第八圈，SV（染色体结构变异）信息，绿色表示染色体缺失（DEL），红色表示染色体插入（INS，本示例图中无INS结构）；圈内红色连线代表染色体内部易位（ITX），绿色连线代表染色体间易位（CTX，本示例图中无CTX结构），蓝色连线代表染色体倒位（INV）。
圈图右侧为图例信息，包含了对圈图中各部分内容的注释。
圈图左侧为一些简单的信息，包含了样本名称、参考基因组名称、参考基因组序列总长、参考基因组GC含量，重测序变异检测结果的简要统计，包含碱基替换类型数量、碱基转换/颠换比例、碱基插入缺失类型数量、大片段变异类型数量等。
数据预处理
在正式绘制圈图之前，我们需要对输入文件作一些基本的统计。例如，首先需要统计参考基因组的长度、GC含量、基因位置信息等，以便在基因组圈图中定位。vcf文件中的变异位点（SNP、InDel）是精确到序列中具体的碱基位置的，由于SNP等过多，直接展示所有的点会非常难处理，此处就根据指定滑窗大小按区段（这里默认使用2000bp）统计SNP、InDel的频数（密度）用作绘图，并在同时统计碱基替换类型、插入/缺失类型等，作为图例注释添加在圈图两侧。其它文件亦是如此。
#指定文件
sample_name <- 'Bacillus_subtilis' #测序样本名称
ref_name <- 'Bacillus_subtilis str168' #参考基因组名称
genome_gff <- 'Bacillus_subtilis.str168.gff' #参考基因组 gff 注释文件
snp_vcf <- 'Bacillus_subtilis.snp.vcf' #SNP 检测结果 vcf 文件
indel_vcf <- 'Bacillus_subtilis.indel.vcf' #InDel 检测结果 vcf 文件
cnv_cnv <- 'Bacillus_subtilis.cnv' #cnv 检测结果文件
sv_ctx <- 'Bacillus_subtilis.ctx' #sv 检测结果文件
depth_base_stat <- 'Bacillus_subtilis.depth_base.txt' #测序深度、碱基含量统计结果文件
seq_split <- 2000 #滑窗大小，与 depth_base_stat 中使用的滑窗大小对应
out_dir = 'output' #生成一个目录，用于存放结果文件
if (!file.exists(out_dir)) dir.create(out_dir)
参考基因组范围、GC含量、测序深度统计（第一、二、三圈数据）
读取测序深度、碱基含量统计结果文件“Bacillus_subtilis.depth_base.txt”，据此统计参考基因组总长度、GC含量等信息。同时统计该文件中的平均测序深度，并根据非0值深度区域的长度大致估算测序reads在参考序列中的总覆盖范围。
##参考基因组长度、GC 统计 & 测序覆盖度、深度统计depth_base <- read.delim(depth_base_stat, stringsAsFactors = FALSE)
genome_size <- sum(depth_base$seq_end - depth_base$seq_start + 1)
genome_GC <- round(mean(depth_base$GC), 2)
depth_exist <- subset(depth_base, depth != 0)
coverage <- round(100 * sum(depth_exist$seq_end - depth_exist$seq_start + 1) / genome_size, 2)
average_depth <- round(mean(depth_base$depth), 0)
seq_stat <- NULL
for (seq_id in unique(depth_base$seq_ID)) seq_stat <- rbind(seq_stat, c(seq_id, 1, max(subset(depth_base, seq_ID == seq_id)$seq_end)))
seq_stat <- data.frame(seq_stat, stringsAsFactors = FALSE)
colnames(seq_stat) <- c('seq_ID', 'seq_start', 'seq_end')
rownames(seq_stat) <- seq_stat$seq_ID
seq_stat$seq_start <- as.numeric(seq_stat$seq_start)
seq_stat$seq_end <- as.numeric(seq_stat$seq_end)
write.table(seq_stat, str_c(out_dir, '/', sample_name, '.genome_stat.txt'), row.names = FALSE, sep = '\t', quote = FALSE)
数据框“seq_stat”中记录了参考基因组序列名称，起始和终止位置，在下文绘制圈图时用于确定全局边界范围。
其它的统计结果，用作图例注释。
参考基因组CDS、rRNA、tRNA位置（第四圈数据）
根据参考基因组gff注释文件“Bacillus_subtilis.str168.gff”中的信息，获取参考基因组序列中蛋白编码区（CDS），以及主要的非编码RNA区（rRNA、tRNA）所在的位置信息。
##参考基因组基因信息，CDS & rRNA & tRNA
gene <- read.delim(genome_gff, header = FALSE, stringsAsFactors = FALSE, comment.char = '#')[c(1, 3, 4, 5, 7)]
gene <- subset(gene, V3 %in% c('CDS', 'rRNA', 'tRNA'))[c(1, 3, 4, 2, 5)]
names(gene) <- c('seq_ID', 'seq_start', 'seq_end', 'type', 'strand')
gene[which(gene$type == 'CDS'),'type'] <- 1
gene[which(gene$type == 'rRNA'),'type'] <- 2
gene[which(gene$type == 'tRNA'),'type'] <- 3
gene$type <- as.numeric(gene$type)
gene <- list(subset(gene, strand == '-')[-5], subset(gene, strand == '+')[-5])
读取文件后，根据CDS、rRNA、tRNA所在基因组正链还是负链，生成两个数据框，分别存储基因组正链、负链中的CDS、rRNA、tRNA位置信息，并最终将两个数据框存放在一个list中。
如下所示，最终得到列表“gene”，包含两个数据框，数据框[1]记录了基因组负链中的CDS、rRNA、tRNA位置信息，数据框[2]记录了基因组正链中的CDS、rRNA、tRNA位置信息（由于circlize作图时，对于同一个列表中的不同数据框，默认由内圈往外圈依次绘制，因此若想将正链的圈绘制在外侧，需要将正链数据框的顺序排在列表的后面）。对于每个数据框：seq_ID，参考基因组序列id；seq_start，基因在参考序列中的起始位置；seq_end，基因在参考序列中的结束位置；type，1为CDS，2为rRNA，3为tRNA。
SNP频数（密度）统计（第五圈数据）
读取SNP检测结果vcf文件“Bacillus_subtilis.snp.vcf”，选择以滑动窗口的形式，统计基因组序列各区间内的SNP类型和频数并在圈图中展示。
#读取 vcf 文件，统计 SNP 类型
snp <- read.delim(snp_vcf, header = FALSE, colClasses = 'character', comment.char = '#')[c(1, 2, 4, 5)]
snp$V2 <- as.numeric(snp$V2)
snp$change <- str_c(snp$V4, snp$V5)
change <- which(snp$change == 'AT')
snp[change,'type1'] <- 'A>T|T>A'; snp[change,'type2'] <- 'tv'
change <- which(snp$change == 'AG')
snp[change,'type1'] <- 'A>G|T>C'; snp[change,'type2'] <- 'ti'
change <- which(snp$change == 'AC')
snp[change,'type1'] <- 'A>C|T>G'; snp[change,'type2'] <- 'tv'
change <- which(snp$change == 'TA')
snp[change,'type1'] <- 'A>T|T>A'; snp[change,'type2'] <- 'tv'
change <- which(snp$change == 'TG')
snp[change,'type1'] <- 'A>C|T>G'; snp[change,'type2'] <- 'tv'
change <- which(snp$change == 'TC')
snp[change,'type1'] <- 'A>G|T>C'; snp[change,'type2'] <- 'ti'
change <- which(snp$change == 'GA')
snp[change,'type1'] <- 'G>A|C>T'; snp[change,'type2'] <- 'ti'
change <- which(snp$change == 'GT')
snp[change,'type1'] <- 'G>T|C>A'; snp[change,'type2'] <- 'tv'
change <- which(snp$change == 'GC')
snp[change,'type1'] <- 'G>C|C>G'; snp[change,'type2'] <- 'tv'
change <- which(snp$change == 'CA')
snp[change,'type1'] <- 'G>T|C>A'; snp[change,'type2'] <- 'tv'
change <- which(snp$change == 'CT')
snp[change,'type1'] <- 'G>A|C>T'; snp[change,'type2'] <- 'ti'
change <- which(snp$change == 'CG')
snp[change,'type1'] <- 'G>C|C>G'; snp[change,'type2'] <- 'tv'
snp_ti <- length(which(snp$type2 == 'ti'))
snp_tv <- length(which(snp$type2 == 'tv'))
snp_at <- length(which(snp$type1 == 'A>T|T>A'))
snp_ag <- length(which(snp$type1 == 'A>G|T>C'))
snp_ac <- length(which(snp$type1 == 'A>C|T>G'))
snp_ga <- length(which(snp$type1 == 'G>A|C>T'))
snp_gt <- length(which(snp$type1 == 'G>T|C>A'))
snp_gc <- length(which(snp$type1 == 'G>C|C>G'))
读取vcf文件中的主要信息后，并在每行标记SNP类型，获得统计结果数据框“snp”。第一列（V1），参考基因组序列id；第二列（V2），SNP位点对应的参考序列中的位置；第三列（V4），参考序列该位置处的碱基类型；第四列（V5），SNP替换后的碱基类型；第五列（change），SNP替换类型；第六列（type1），SNP替换类型，合并正负链信息；第七列（type2），表明碱基转换（transitions）还是颠换（transversions）。
SNP的主要类型数量统计结果，用作图例注释。
随后根据预设滑窗大小（本文设置2000bp），统计参考基因组序列中各滑窗区段中所对应的SNP频数（密度）。
#统计 SNP 密度snp <- snp[c(1, 2, 5, 6, 7)]
colnames(snp)[1:2] <- c('seq_ID', 'seq_site')
snp_stat <- NULL
seq_ID <- unique(snp$seq_ID)
for (seq_ID_n in seq_ID) {
snp_subset <- subset(snp, seq_ID == seq_ID_n)
seq_end <- seq_split
snp_num <- 0
for (i in 1:nrow(snp_subset)) {
if (snp_subset[i,'seq_site'] <= seq_end) snp_num <- snp_num + 1
else {
snp_stat <- rbind(snp_stat, c(seq_ID_n, seq_end - seq_split + 1, seq_end, snp_num))
seq_end <- seq_end + seq_split
snp_num <- 0
while (snp_subset[i,'seq_site'] > seq_end) {
snp_stat <- rbind(snp_stat, c(seq_ID_n, seq_end - seq_split + 1, seq_end, snp_num))
seq_end <- seq_end + seq_split
}
snp_num <- snp_num + 1
}
}
while (seq_end < seq_stat[seq_ID_n,'seq_end']) {
snp_stat <- rbind(snp_stat, c(seq_ID_n, seq_end - seq_split + 1, seq_end, snp_num))
seq_end <- seq_end + seq_split
snp_num <- 0
}
snp_stat <- rbind(snp_stat, c(seq_ID_n, seq_end - seq_split + 1, seq_stat[seq_ID_n,'seq_end'], snp_num))
}
snp_stat <- data.frame(snp_stat, stringsAsFactors = FALSE)
names(snp_stat) <- c('seq_ID', 'seq_start', 'seq_end', 'snp_num')
snp_stat$seq_start <- as.numeric(snp_stat$seq_start)
snp_stat$seq_end <- as.numeric(snp_stat$seq_end)
snp_stat$snp_num <- as.numeric(snp_stat$snp_num)
write.table(snp_stat, str_c(out_dir, '/', sample_name, '.snp_stat.txt'), row.names = FALSE, sep = '\t', quote = FALSE)
得到的数据框“snp_stat”中，统计了参考基因组序列中各滑窗区段中所对应的SNP频数（密度）信息。seq_ID，参考基因组序列id；seq_start，各滑窗在参考序列中的起始位置；seq_end，各滑窗在参考序列中的结束位置；snp_num，各滑窗内的SNP数量总计。
InDel频数（密度）统计（第六圈数据）
InDel的统计方法同上述SNP。
首先将“Bacillus_subtilis.indel.vcf”以表格形式读入，并首先在整体水平统计InDel类型和数量信息。
#读取 vcf 文件，统计 InDel 长度
indel <- read.delim(indel_vcf, header = FALSE, colClasses = 'character', comment.char = '#')[c(1, 2, 4, 5)]
indel$V2 <- as.numeric(indel$V2)
indel$length <- str_length(indel[ ,4]) - str_length(indel[ ,3])
indel_insert <- length(which(indel$length > 0))
indel_delet <- length(which(indel$length < 0))
读取vcf文件中的主要信息后，并在每行标记InDel类型，获得统计结果数据框“indel”。第一列（V1），参考基因组序列id；第二列（V2），InDel位点对应的参考序列中的位置；第三列（V4），参考序列该位置处的碱基类型；第四列（V5），InDel后的碱基类型；第五列（length），InDel长度，正值为插入，负值为缺失，数值绝对值代表插入/缺失的碱基长度。
InDel的主要类型数量统计结果，用作图例注释。
随后根据预设滑窗大小（本文设置2000bp），统计参考基因组序列中各滑窗区段中所对应的InDel频数（密度）。
#统计 InDel 密度indel <- indel[c(1, 2, 5)]
colnames(indel)[1:2] <- c('seq_ID', 'seq_site')
indel_stat <- NULL
seq_ID <- unique(indel$seq_ID)
for (seq_ID_n in seq_ID) {
indel_subset <- subset(indel, seq_ID == seq_ID_n)
seq_end <- seq_split
indel_num <- 0
for (i in 1:nrow(indel_subset)) {
if (indel_subset[i,'seq_site'] <= seq_end) indel_num <- indel_num + 1
else {
indel_stat <- rbind(indel_stat, c(seq_ID_n, seq_end - seq_split + 1, seq_end, indel_num))
seq_end <- seq_end + seq_split
indel_num <- 0
while (indel_subset[i,'seq_site'] > seq_end) {
indel_stat <- rbind(indel_stat, c(seq_ID_n, seq_end - seq_split + 1, seq_end, indel_num))
seq_end <- seq_end + seq_split
}
indel_num <- indel_num + 1
}
}
while (seq_end < seq_stat[seq_ID_n,'seq_end']) {
indel_stat <- rbind(indel_stat, c(seq_ID_n, seq_end - seq_split + 1, seq_end, indel_num))
seq_end <- seq_end + seq_split
indel_num <- 0
}
indel_stat <- rbind(indel_stat, c(seq_ID_n, seq_end - seq_split + 1, seq_stat[seq_ID_n,'seq_end'], indel_num))
}
indel_stat <- data.frame(indel_stat, stringsAsFactors = FALSE)
names(indel_stat) <- c('seq_ID', 'seq_start', 'seq_end', 'indel_num')
indel_stat$seq_start <- as.numeric(indel_stat$seq_start)
indel_stat$seq_end <- as.numeric(indel_stat$seq_end)
indel_stat$indel_num <- as.numeric(indel_stat$indel_num)
write.table(indel_stat, str_c(out_dir, '/', sample_name, '.indel_stat.txt'), row.names = FALSE, sep = '\t', quote = FALSE)
得到的数据框“indel_stat”中，统计了参考基因组序列中各滑窗区段中所对应的InDel频数（密度）信息。seq_ID，参考基因组序列id；seq_start，各滑窗在参考序列中的起始位置；seq_end，各滑窗在参考序列中的结束位置；indel_num，各滑窗内的InDel数量总计。
CNV统计（第七圈数据）
读取CNV变异检测结果文件“Bacillus_subtilis.cnv”，获取CNV位置、类型和数量。
##读取 CNV 检测结果
cnv <- read.delim(cnv_cnv, header = FALSE, stringsAsFactors = FALSE)[1:2]
for (i in 1:nrow(cnv)) {
str_split <- unlist(str_split(cnv[i,2], '[:-]'))
cnv[i,'seq_ID'] <- str_split[1]
cnv[i,'seq_start'] <- str_split[2]
cnv[i,'seq_end'] <- str_split[3]
}
cnv <- cnv[c(3:5, 1)]
cnv$seq_start <- as.numeric(cnv$seq_start)
cnv$seq_end <- as.numeric(cnv$seq_end)
cnv_dup <- length(which(cnv[[4]] == 'duplication'))
cnv_del <- length(which(cnv[[4]] == 'deletion'))
读取CNV变异检测结果文件“Bacillus_subtilis.cnv”中的CNV位置和类型信息，得到数据框“cnv”。第一列（seq_ID），参考基因组序列id；第二列（seq_start），CNV在参考序列中的起始位置；第三列（seq_end），CNV在参考序列中的结束位置；第四列（V1），CNV类型，duplication为大片段重复，deletion为大片段缺失。
CNV的主要类型数量统计结果，用作图例注释。
SV统计（第八圈数据）
最后统计染色体结构变异（SV）类型数量。
##读取 SV 检测结果sv <- read.delim(sv_ctx, header = FALSE, stringsAsFactors = FALSE, comment.char = '#')[c(1, 2, 4, 5, 7)]
sv_in_de <- subset(sv, sv[[5]] %in% c('INS', 'DEL'))[c(1, 2, 4, 5)]
sv_ITX <- subset(sv, sv[[5]] %in% 'ITX')
sv_CTX <- subset(sv, sv[[5]] %in% 'CTX')
sv_INV <- subset(sv, sv[[5]] %in% 'INV')
sv_ins <- length(which(sv[[5]] == 'INS'))
sv_del <- length(which(sv[[5]] == 'DEL'))
sv_itx <- nrow(sv_ITX)
sv_ctx <- nrow(sv_CTX)
sv_inv <- nrow(sv_INV)
读取SV变异检测结果文件“Bacillus_subtilis.ctx”中的SV位置和类型信息，得到数据框“sv”。并根据SV类型，分割为不同的子数据集。
对于SV插入/缺失（数据框“sv_in_de”）：第一列（V1），参考基因组序列id；第二列（V2），SV在参考序列中的起始位置；第三列（V5），SV在参考序列中的结束位置；第四列（V7），SV类型，DEL为缺失（deletion），INS为插入（insertion）。
对于SV倒位和易位（数据框“sv_ITX”、“sv_CTX”、“sv_INV”）：第一列（V1），SV倒位/易位前所位于的参考基因组序列id；第二列（V2），SV倒位/易位前的位置；第三列（V4），SV倒位/易位后所位于的参考基因组序列id；第四列（V5），SV倒位/易位后的位置；第五列（V7），SV类型，ITX为染色体内部的易位（intra-chromosomal translocation），CTX为染色体之间的易位（inter-chromosomal translocation），INV为倒位（inversion）。
其它统计信息，用作图例注释。
基因组重测序变异圈图绘制
然后接下来就是使用circlize包绘制基因组圈图了。
注：绘图细节部分的调整很多也很繁琐，此处默认使用的各细节参数设置（如字体大小等）均根据示例文件而来。若是大家后续更换为自己的数据时，还需多加调试参数了。
基因组全局（第一圈）
首先定义整体的布局。
library(circlize) #绘制圈图library(grid) #设置画板
pdf(str_c(out_dir, '/', sample_name, '.circlize.pdf'), width = 14, height = 8)
circle_size = unit(1, 'snpc')
circos.par(gap.degree = 2, start.degree = 90)
#第一圈，染色体区域
circos.genomicInitialize(seq_stat, plotType = c('axis', 'labels'), major.by = 250000, track.height = 0.05)
circos.genomicTrackPlotRegion(
seq_stat, track.height = 0.05, stack = TRUE, bg.border = NA,
panel.fun = function(region, value, ...) {
circos.genomicRect(region, value, col = '#049a0b', border = NA, ...)
} )
预先生成pdf文件，将图片绘制到pdf中。我这里设置了一个14×8大小的pdf。
circos.genomicInitialize()首先根据参考基因组的长度统计结果，确定了基因组边界范围，以及绘制刻度线。circos.track()绘制参考基因组序列区块，并添加序列id。
此时若使用“dev.off()”关闭pdf画板，则结果如下。
参考基因组GC含量曲线（第二圈）
基因组范围确定后，添加GC含量信息。在前面我们已经按滑窗统计了参考基因组中的平均GC含量信息，继续添加至圈图中。
#第二圈，GC% 含量图
circos.genomicTrack(
depth_base[c(1:3, 5)], track.height = 0.08, bg.col = '#EEEEEE6E', bg.border = NA,
panel.fun = function(region, value, ...) {
circos.genomicLines(region, value, col = 'blue', lwd = 0.35, ...)
circos.lines(c(0, max(region)), c(genome_GC, genome_GC), col = 'blue2', lwd = 0.15, lty = 2)
circos.yaxis(labels.cex = 0.2, lwd = 0.1, tick.length = convert_x(0.15, 'mm'))
} )
使用circos.genomicTrack()，读取数据框“depth_base”中第1（序列id）、2（滑窗起始位置）、3（滑窗结束位置）、5（GC含量百分比）列的结果将GC含量信息添加在图中。track.height控制圈图区域高度，bg.col设置背景色浅灰色，bg.border设置边框颜色无。定义的函数panel.fun()使用循环的方式，依次将各滑窗的位置信息、GC含量信息添加在圈图中；circos.genomicLines()绘制GC曲线，col定义颜色蓝色，lwd设置曲线宽度；circos.lines()在途中额外添加一条虚线弧线，即基因组的平均GC含量，lty = 2即绘制为虚线；circos.yaxis()用于在边界处添加标尺，内部各参数可控制标尺大小、线长等。
此时若使用“dev.off()”关闭pdf画板，则结果如下（展示了部分结果）。
测序覆盖度、深度柱形图（第三圈）
添加测序覆盖度和深度信息，用于反应本次测序数据量，以及各变异检测结果是否具有合理的依据。在前面我们已经按滑窗统计了reads的平均覆盖深度，继续添加至圈图中。
#第三圈，覆盖度 & 深度
circos.genomicTrack(
depth_base[1:4], track.height = 0.08, ylim = c(0, (max(depth_base$depth) + 1)), bg.col = '#EEEEEE6E', bg.border = NA,
panel.fun = function(region,value, ...) {
circos.genomicRect(region, value, ytop.column = 1, ybottom = 0, border = 'white', lwd = 0.02, col = 'red', ...)
circos.lines(c(0, max(region)), c(average_depth, average_depth), col = 'red3', lwd = 0.15, lty = 2)
circos.yaxis(labels.cex = 0.2, lwd = 0.1, tick.length = convert_x(0.15, 'mm'))
} )
使用circos.genomicTrack()，读取数据框“depth_base”中第1（序列id）、2（滑窗起始位置）、3（滑窗结束位置）、4（reads覆盖深度信息）列的结果将reads覆盖深度信息添加在图中。track.height控制圈图区域高度，ylim确定了最大深度，bg.col设置背景色浅灰色，bg.border设置边框颜色无。定义的函数panel.fun()使用循环的方式，依次将各滑窗的reads覆盖深度信息添加在圈图中；circos.genomicRect()绘制柱形图，col定义柱形图颜色红色；其余主要参数和上述绘制GC含量曲线时的类似，circos.lines()在途中额外添加一条虚线弧线，即整体水平的平均测序深度，lty = 2即绘制为虚线；circos.yaxis()用于在边界处添加标尺，内部各参数可控制标尺大小、线长等。
此时若使用“dev.off()”关闭pdf画板，则结果如下（展示了部分结果）。
参考基因组CDS、rRNA、tRNA区域（第四圈）
将参考基因组序列中蛋白编码区（CDS），以及主要的非编码RNA区（rRNA、tRNA）添加在圈图中，以用于反应各变异位置所对应的基因水平。
#第四圈，CDS & rRNA & tRNA
color_assign <- colorRamp2(breaks = c(1, 2, 3), col = c('#00ADFF', 'orange', 'green2'))
circos.genomicTrackPlotRegion(
gene, track.height = 0.12, stack = TRUE, bg.border = NA,
panel.fun = function(region, value, ...) {
circos.genomicRect(region, value, col = color_assign(value[[1]]), border = NA, ...)
} )
colorRamp2()使用了3种颜色，分别表示CDS、rRNA、tRNA。
使用circos.genomicTrackPlotRegion()，读取列表“gene”中的两个数据框（分别记录了基因组正链、负链中的基因信息），并在圈图中绘制各基因的所处位置。track.height控制圈图区域高度，bg.border设置边框颜色无。定义的函数panel.fun()使用循环的方式，依次将各CDS、rRNA、tRNA区域的位置信息添加在圈图中；其中，正链绘制在外圈，负链绘制在内圈（可以看到，circlize作图时，对于同一个列表中的不同数据框，默认由内圈往外圈依次绘制，因此若想将正链的圈绘制在外侧，需要将正链数据框的顺序排在列表的后面）；CDS、rRNA、tRNA分别以3种颜色表示。
此时若使用“dev.off()”关闭pdf画板，则结果如下（展示了部分结果）。
SNP、InDel频数密度圈（第五、六圈）
添加重测序变异检测信息。首先我们将上文中基于SNP、InDel检测结果vcf文件，按滑窗（2000bp）统计得到的SNP、InDel频数密度信息添加在圈图中。
#第五圈，SNP 密度
value_max <- max(snp_stat$snp_num)
colorsChoice <- colorRampPalette(c('white', '#245B8E'))
color_assign <- colorRamp2(breaks = c(0:value_max), col = colorsChoice(value_max + 1))
circos.genomicTrackPlotRegion(
snp_stat, track.height = 0.08, stack = TRUE, bg.border = NA,
panel.fun = function(region, value, ...) {
circos.genomicRect(region, value, col = color_assign(value[[1]]), border = NA, ...)
} )
#第六圈，InDel 密度
value_max <- max(indel_stat$indel_num)
colorsChoice <- colorRampPalette(c('white', '#7744A4'))
color_assign <- colorRamp2(breaks = c(0:value_max), col = colorsChoice(value_max + 1))
circos.genomicTrackPlotRegion(
indel_stat, track.height = 0.08, stack = TRUE, bg.border = NA,
panel.fun = function(region, value, ...) {
circos.genomicRect(region, value, col = color_assign(value[[1]]), border = NA, ...)
} )
SNP和InDel密度信息的绘制方法一致，这里合并在一起说明了。
首先分别统计所有滑窗中，SNP、InDel密度的最大值；并使用colorRampPalette()，设定一个颜色范围；最后使用colorRamp2()，将预设颜色范围的最大值规定为SNP、InDel密度的最大值，颜色范围的最小值规定为0，以便以渐变色的形式表示展示SNP、InDel密度的大小。
之后使用circos.genomicTrackPlotRegion()，根据滑窗统计结果“snp_stat”和“indel_stat”中的信息，在圈图中绘制绘制出对应的滑窗位置。track.height控制圈图区域高度，bg.border设置边框颜色无。定义的函数panel.fun()使用循环的方式，依次将各滑窗添加在圈图中；circos.genomicRect()用于绘制滑窗区块，参数col = color_assign()根据滑窗统计结果中的密度数值，识别渐变色范围给滑窗上色（根据SNP、InDel密度不同数值，自动匹配渐变色范围）。
此时若使用“dev.off()”关闭pdf画板，则结果如下（展示了部分结果）。
CNV、SV变异圈（第七、八圈）
添加CNV和SV信息。上文已经将CNV和SV的位置信息读取到R中，并标记了变异类型，使用读取的数据继续做图即可。
#第七圈，CNV 变异；duplication、deletion，1 插入、2 缺失
for (i in 1:nrow(cnv)) {
if (cnv[i,4] == 'duplication') cnv[i,4] <- 1
if (cnv[i,4] == 'deletion') cnv[i,4] <- 2
}
cnv[[4]] <- as.numeric(cnv[[4]])
color_assign <- colorRamp2(breaks = c(1, 2), color = c('#FB564A', '#74C476'))
circos.genomicTrackPlotRegion(
cnv, track.height = 0.06, stack = TRUE, bg.border = NA,
panel.fun = function(region, value, ...) {
circos.genomicRect(region, value, col = color_assign(value[[1]]), border = NA, ...)
} )
#第八圈，SV 变异；INS、DEL，1 插入、2 缺失
for (i in 1:nrow(sv_in_de)) {
if (sv_in_de[i,4] == 'INS') sv_in_de[i,4] <- 1
if (sv_in_de[i,4] == 'DEL') sv_in_de[i,4] <- 2
}
sv_in_de[[4]] <- as.numeric(sv_in_de[[4]])
color_assign <- colorRamp2(breaks = c(1, 2), color = c('#FB564A', '#74C476'))
circos.genomicTrackPlotRegion(
sv_in_de, track.height = 0.06, stack = TRUE, bg.border = NA,
panel.fun = function(region, value, ...) {
circos.genomicRect(region, value, col = color_assign(value[[1]]), border = NA, ...)
} )
#内部连线分别为
#ITX，内易位
circos.genomicLink(sv_ITX[c(1, 2, 2)], sv_ITX[c(3, 4, 4)], col = '#FF46346E', lwd = 0.5)
#CTX，间易位
circos.genomicLink(sv_CTX[c(1, 2, 2)], sv_CTX[c(3, 4, 4)], col = '#66CB716E', lwd = 0.5)
#INV，倒位
circos.genomicLink(sv_INV[c(1, 2, 2)], sv_INV[c(3, 4, 4)], col = '#53C6FF6E', lwd = 0.5)
对于CNV的插入和缺失类型以及SV的插入和缺失类型，二者的绘制方法一致。首先分别在上文所得的统计结果数据框“cnv”和“sv_in_de”中识别二者的变异类型（1，插入；2，缺失），分别以不同的颜色表示。然后使用circos.genomicTrackPlotRegion()读取变异的位置信息，并添加在圈图中。定义的函数panel.fun()使用循环的方式，依次添加变异位置；circos.genomicRect()用于绘制插入/缺失的变异区块，col = color_assign参数为不同变异类型指定颜色。
对于SV的ITX、CTX、INV类型，以连线的形式表示，绘制在圈图的最内侧。circos.genomicLink()绘制圈图连线，前两个参数即指定连线起始和终止的位置，col设置连线颜色，lwd设置连线宽度。这里使用的数据框“sv_ITX”、“sv_CTX”、“sv_INV”均为上文所得，记录了连线位置（SV变异前后位置）信息。
此时若使用“dev.off()”关闭pdf画板，则结果如下（展示了部分结果）。
添加图例
圈图主体部分，到了上一步就已经完成了。
circlize绘制的基因组圈图是不包含图例的，在circlize文档中，提到可通过ComplexHeatmap包Legend()函数，为circlize图添加图例。
我们继续使用Legend()函数为圈图添加图例。同时使用grid包中的函数确定图例位置，并最终和圈图组合在一起。
library(ComplexHeatmap) #绘制图例#GC 含量曲线图例
gc_legend <- Legend(
at = 1, labels = c(str_c('GC % ( Average: ', genome_GC, ' % )')), labels_gp = gpar(fontsize = 8),
grid_height = unit(0.5, 'cm'), grid_width = unit(0.5, 'cm'), type = 'lines', background = '#EEEEEE6E',
legend_gp = gpar(col = 'blue', lwd = 0.5))
#测序深度、覆盖度图例
depth_legend <- Legend(
at = 1, labels = str_c(' Depth ( average: ', average_depth, ' X )'), labels_gp = gpar(fontsize = 8),
title = str_c('Coverage: ', coverage , ' %'), title_gp = gpar(fontsize = 9),
grid_height = unit(0.4, 'cm'), grid_width = unit(0.4, 'cm'), type = 'points', pch = NA, background = 'red')
#CDS & rRNA & tRNA 图例
gene_legend <- Legend(
at = c(3, 2, 1), labels = c(' CDS', ' rRNA', ' tRNA'), labels_gp = gpar(fontsize = 8),
title = 'CDS | rRNA | tRNA', title_gp = gpar(fontsize = 9),
grid_height = unit(0.4, 'cm'), grid_width = unit(0.4, 'cm'), type = 'points', pch = NA, background = c('#00ADFF', 'orange', 'green2'))
#SNP 密度图例
snp_legend <- Legend(
at = round(seq(0, max(snp_stat$snp_num), length.out = 6), 0), labels_gp = gpar(fontsize = 8),
col_fun = colorRamp2(round(seq(0, max(snp_stat$snp_num), length.out = 6), 0), colorRampPalette(c('white', '#245B8E'))(6)),
title_position = 'topleft', title = 'SNP density', legend_height = unit(4, 'cm'), title_gp = gpar(fontsize = 9))
#InDel 密度图例
indel_legend <- Legend(
at = round(seq(0, max(indel_stat$indel_num), length.out = 6), 0), labels_gp = gpar(fontsize = 8),
col_fun = colorRamp2(round(seq(0, max(indel_stat$indel_num), length.out = 6), 0), colorRampPalette(c('white', '#7744A4'))(6)),
title_position = 'topleft', title = 'InDel density', legend_height = unit(4, 'cm'), title_gp = gpar(fontsize = 9))
#CNV 图例
cnv_legend <- Legend(
at = c(1, 2), labels = c(' Duplication', ' Deletion'), labels_gp = gpar(fontsize = 8), title = 'CNV type', title_gp = gpar(fontsize = 9),
grid_height = unit(0.5, 'cm'), grid_width = unit(0.5, 'cm'), type = 'points', pch = NA, background = c('#FB564A', '#74C476'))
#SV 图例
sv_legend <- Legend(
at = c(1, 2, 3, 4, 5), labels = c(' INS', ' DEL', ' ITX',' CTX', ' INV'),
labels_gp = gpar(fontsize = 8), grid_height = unit(0.5, 'cm'), grid_width = unit(0.5, 'cm'),
type = c('points', 'points', 'lines', 'lines', 'lines'), pch = NA, background = c('#FB564A', '#74C476', NA, NA, NA),
legend_gp = gpar(col = c(NA, NA, '#FF4634', '#66CB71', '#53C6FF'), lwd = 1),
title = 'SV type', title_gp = gpar(fontsize = 9))
#左侧统计总览（涵括了上文大部分的统计概况，例如 SNP 替换类型统计等）
stat_legend <- Legend(
at = 1, labels = '1', labels_gp = gpar(fontsize = 0), title_gp = gpar(fontsize = 9),
grid_height = unit(0, 'cm'), grid_width = unit(0, 'cm'), type = 'points', pch = NA, background = NA,
title = str_c('Sample: ', sample_name, '\nRefer species: ', ref_name, '\nRefer size: ', genome_size, ' bp\nRefer GC: ', genome_GC, ' %\n\n\nTotal SNP: ', snp_ti + snp_tv, '\nTransitions: ', snp_ti, '\nTransversions: ', snp_tv, '\nTi/Tv: ', round(snp_ti / snp_tv, 2), '\nA>T|T>A: ', snp_at, '\nA>G|T>C: ', snp_ag, '\nA>C|T>G: ', snp_ac, '\nG>A|C>T: ', snp_ga, '\nG>T|C>A: ', snp_gt, '\nG>C|C>G: ', snp_gc, '\n\n\nTotal InDel: ', indel_insert + indel_delet, '\nInsert: ', indel_insert, '\nDelet: ', indel_delet, '\n\n\nTotal CNV: ', cnv_dup + cnv_del, '\nDuplication: ', cnv_dup, '\nDeletion: ', cnv_del, '\n\n\nTotal SV: ', sv_ins + sv_del + sv_itx + sv_ctx + sv_inv, '\nINS: ', sv_ins,'\nDEL: ', sv_del, '\nITX: ', sv_itx, '\nCTX: ', sv_ctx, '\nINV: ', sv_inv))
#最后添加图例，图例在图中的存放位置自己看着调吧
y_coord <- 0.8
x_coord <- 0.87
pushViewport(viewport(x = x_coord + 0.011, y = y_coord))
grid.draw(gc_legend)
y_coord <- y_coord - 0.06
upViewport()
pushViewport(viewport(x = x_coord + 0.005, y = y_coord))
grid.draw(depth_legend)
y_coord <- y_coord - 0.1
upViewport()
pushViewport(viewport(x = x_coord - 0.0063, y = y_coord))
grid.draw(gene_legend)
y_coord <- y_coord - 0.19
upViewport()
pushViewport(viewport(x = x_coord - 0.0205, y = y_coord))
grid.draw(snp_legend)
y_coord <- y_coord
upViewport()
pushViewport(viewport(x = x_coord + 0.0505, y = y_coord))
grid.draw(indel_legend)
y_coord <- y_coord - 0.195
upViewport()
pushViewport(viewport(x = x_coord - 0.018, y = y_coord))
grid.draw(cnv_legend)
y_coord <- y_coord
upViewport()
pushViewport(viewport(x = x_coord + 0.04, y = y_coord - 0.0373))
grid.draw(sv_legend)
y_coord <- y_coord
upViewport()
pushViewport(viewport(x = 0.12, y = 0.5))
grid.draw(stat_legend)
upViewport()
#最后清除预设样式并关闭画板，以免影响继续作图
circos.clear()
dev.off()
图例添加完毕后，dev.off()关闭pdf画板后，最终效果图如下。
R包circlize绘制弦状图示例（多组分信息样式，从头绘制）