使用Bioconductor画染色体示意图ideograms

之前介绍过如何使用circos来画染色体示意图ideograms。那么,在bioconductor中是否可以实现呢?肯定的说,是可以的,但是其生成的图像却有些差强人意。不过,总好过没有,而且,它还有很大的成长空间。

在Bioconductor中,用于画染色体示意图的有三种工具,它们分别是:quantsmooth,GenomeGraphs,以及ggbio。

首先介绍quantsmooth。我们先了解一下它是如何画单个的ideogram的。

> library(quantsmooth)
> plot(c(0,4),c(0,3),type="n",xaxt="n",yaxt="n",xlab="",ylab="")
> drawSimpleChrom(2,3,fill=c("p","q3"),col=c("red","blue"),orientation="v")

drawSimpleChrom

如果我们想简单地直接画出人类1号染色体ideogram应该怎么办呢?

> grid.newpage()
> grid.chromosome(1,units="bases",height=0.15)

grid.chromosome
有了grid的帮助,我们就可以使用viewport()来准确地控制ideogram所放的位置,然后生成比较复杂的示意图了。

那么如果我们想一次性把所有的染色体全画出来应该怎么办呢?

> # construct genomic positions
> CHR<-sample(22,40,replace=TRUE)  # Chromosomes
> MapInfo<-lengthChromosome(CHR,"bases")*runif(length(CHR)) # position on chromosome
> chrompos<-prepareGenomePlot(data.frame(CHR,MapInfo),paintCytobands = TRUE, organism="hsa")

prepareGenomePlot

那么GenomeGraphs包又有什么特点呢?

> library(GenomeGraphs)
> minbase <- 180292097
> maxbase <- 180492096
> ideog <- makeIdeogram(chromosome = 1)
> ideog2 <- makeIdeogram(chromosome = 2)
> ideog3 <- makeIdeogram(chromosome = 3)
> ideog4 <- makeIdeogram(chromosome = 4)
> gdPlot(list("1"= ideog, "2" = ideog2, "3" =ideog3, "4"=ideog4 ), 
+  minBase = minbase, maxBase = maxbase)

gdPlot
很明显,其绘制的示意图无论是形状,还是方式都与前者有很大的不同。它的绘制方式,似乎也限制了它的使用,不那么随心所欲。但是它却很容易实现对齐。

最后介绍ggbio。ggbio是一个非常强大的工具,从名字就可以看出它与ggplot2有着很大的关联。它甚至可以部分实现circos的功能(但是效果比较粗陋)。
我们还是先画一段chr1。

> library(ggbio)
> library(biovizBase)
> data(hg19IdeogramCyto)
> plotIdeogram(hg19IdeogramCyto, "chr1", aspect.ratio = 1/20, zoom.region = c(1e+07,
+ 5e+07))

plotIdeogram

我们可以通过biovizBase来调取更多物种的ideograms。

> library(biovizBase)
> obj <- getIdeogram()
Please specify genome 
 
  1: hg19       2: hg18       3: hg17       4: hg16       5: vicPac1    6: dasNov3    7: otoGar3    8: papHam1    9: papAnu2   10: felCat5
 11: felCat4   12: felCat3   13: panTro4   14: panTro3   15: panTro2   16: panTro1   17: bosTau7   18: bosTau6   19: bosTau4   20: bosTau3
 21: bosTau2   22: canFam3   23: canFam2   24: canFam1   25: turTru2   26: loxAfr3   27: nomLeu2   28: nomLeu1   29: gorGor3   30: cavPor3
 31: eriEur1   32: equCab2   33: equCab1   34: dipOrd1   35: triMan1   36: calJac3   37: calJac1   38: pteVam1   39: myoLuc2   40: mm10   
 41: mm9       42: mm8       43: mm7       44: micMur1   45: hetGla2   46: hetGla1   47: monDom5   48: monDom4   49: monDom1   50: ponAbe2
 51: ailMel1   52: susScr3   53: susScr2   54: ochPri2   55: ornAna1   56: oryCun2   57: rn5       58: rn4       59: rn3       60: rheMac3
 61: rheMac2   62: proCap1   63: oviAri1   64: sorAra1   65: choHof1   66: speTri2   67: saiBol1   68: sorAra1   69: sarHar1   70: echTel1
 71: tupBel1   72: macEug2   73: cerSim1   74: gadMor1   75: melUnd1   76: galGal4   77: galGal3   78: galGal2   79: latCha1   80: fr3    
 81: fr2       82: fr1       83: petMar2   84: petMar1   85: anoCar2   86: anoCar1   87: oryLat2   88: geoFor1   89: oreNil2   90: chrPic1
 91: gasAcu1   92: tetNig2   93: tetNig1   94: melGal1   95: xenTro3   96: xenTro2   97: xenTro1   98: taeGut1   99: danRer7  100: danRer6
101: danRer5  102: danRer4  103: danRer3  104: ci2      105: ci1      106: braFlo1  107: strPur2  108: strPur1  109: apiMel2  110: apiMel1
111: anoGam1  112: droAna2  113: droAna1  114: droEre1  115: droGri1  116: dm3      117: dm2      118: dm1      119: droMoj2  120: droMoj1
121: droPer1  122: dp3      123: dp2      124: droSec1  125: droSim1  126: droVir2  127: droVir1  128: droYak2  129: droYak1  130: caePb2 
131: caePb1   132: cb3      133: cb1      134: ce10     135: ce6      136: ce4      137: ce2      138: caeJap1  139: caeRem3  140: caeRem2
141: priPac1  142: aplCal1  143: sacCer3  144: sacCer2  145: sacCer1  
 
Selection: 40
Loading...
Done
> head(obj)
GRanges with 6 ranges and 2 metadata columns:
      seqnames               ranges strand |     name gieStain
         <Rle>            <IRanges>  <Rle> | <factor> <factor>
  [1]     chr1 [       0,  8840440]      * |      qA1  gpos100
  [2]     chr1 [ 8840440, 12278390]      * |      qA2     gneg
  [3]     chr1 [12278390, 20136559]      * |      qA3   gpos33
  [4]     chr1 [20136559, 22101102]      * |      qA4     gneg
  [5]     chr1 [22101102, 30941543]      * |      qA5  gpos100
  [6]     chr1 [30941543, 43219933]      * |       qB     gneg
  ---
  seqlengths:
    chr1 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19  chr2  chr3  chr4  chr5  chr6  chr7  chr8  chr9  chrX  chrY
      NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
> p <- lapply(1:19, function(.ele) plotIdeogram(obj, paste("chr",.ele,sep="")))
> alignPlots(p)

alignPlots

还可以绘制类似circos的图像。
Screen shot 2013-03-12 at 11.32.29 PM

在这里,我只是简单介绍了可以实现ideograms的包,并没有对每个函数进行深入的讲解。但我相信,有一这些,加上对原代码的仔细研读,一定会绘制出你想要的ideograms的。

13 thoughts on “使用Bioconductor画染色体示意图ideograms

  1. 博主,我对于quantsmooth包的prepareGenomePlot函数有一些问题,我想画小鼠的染色体图,前面数据已改为小鼠的chrompos,但是当我改organism=”mmu”后,就无法画出paintCytobands = TRUE的效果了,求解答。
    我的目的是把一套差异基因MAP到全基因组位置上,用染色体图呈现,不想要circos,圆圈不直观,博主有不有什么更好的建议?

  2. 博主,我还有个问题,quantsmooth包按上面的代码进行后,我发现通过prepareGenomePlot获得的chrompos[1:10,]与data.frame(CHR,MapInfo)[1:10,]的CHR发生了改变,但是位置信息没有变化,也就是说他完全把位置信息弄错了,这算BUG吗?

  3. 您好,請問ggbio及quantsmooth是否可以畫whole genome並加上cytoband呢?
    若可以的話不知道有沒有sample呢?
    謝謝您~

  4. 您好,请问我现在有组转录组分析结果,是同一材料对不同胁迫相应。 我想把这两组画在20条染色体上,基因组已经测序,基因的物理位置也知道。我想上调画在染色体上面,下调基因在下面,不同情况用不同的颜色,请问以你的经验,用什么包或者如何才能实现?多谢了。

  5. 博主,这个画出来的染色体的宽度都是一样的,可以调节么?我用的Gviz也是,在同一张大小的画布上,每一条染色体的宽度都是一样的~求助,谢谢!

发表评论

电子邮件地址不会被公开。 必填项已用*标注