karyoploteR-染色体全局可视化

六月份的一个学徒自己有甲基化项目数据要处理,所以开始自学一些R包了

下面是2021六月份学徒的投稿

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/

文末友情推荐

与十万人一起学生信,你值得拥有下面的学习班:

(0)

相关推荐

  • ​​​​R语言学习笔记(五)——曼哈顿图

    iJournal 学术期刊信息查询 386篇原创内容 Official Account  ↑ ↑  关注iJournal,选刊快人一步  ↑ ↑  iJournal后台回复"2021学科&q ...

  • 使用biopython可视化染色体和基因元件

    基因组结构元件的可视化有多种方式,比如IGV等基因组浏览器中以track为单位的展示形式,亦或以circos为代表的圈图形式,比如在细胞器基因组组装中,基因元件常用圈图形式展示,示例如下 在biopy ...

  • 用人事物三个视角来串通知识、技术、技能!对学习构建全局认知!

     认知框架 2020-05-24优质教育领域创作者 所属专栏:建构符合认知机制的一套学习方法论,用它重构大脑里的认知框架! 用人事物来搞定知识.技术.技能.知识是物的层面,技术是事的层面,技能是人的 ...

  • 这是我看过最走心的数据可视化文章了,看完血赚!

    注:本文转载自 https://antv.alipay.com/zh-cn/vis/blog/vis-introduce.html 一.什么是数据可视化 科学可视化(Scientific Visual ...

  • 麻省理工科学家提出新冠基因可能整合入宿主染色体

    通知 备用:如果以后读者无法再阅读到本公号的更新,可以去网易搜索同名的网易号.海外的读者也可以下载Telegram , 然后再手机浏览器中打开链接"https://t.me/joinchat ...

  • 人类源流——人类Y染色体1

    分子人类学产生于二十世纪60年代,它是分子生物学与人类学交叉产生的边缘学科.从60年代开始,一些分子生物学家逐步将分子生物学技术引入人类学研究领域,试图通过研究人类DNA中所蕴藏的遣传信息来揭示整个人 ...

  • 人类源流——人类Y染色体2

    5.现代智人与古人类混血 约20万年前,智人的'亚当'和'夏娃'产生.粒线体DNA与化石证明现代人类大约于20万年前起源于东非.与其他动物相比,人具有高度发达的大脑,具有抽象思维.语言.自我意识以及解 ...

  • 人类源流——人类Y染色体3

    6.Y-C单倍群 Y染色体C-M130单倍群是Y-CF两个分支中的一支,C的地位与F相当.Y-C人群发现于除非洲以外的各个大陆古代人群中,是中亚.西伯利亚.北美和大洋洲一些土著部落的主流单倍群.在早期 ...

  • 人类源流——人类Y染色体5

    13.Y-Q单倍群 Y-Q系是远古时代唯一一个在全球范围内进行过扩张迁徙的单倍群,几乎走到过地球上的每个角落.除美洲的印第安人在500年前欧洲人进入之前因缺少其他人种竞争,并随着欧洲殖民者到来,同时带 ...

  • 人类源流——人类Y染色体6

    14.Y- R单倍群 印欧人种在分子人类学上主要属于Y-R系.分子人类学越来越精细的研究,正在越来越清晰地勾画出人类种族迁徙和文明发展史.虽然Y-R系诞生很早,大约27000年前就已经出现于中亚P集团 ...