?想到一個(gè)單詞,karyotype: 核型; 染色體組型.
?基因組瀏覽器(ucsc genome browser, igv這類(lèi))可以交互式地展示基因組數(shù)據(jù)的可視化結(jié)果,染色體的區(qū)段信息等等丽蝎。
?有時(shí)我們也需要根據(jù)自己的數(shù)據(jù)畫(huà)一個(gè)這類(lèi)的圖,該怎么辦呢?這兩天就遇到這個(gè)問(wèn)題了。組里面的一位前輩讓我根據(jù)一些CNV的位置信息,畫(huà)出它們?cè)谌旧w上面的示意圖。如下:
?起初想了下circos圖似乎也能表示,比如某個(gè)區(qū)域一個(gè)樣本缺失記為1,兩個(gè)樣本缺失記為2......最后畫(huà)柱狀圖的圈。但想想又不太合適,因?yàn)闆](méi)人會(huì)在circos圖里面只畫(huà)一兩個(gè)圈,太少了节腐,而且一個(gè)圈畫(huà)24條染色體太緊湊,CNV的位置看不清楚掷空。所以最后還是決定染色體全都鋪開(kāi)展示酿傍。找了很久才找到這個(gè)R包——karyoploteR亏较,還比較新莺褒,是17年發(fā)的文章,相關(guān)的教程也很少(只找到一兩篇)雪情。今天只用到了這個(gè)包里面的幾個(gè)小功能遵岩,下面開(kāi)始介紹。
1. 安裝
這是我第一次用Bioconductor上的R包巡通,所以先簡(jiǎn)單看了一下安裝Bioconductor上R包的方法尘执。我知道的是這么幾種:install.packages()、source()+biocLite()宴凉、BiocManager::install()誊锭。
R語(yǔ)言 Bioconductor包的安裝指南
不用biocLite安裝Bioconductor包——Y叔叔 biobabble
再結(jié)合官網(wǎng)說(shuō)明,安裝問(wèn)題就解決了: http://bioconductor.org/packages/release/bioc/html/karyoploteR.html
(library(karyoploteR)的時(shí)候弥锄,可能會(huì)顯示沒(méi)有某些包,補(bǔ)上就好了)
2. 使用
快速入門(mén): http://bioconductor.org/packages/release/bioc/vignettes/karyoploteR/inst/doc/karyoploteR.html
詳細(xì)手冊(cè): http://bioconductor.org/packages/release/bioc/manuals/karyoploteR/man/karyoploteR.pdf
thus can be used to plot genomic data coming from anywhere. The only exception to this are the ideograms cytobands, that by default are plotted using predownloaded data from UCSC.
可用于繪制來(lái)自任何地方的基因組數(shù)據(jù)温治。唯一的例外是染色體G顯帶信息图仓,默認(rèn)情況下使用來(lái)自UCSC預(yù)先下載的數(shù)據(jù)繪制救崔。
2.1 kpPlotRegions
這個(gè)繪圖函數(shù)剛好可以用來(lái)畫(huà)我需要的圖六孵,即CNV的位置信息。
舉個(gè)例子:
gains <- toGRanges(data.frame(chr=c("chr17"), start=c(4000000),end=c(80000000)))
losses <- toGRanges(data.frame(chr=c("chr3"), start=c(80000000),end=c(170000000)))
kp <- plotKaryotype(genome="hg19")
kpPlotRegions(kp, gains, col="#FFAACC")
kpPlotRegions(kp, losses, col="#CCFFAA")
data.frame(): 創(chuàng)建一個(gè)數(shù)據(jù)框;
toGRanges(): 將接收到的表示染色體位置的“表格”轉(zhuǎn)換為R能識(shí)別的格式——GRanges
;
用法
: toGRanges(A, ..., genome=NULL), A可以是(一個(gè)或幾個(gè))數(shù)據(jù)框, BED格式的文件, GRanges對(duì)象;
GRanges
: 存放基因組位置和相關(guān)注釋的向量。 向量中的每個(gè)元素由序列名稱主巍,區(qū)間孕索,鏈和可選數(shù)據(jù)列(例如score,GC含量等)組成散怖。
plotKaryotype(): 創(chuàng)建一個(gè)只含有染色體核型和染色體名稱的圖(最原始的圖);
用法
:
- plotKaryotype(genome="hg19", plot.type=1, chromosomes=c("chr1","chr2"));
- genome可以是UCSC中認(rèn)可的基因組名稱镇眷,比如hg19, mm10等等欠动,也可以是GRanges對(duì)象;
- plot.type有1...7可以選钝的,默認(rèn)是1硝桩,具體含義如下:
- plot.type=1 所有染色體依次水平鋪開(kāi),每個(gè)染色體上面有一個(gè)畫(huà)圖面板啼肩;
- plot.type=2 所有染色體依次水平鋪開(kāi)祈坠,每個(gè)染色體上下各有一個(gè)畫(huà)圖面板赦拘;
- plot.type=3 所有染色體在一條水平線上鋪開(kāi)芬沉,每個(gè)染色體上下各有一個(gè)畫(huà)圖面板;
......
其它的幾個(gè)也好理解蹋艺,但只能水平畫(huà)染色體算不算一個(gè)小缺點(diǎn)呢?
- chromosomes=c()選擇需要展示的染色體民效,默認(rèn)情況下全都展示;
- zoom = A可以用來(lái)定義顯示區(qū)間畏邢,A可以是數(shù)據(jù)框, BED格式的文件, GRanges對(duì)象州叠,但只有第一個(gè)區(qū)間會(huì)被使用;
- cytobands=A可以用來(lái)畫(huà)自定義的染色體條帶(涉及到細(xì)胞遺傳學(xué)的概念),A可以是GRanges對(duì)象,數(shù)據(jù)框, BED交煞。不過(guò)表示條帶的這一列必須以 "gieStain"命名素征,值為gneg 萝挤,gpos25 ,gpos50怜珍,gpos75今豆,gpos100呆躲,acen(#red, centromere)插掂,stalk(#indented region, region with repeats)箩祥,gvar等。在查找這一部分內(nèi)容時(shí)谢揪,意外發(fā)現(xiàn)MATLAB可以畫(huà)染色體色帶圖,學(xué)過(guò)數(shù)學(xué)建模的小伙伴可以了解一下:https://www.mathworks.com/help/bioinfo/ref/chromosomeplot.html
- main=""可以定義一個(gè)圖片名稱;
以上就是plotKaryotype()的用法患民。再來(lái)看看kpPlotRegions()匹颤。
kpPlotRegions(): 根據(jù)GRanges對(duì)象提供的位置信息,在染色體上畫(huà)矩形赦肃。
用法
:
- kpPlotRegions(karyoplot, data, data.panel=1, r0=NULL, r1=NULL, col="black", border=NULL, avoid.overlapping=TRUE, num.layers=NULL, layer.margin=0.05);
- karyoplot: plotKaryotype()的返回值,也就是只含有染色體核型和染色體名稱的圖厅各;
- data=A表示區(qū)間讯检,GRanges對(duì)象,數(shù)據(jù)框, BED都行投放;
- data.panel=1表示在染色體的上面板畫(huà)涝桅,2表示在染色體的下面板畫(huà)冯遂,具體取值要結(jié)合plotKaryotype()函數(shù)中的plot.type參數(shù)來(lái)選擇蛤肌;
- r0, r1表示針對(duì)這一行畫(huà)圖命令,利用面板(范圍是0-1)的哪一部分炒俱,比如r0=0, r1=0.5就是只利用一半面板來(lái)畫(huà)圖爪膊;
- col, border表示填充顏色和邊界顏色;
- avoid.overlapping為真僵芹,若區(qū)間有重疊則堆疊起來(lái)小槐,為假若區(qū)間有重疊則取并集只顯示一個(gè)矩形,默認(rèn)為真件豌;
- num.layers表示將一個(gè)大面板(panel)劃分為幾個(gè)小面板(layers)來(lái)畫(huà)矩形茧彤,默認(rèn)為NULL曾掂,會(huì)根據(jù)最大重疊次數(shù)來(lái)自動(dòng)劃分panel珠洗;
- layer.margin小面板之間的間隔若专;
2.2 實(shí)戰(zhàn)
這一個(gè)繪圖函數(shù)應(yīng)該講得比較清楚了许蓖,下面來(lái)應(yīng)用一下自阱。
dels <- read.table("del.txt",header = T)
dups <- read.table("dup.txt",header = T)
gains <- toGRanges(dups)
losses <- toGRanges(dels)
my_kp <- plotKaryotype(genome="hg19",plot.type = 2, chromosomes=c("chr1","chr2"),main="my test")
kpPlotRegions(my_kp, gains, col="blue",border = "white",r0 = 0, r1 = 1,num.layers = 5)
kpPlotRegions(my_kp, losses, col="red",border = "white",data.panel = 2)
到這兒這一部分就結(jié)束了,這個(gè)R包的其他繪圖函數(shù)下次再講赃额。
?另外爬早,看完這些可能會(huì)有一個(gè)疑問(wèn)筛严,其他的一些物種可以畫(huà)嗎桨啃?按照官方的說(shuō)明照瘾,在UCSC上能找到的應(yīng)該都能畫(huà)丧慈。不過(guò)逃默,一些植物物種應(yīng)該就只能根據(jù)自己的數(shù)據(jù)繪制核型圖了完域。