| 查看: 1023 | 回复: 0 | |||
[交流]
线性基因组可视化神来之笔---karyoploteR
|
|
线性基因组可视化神来之笔---karyoploteR --光俊 随着人类基因组逐渐被破译，一张生命之图将被绘制，人们的生活也将发生巨大变化。今天我们给大家介绍一个包，对于性基因组可视化来说，可谓是神来之笔。 记得多年前，刚接触高通量测序，为了画一张覆盖度图，花了一周时间写代码。那时用的办法非常笨，不过总算能够完成操作。我们可以描述一下大概的过程：我们是从bam文件出发，先把bam文件转换为bed文件，有了bed文件以后，在对染色体划分区段，比如10kb为一段，然后计算每个区段的覆盖深度。然后在应用复杂的R代码，先把23挑染色体画出来，然后再把每个区段的深度plot上去。现在想想，听起来就觉得复杂。现在好了，这些复杂的工作，一个命令就搞定，这就是我们今天要介绍的karyoploteR包。 我们只要准备简单的输入文件，就能完成漂亮的基因图的绘制。比如我们做GEO和TCGA数据挖掘的时候，得到了差异表达的表格，我们只需要利用差异表达里面基因名字和logFC，就能完成一张漂亮的基因组图。 再比如，我们手头有CNV或者基因融合的结果，我们就可以利用最原始的数据，绘制基因组图。 更神奇的是，我们可以直接输入bam文件，就能得到漂亮的圈图。是不是很神奇，以前公司做项目，需要两三千的图形，现在只要几个命令就能搞定，想学的赶紧加入我们生信自学网。 下面我们看一段脚本，非常简单将bam直接可视化。 setwd("C:\\Users\\lexb4\\Desktop\\karyoploteR\\12.kpPlotBAMDensity"  library(karyoploteR) bamFile <- "input.bam" window.size <- 1e4 #compute the density with 10kb windows pdf(file="kpPlotBAMDensity.pdf",width=10,height=7) kp <- plotKaryotype(genome="hg38" kp <- kpPlotBAMDensity(kp, data = bamFile, window.size = window.size, r0=0, r1=1, ymax=600, col="darkorchid" kpAddBaseNumbers(kp, tick.dist=10000000, minor.tick.dist=1000000) dev.off() 接下来我们看看karyoploteR的效果图，如下： |
回复此楼» 猜你喜欢
垃圾破二本职称评审标准
已经有18人回复
-
职称评审没过，求安慰
已经有53人回复
-
毕业后当辅导员了，天天各种学生超烦
已经有5人回复
-
26申博自荐
已经有3人回复
-
A期刊撤稿
已经有4人回复
-
投稿Elsevier的Neoplasia杂志，到最后选publishing options时页面空白，不能完成投稿
已经有22人回复
-
EST投稿状态问题
已经有7人回复