karyoploteR-染色体全局可视化
六月份的一个学徒自己有甲基化项目数据要处理,所以开始自学一些R包了
rm(list = ls())
options(stringsAsFactors = F)
#BiocManager::install("karyoploteR")
library(karyoploteR)
输入格式
先建立data frame 之后处理成toGRanges格式 这里先绘制常用的hg19参考基因组(plotKaryotype,人类) 再用颜色标志出特定的位点(kpPlotRegions)
a <- data.frame(chr=c("chr1", "chr5", "chr17", "chr22"), start=c(1, 1000000, 4000000, 1),end=c(5000000, 3200000, 80000000, 1200000))
head(a)
gains <- toGRanges(a)
class(gains)
losses <- toGRanges(data.frame(chr=c("chr3", "chr9", "chr17"), start=c(80000000, 20000000, 1),end=c(170000000, 30000000, 25000000)))
# 需要一起run
kp <- plotKaryotype(genome="hg19")
kpPlotRegions(kp, gains, col="#FFAACC")
三个参数
plotKaryotype 常用有三个参数:genome, chromosomes and plot.type genome 参考基因组(默认hg19) chromosomes 可展示特定染色体(默认"autosomal") plot.type 位点展示的位置(above, below or on the ideograms)
# chromosomes=c("chr1", "chr2", "chr3")
kp <- plotKaryotype(genome="hg19", plot.type = 1, main="The mm10 genome")
kpPlotRegions(kp, gains, col="#FFAACC")
plot.type
plot.type=1 A plot with stacked horizontal ideograms and a single data panel on top of them plot.type=2 A plot with stacked horizontal ideograms and two data panels, one above and one below them plot.type=3 A plot with horizontal ideograms all in one level, and two data panels, one above and one below them plot.type=4 A plot with horizontal ideograms all in one level, and one data panel above them plot.type=5 A plot with horizontal ideograms all in one level, and one data panel below them 本质上调节参考基因组和数据面板位置 kpDataBackground 模拟背景,更好的理解参数
kp <- plotKaryotype(chromosomes = c("chr1", "chr2"), plot.type = 1)
kpDataBackground(kp, data.panel = 1)
kp <- plotKaryotype(chromosomes = c("chr1", "chr2"), plot.type = 2)
kpDataBackground(kp, data.panel = 1)
kpDataBackground(kp, data.panel = 2)
plot.type=2 上下有个数据面板
kp <- plotKaryotype(chromosomes=c("chr1", "chr2"), plot.type=2)
# 分开设置,首先针对上面的数据面板
kpDataBackground(kp,data.panel = 1) # 设置背景
# 设置左端y轴
kpAxis(kp)
# 添加右端的y轴
kpAxis(kp,side = 2)
# 设置下面的数据面板
# 针对data.panel=2设置了两次,下面的数据面板也分隔成两份
# r1,r0控制比例
kpDataBackground(kp, r1=0.47, data.panel=2)
#Changing the limits and having more ticks, with a smaller font size
kpAxis(kp, r1=0.47, ymin=-5000, ymax = 5000, numticks = 5, cex=0.5, data.panel=2)
#and a different scale on the right
kpAxis(kp, r1=0.47, ymin=-2, ymax = 2, numticks = 3, cex=0.5, data.panel=2, side=2)
kpDataBackground(kp, r0=0.53, data.panel=2)
#Changing the colors and labels and tick positions
kpAxis(kp, r0=0.53, tick.pos = c(0.3, 0.6, 1), labels = c("A", "B", "C"), col="#66AADD",cex=0.5, data.panel=2) # 只设置了左边
添加图形
添加点kpPoints和方形kpRect
kp <- plotKaryotype(chromosomes=c("chr1"))
kpDataBackground(kp)
kpAxis(kp)
kpPoints(kp, chr="chr1", x=30000000, y=0.2)
kpRect(kp, chr="chr1", x0=100000000, x1=120000000, y0=0.2, y1=0.4)
一共10种 kpAbline Draws horizontal and vertical lines spanning the whole avaliable space kpArrows Draws arrows. Using the code parameter, it’s possible to specify where the arrowhead should be placed. kpBars Draws vertical bars. If y0 is present, bars span from y0 to y1, if it’s ommited, bars span from 0 to y1 kpHeatmap Draws a heatmap-like representation of the data. Specifically, for each data point a rectangle is drawn with the color determined by the y value kpLines Draws straight lines joining the data points. kpPoints Draws points (or other shaes using pch) in the data points. kpPolygon Draws a closed polygon joining all data points kpRect Draws a rectangle at the specified points. kpSegments Draws segments joining the specified points kpText Draws text labels at the specified positions. Text labels are provided via the labels argument.
pp <- getDefaultPlotParams(plot.type = 1)
pp$data1height=600
tr.i <- 1/11
tr.o <- 1/10
kp <- plotKaryotype(chromosomes=c("chr1"), plot.params = pp)
dd <- toGRanges(data.frame(chr="chr1", start=end(kp$genome[1])/50*(0:49), end=end(kp$genome[1])/50*(1:50)))
mcols(dd) <- data.frame(y=((sin(start(dd)) + rnorm(n=50, mean=0, sd=0.1))/5)+0.5)
tn <- 0
kpDataBackground(kp, r0=tr.o*tn, r1=tr.o*tn+tr.i)
kpPoints(kp, dd, r0=tr.o*tn, r1=tr.o*tn+tr.i, pch=".", cex=2)
kpRect(kp, chr="chr1", x0=5000000, x1=45000000, y0=0.2, y1=0.8, r0=tr.o*tn, r1=tr.o*tn+tr.i, col="#EEEEEE", border="#666666")
kpText(kp, chr="chr1", x=25000000, y=0.5, col="red", r0=tr.o*tn, r1=tr.o*tn+tr.i, labels="kpPoints", cex=0.7)
tn <- 1
kpDataBackground(kp, r0=tr.o*tn, r1=tr.o*tn+tr.i)
kpLines(kp, dd, r0=tr.o*tn, r1=tr.o*tn+tr.i, pch=".", cex=2)
kpRect(kp, chr="chr1", x0=5000000, x1=45000000, y0=0.2, y1=0.8, r0=tr.o*tn, r1=tr.o*tn+tr.i, col="#EEEEEE", border="#666666")
kpText(kp, chr="chr1", x=25000000, y=0.5, col="red", r0=tr.o*tn, r1=tr.o*tn+tr.i, labels="kpLines", cex=0.7)
tn <- 2
kpDataBackground(kp, r0=tr.o*tn, r1=tr.o*tn+tr.i)
kpBars(kp, dd, y1=dd$y, r0=tr.o*tn, r1=tr.o*tn+tr.i, col="#AAFFAA", border="#66DD66")
kpRect(kp, chr="chr1", x0=5000000, x1=45000000, y0=0.2, y1=0.8, r0=tr.o*tn, r1=tr.o*tn+tr.i, col="#EEEEEE", border="#666666")
kpText(kp, chr="chr1", x=25000000, y=0.5, col="red", r0=tr.o*tn, r1=tr.o*tn+tr.i, labels="kpBars", cex=0.7)
tn <- 3
kpDataBackground(kp, r0=tr.o*tn, r1=tr.o*tn+tr.i)
kpRect(kp, dd, y0=dd$y-0.3, y1=dd$y, r0=tr.o*tn, r1=tr.o*tn+tr.i, col="#AAAAFF", border="#6666DD")
kpRect(kp, chr="chr1", x0=5000000, x1=45000000, y0=0.2, y1=0.8, r0=tr.o*tn, r1=tr.o*tn+tr.i, col="#EEEEEE", border="#666666")
kpText(kp, chr="chr1", x=25000000, y=0.5, col="red", r0=tr.o*tn, r1=tr.o*tn+tr.i, labels="kpRect", cex=0.7)
tn <- 4
kpDataBackground(kp, r0=tr.o*tn, r1=tr.o*tn+tr.i)
kpText(kp, dd, labels=as.character(1:50), cex=0.5, r0=tr.o*tn, r1=tr.o*tn+tr.i, col="#DDAADD")
kpRect(kp, chr="chr1", x0=5000000, x1=45000000, y0=0.2, y1=0.8, r0=tr.o*tn, r1=tr.o*tn+tr.i, col="#EEEEEE", border="#666666")
kpText(kp, chr="chr1", x=25000000, y=0.5, col="red", r0=tr.o*tn, r1=tr.o*tn+tr.i, labels="kpText", cex=0.7)
tn <- 5
kpDataBackground(kp, r0=tr.o*tn, r1=tr.o*tn+tr.i)
kpSegments(kp, dd, y0=dd$y-0.3, y1=dd$y, r0=tr.o*tn, r1=tr.o*tn+tr.i)
kpRect(kp, chr="chr1", x0=5000000, x1=45000000, y0=0.2, y1=0.8, r0=tr.o*tn, r1=tr.o*tn+tr.i, col="#EEEEEE", border="#666666")
kpText(kp, chr="chr1", x=25000000, y=0.5, col="red", r0=tr.o*tn, r1=tr.o*tn+tr.i, labels="kpSegments", cex=0.7)
tn <- 6
kpDataBackground(kp, r0=tr.o*tn, r1=tr.o*tn+tr.i)
kpArrows(kp, dd, y0=dd$y-0.3, y1=dd$y, r0=tr.o*tn, r1=tr.o*tn+tr.i, length=0.04)
kpRect(kp, chr="chr1", x0=5000000, x1=45000000, y0=0.2, y1=0.8, r0=tr.o*tn, r1=tr.o*tn+tr.i, col="#EEEEEE", border="#666666")
kpText(kp, chr="chr1", x=25000000, y=0.5, col="red", r0=tr.o*tn, r1=tr.o*tn+tr.i, labels="kpArrows", cex=0.7)
tn <- 7
kpDataBackground(kp, r0=tr.o*tn, r1=tr.o*tn+tr.i)
kpHeatmap(kp, dd, r0=tr.o*tn+tr.i/4, r1=tr.o*tn+tr.i-tr.i/4, colors = c("green", "black", "red"))
kpRect(kp, chr="chr1", x0=5000000, x1=45000000, y0=0.2, y1=0.8, r0=tr.o*tn, r1=tr.o*tn+tr.i, col="#EEEEEE", border="#666666")
kpText(kp, chr="chr1", x=25000000, y=0.5, col="red", r0=tr.o*tn, r1=tr.o*tn+tr.i, labels="kpHeatmap", cex=0.7)
tn <- 8
kpDataBackground(kp, r0=tr.o*tn, r1=tr.o*tn+tr.i)
kpPolygon(kp, dd, r0=tr.o*tn, r1=tr.o*tn+tr.i)
kpRect(kp, chr="chr1", x0=5000000, x1=45000000, y0=0.2, y1=0.8, r0=tr.o*tn, r1=tr.o*tn+tr.i, col="#EEEEEE", border="#666666")
kpText(kp, chr="chr1", x=25000000, y=0.5, col="red", r0=tr.o*tn, r1=tr.o*tn+tr.i, labels="kpPolygon", cex=0.7)
tn <- 9
kpDataBackground(kp, r0=tr.o*tn, r1=tr.o*tn+tr.i)
kpAbline(kp, h=c(0.25, 0.5, 0.75), v=start(dd), r0=tr.o*tn, r1=tr.o*tn+tr.i)
kpRect(kp, chr="chr1", x0=5000000, x1=45000000, y0=0.2, y1=0.8, r0=tr.o*tn, r1=tr.o*tn+tr.i, col="#EEEEEE", border="#666666")
kpText(kp, chr="chr1", x=25000000, y=0.5, col="red", r0=tr.o*tn, r1=tr.o*tn+tr.i, labels="kpAbline", cex=0.7)
高级功能
https://bernatgel.github.io/karyoploter_tutorial/ Plot Markers 将文本标签(text labels)在基因组(基因或任何其他特征)上进行定位,避免标签重叠 Plot Regions Plot Density 执行计算基因组上特征的密度计算并画图 Plot Coverage Plot BAM Coverage 绘制 BAM 文件中的 reads 密度 Plot Horizon plots Manhattan plots
Plot Density
library(BiocManager)
#install("BSgenome.Hsapiens.UCSC.hg19")
library(BSgenome.Hsapiens.UCSC.hg19)
regions <- createRandomRegions(nregions=10000, length.mean = 1e6, mask=NA, non.overlapping = FALSE)
# 如果使自己的矩阵,也是从data frame开始,转成toGRanges
a <- data.frame(regions)
head(a)
gains <- toGRanges(a)
kp <- plotKaryotype()
kpPlotDensity(kp, data=gains)
例子
资料上作者结合11个例子绘图
Gene expression results from DESeq2
数据准备
load(file = 'input_expressresult')
# 如果使自己的矩阵,也是从data frame开始,转成toGRanges
a <- data.frame(dm.genes)
head(a)
dm.genes <- toGRanges(a)
kpPlotMarkers标志
# na.last:If `TRUE', missing values in the data are put last
ordered <- dm.genes[order(dm.genes$padj, na.last = TRUE),]
kp <- plotKaryotype(genome="dm6")
# 取p值最小前10基因
kp <- kpPlotMarkers(kp, ordered[1:10], labels = names(ordered[1:10]), text.orientation = "horizontal")
kpPoints添加散点
filtered.dm.genes <- dm.genes[!is.na(dm.genes$padj)]
log.pval <- -log10(filtered.dm.genes$padj)
# 添加一列log.pval
mcols(filtered.dm.genes)$log.pval <- log.pval
filtered.dm.genes
sign.genes <- filtered.dm.genes[filtered.dm.genes$padj < 0.05,]
range(sign.genes$log.pval)
kp <- plotKaryotype(genome="dm6")
kpPoints(kp, data=sign.genes, y=sign.genes$log.pval,ymax=max(sign.genes$log.pval))
kpAxis添加y轴
fc.ymax <- ceiling(max(abs(range(sign.genes$log2FoldChange))))
fc.ymin <- -fc.ymax
cex.val <- sqrt(sign.genes$log.pval)/3
kp <- plotKaryotype(genome="dm6")
# cex点的大小,根据p值;y轴位置根据log2FoldChange
kpPoints(kp, data=sign.genes, y=sign.genes$log2FoldChange, cex=cex.val, ymax=fc.ymax, ymin=fc.ymin)
# y轴最大最小值
kpAxis(kp, ymax=fc.ymax, ymin=fc.ymin)
kpAddLabels(kp, labels = "log2 Fold Change", srt=90, pos=1, label.margin = 0.02, ymax=fc.ymax, ymin=fc.ymin)
col.over <- "#FFBD07AA"
col.under <- "#00A6EDAA"
points.top <- 0.8
top.genes <- ordered[1:20]
# 分配颜色
sign.col <- rep(col.over, length(sign.genes))
sign.col[sign.genes$log2FoldChange<0] <- col.under
kp <- plotKaryotype(genome="dm6", plot.type=2)
# Data panel 1
# 上面散点图占比80%
kpPoints(kp, data=sign.genes, y=sign.genes$log2FoldChange, cex=cex.val, ymax=fc.ymax, ymin=fc.ymin, r1=points.top, col=sign.col)
kpAxis(kp, ymax=fc.ymax, ymin=fc.ymin, r1=points.top)
kpAddLabels(kp, labels = "log2 FC", srt=90, pos=1, label.margin = 0.04, ymax=fc.ymax, ymin=fc.ymin, r1=points.top)
gene.mean <- start(top.genes) + (end(top.genes) - start(top.genes))/2
# kpSegments标志线倾斜度 x0, y0, x1 and y1 it will plot segments going from (x0, y0) to (x1, y1).
kpSegments(kp, chr=as.character(seqnames(top.genes)), x0=gene.mean, x1=gene.mean, y0=top.genes$log2FoldChange, y1=fc.ymax, ymax=fc.ymax, ymin=fc.ymin, r1=points.top)
kpPlotMarkers(kp, top.genes, labels = names(top.genes), text.orientation = "horizontal", r0=points.top)
#Data panel 2
kp <- kpPlotDensity(kp, data=dm.genes, window.size = 10e4, data.panel = 2)
基本参数
参考:
https://bioconductor.org/packages/devel/bioc/vignettes/karyoploteR/inst/doc/karyoploteR.html#introduction https://bernatgel.github.io/karyoploter_tutorial//Tutorial/PlotDensity/PlotDensity.html https://bernatgel.github.io/karyoploter_tutorial/
文末友情推荐
与十万人一起学生信,你值得拥有下面的学习班: