写在2018年的第一天

写在2018年的第一天

先说一件比较尴尬的事，昨天晚上肚子不舒服，于是蹲了好久，等出来才发现一年就这样过去了。时间节点，如情人节、除夕夜，自己的生日，对我而言并没有多大意义（女朋友的生日除外）。

我们并不需要专门找一天来下定决心重新开始，如果这件事情非常重要你早已经开始，这件事情无关紧要你就会一推再推。但是我们也的确需要这么一天来回顾过去，展望未来。大脑需要一种仪式感，一种暗示，或者说是一种欺骗吧。

我的2017年

2017年刚开始时，我还是研一新生，还在迷惘未来何去何从。当时的情况是一轮轮转的老板要离开原单位去其他学校，先是告诉我可以继续留下来，后来却否定了这一可能性。二轮的老板看我情况太复杂，也不知道该不该留我，让我和所里继续抗争着。后来我去做第三轮，跑到一个院士下面说搞生信，院士让我先留下学习这2周。实际情况是，我不太懂如何和别人打交道，一直觉得自己格格不入。他们的主力基本走完了，留下的工作人员大部分工作是在改图，原本负责带我的老师正在准备投递文章。我每天早点7点30出发，挤地铁过去，边学边翻译bioconductor的教程。由于那边只供应午饭，所以每天晚上7点左右就回来，半路点外卖，到寝室开吃，然后晚上和大学同学打打DOTA2。有一天，我听着《春风十里》，出了地铁口，看着夜空，吹过的峰让我又紧了紧外衣，不知何去何从。最后我没有在院士那里留下来，我知道问题在我。我得感谢现在的老板，给了我一个去处，又因为实验室还不够大，所以让我继续去自学生信。

为了搞定师姐的一个EMS诱变群体的基因定位项目，我接触了生信媛的幕后大大，也是一位自学生信的师姐，开始在生信媛发我的学习笔记，成为生信媛的一名小编。在2017年春季的某一天，生信媛的大部分小伙伴跑去幕后大大的家里蹭饭，通过那次聚会中我才认识了那些数字信息后的物理意义上的人。或许是因为我稳定的更新频率，或许是我把统计学的作业发在生信媛，生信媛的订阅迎来了第一个1000人。后来也没有刻意转发到各种生信相关群里，每天增加的人也特别的稳定，到10月份达到了5000人。

当然生信媛的微信推送文章一般都会同步到生信技能树论坛，也就顺理成章加入生信技能树论坛的作者群。某一天在承包下生信工作资讯速递这个栏目后，我接触了生信技能树的Jimmy。又由于自己写的比较勤快，让Jimmy看到自己年轻时候的样子，于是把我拉入了生信技能树公众号的小编群，认识了更多的人，如果子、思考问题的熊、Panda姐等。Jimmy有很多想法，比如说写书，做视频课，搞论坛，结果还是先开展了武汉的线下培训。由于线下培训，我又认识了一批人，比如说武汉的线下团队。并且还让Jimmy从澳门给我带来了一台macbook pro，圆了我多年的梦想，也感谢我那个土豪师兄借钱。

因此，在物理世界中，我不懂如何和陌生人打交道。在网络世界中，其实也没有好到哪里去。我唯一交流的方式，就是把自己的学习过程记录下来，写成笔记发到互联网上，当别人看到我的文字时，就知道我差不多是什么样的人。我和陌生人的寒暄方式就是提问和解答，有时候我问别人，有时候别人问我。如果有幸能够持续交流下去，那么在我的世界里就多了一个认识的人。如果这些人还刚好能在现实世界里碰面，那么我或许就会多几个朋友。

我的2018年目标

在2017年的年初时候，我曾经许下很多目标

  • 得到的通往财富自由之路每周总结
  • 天天用英语每日笔记
  • 读书做笔记
  • 做完200道生信编程题
  • 跟随jimmy的直播我的基因组

实际情况是这些目标一个我都没有完成。当然这可以从《意志力》这本书找到答案，设下的目标越多，意志力在决定目标的时候的消耗就越多，于是目标就很不可能完成。于是2018年我就放弃制定那么多目标了，我仅仅只有一个目标，就是在2018年争取发一篇文章，发在bioinformatics上。

其他

从typora转移到VNOTE

因为之前一直用的是Typora进行写作，笔记管理的方式就是新建文件夹，也是下图这种情况。

笔记本目录

从轻度写作而言，Typora基本上够用了，但是它并不是一个笔记管理工具。熊的原话是，就本地文件的管理而言，Typora和VNote之间相差一到两个印象笔记。并且从LTF(List, Tag, Filter)角度来看，我只有list，而没有tag，更不要说filter了，因此查找笔记全靠记忆或者靠搜索引擎（我的学习笔记基本上都发布在互联网了）。

为了减少我转移笔记的手工操作，我分析了下VNOTE的笔记管理模式，发现核心文件就两个

  • _v_imaes: 用于存放本地图片
  • _vnote.json: 用于记录文件信息

由于并且我的图片存储习惯和VNOTE一致，也就是为不同笔记本建立单独的图片目录，我的是assets，而VNote是_v_images，因此我在添加笔记本的时候，只做了如下简单的更改。

修改配置

最终结果如下

最终效果

因此，对于使用相对路径的我而言，我并不是从typora转移到VNote，而只是让我的之前的笔记本能够被VNote管理而已。

生信软件工具箱

「BioNano系列」如何进行cmap之间的比对

BioNano以cmap格式存放光学图谱，为了评估基因组的组装质量或者了解光学图谱中冗余情况(高杂合基因组组装结果偏大)，我们就需要进行cmap之间的比较。

CMAP间比对

Solve套件提供了runCharacterize.py脚本封装了RefAligner，用于进行CMAP之间的比对。

1
2
3
4
5
6
7
python2.7 runCharacterize.py \
    -t RefAligner的二进制文件路径 \
    -q 用于比对的CMAP \
    -r 参考CMAP \
    -p Pipeline文件路径\
    -a 参数配置文件.xml \
    -n 线程数,默认4

需要注意的是-p和-a参数的设置。-p是Pipeline的文件位置，比如说我的Solve安装在/opt/biosoft/Solve3.4_06042019a，那么参数设置为 -p /opt/biosoft/Solve3.4_06042019a/Pipeline/06042019。 而-a则是要在/opt/biosoft/Solve3.4_06042019a/RefAligner/8949.9232rel/目录下选择合适的xml文件。比如你的CMAP是Irys平台，那么你可以考虑用optArguments_nonhaplotype_irys.xml.

以最新发表的辣椒的光学图谱为例，该物种有比较高的杂合度，组装结果偏大，我们可以通过自比对来寻找冗余区域, 

1
2
3
4
5
6
7
8
9
# 下载CMAP
wget https://submit.ncbi.nlm.nih.gov/ft/byid/o62junnn/piper_nigrum_no_rcmap_refinefinal1.cmap
# 自比对
python /opt/biosoft/Solve3.4_06042019a/Pipeline/06042019/runCharacterize.py \
	-t /opt/biosoft/Solve3.4_06042019a/RefAligner/8949.9232rel/RefAligner \
	-q piper_nigrum_no_rcmap_refinefinal1.cmap \
	-r piper_nigrum_no_rcmap_refinefinal1.cmap \
	-p /opt/biosoft/Solve3.4_06042019a/Pipeline/06042019 \
	-a /opt/biosoft/Solve3.4_06042019a/RefAligner/8949.9232rel/optArguments_nonhaplotype_saphyr.xml -n 64

最终会在当前文件下生成一个alignRef文件夹，其中结果是q.cmap,r.cmap和xmap的文件可以用于上传到BioNano Access上进行展示。下图就是一个冗余实例，可以把图中较短的图谱删掉

冗余

基因组回帖

为了将基因组回帖到CMAP上，需要先将基因组的fasta格式转成CMAP格式，参数如下

1
perl fa2cmap_multi_color.pl -i 输入FASTA -e 酶1 通道1 [酶2 通道2]

其中一个最重要的参数就是酶切类型。例如我需要将序列回帖到用Nt.BspQI酶切组装的光学图谱上，因此运行参数如下

1
perl /opt/biosoft/Solve3.4_06042019a/HybridScaffold/06042019/scripts/fa2cmap_multi_color.pl -i athaliana.fa -e BspQI 1

最后的athaliana_BSPQI_0kb_0labels.cmap就是模拟酶切的CMAP序列。

之后将模拟酶切的结果回帖到实际的CMAP

1
2
3
4
5
6
7
python /opt/biosoft/Solve3.4_06042019a/Pipeline/06042019/runCharacterize.py \
    -t /opt/biosoft/Solve3.4_06042019a/RefAligner/8949.9232rel/RefAligner \
    -q athaliana_BSPQI_0kb_0labels.cmap \
    -r kbs-mac-74_bng_contigs2017.cmap \
    -p /opt/biosoft/Solve3.4_06042019a/Pipeline/06042019 \
    -a /opt/biosoft/Solve3.4_06042019a/RefAligner/8949.9232rel/optArguments_nonhaplotype_saphyr.xml \
    -n 64

最终会在当前文件下生成一个alignRef文件夹，其中结果是q.cmap,r.cmap和xmap的文件.

R

答扪心自问，meme几何？

在2018-04-06，Y叔推送了一篇文章扪心自问，meme几何？。从一个不到146行的meme.R出发，提了5个问题。让检查我们的ggplot2的理解程度。

从我第一次接触ggplot2开始，至今差不多过去了2年多时间。Y叔推荐的「ggplot2: 数据分析和图形语法」和「R绘图系统」也被我翻得书页泛黄，加上近日又在学习Hadley写的extending-ggplot2，近日终于有所悟，尝试解答Y叔提出的几个问题。

1. meme的输出不是图？为什么能画出图？

这一题，Y叔在R绘图系统的书评中给出了答案。

感谢译者送我的签名版，这是最全面介绍R绘图系统的书，没有之一。Base因为有大量的绘图函数还大量在使用，但做为新人学习，必须是grid系统！因为图是对象，可以操作，只有在需要渲染成静态图片的时候才产生图片。

更加详细一点，核心语句在于imageGrob <- rasterGrob(x)将读取的图像转成了grid的Grob对象，之后在此基础上构建了p，它继承了meme和recordeplot。这里的继承是关键词，也就是第二题的答案所在。

2. 为什么+可以改变图的内容和状态？

解答这一题需要两个关键知识

  1. R语言是一个面向对象编程的语言，里面有一类泛型函数，可以根据你的对象类型自动调用对应的函数。
  2. +是函数。

Y叔为meme对象专门定义了一个泛型函数+.meme, 因此在调用+的时候，也就是调用了+.meme函数。

3. 为什么ggsave能识别meme对象

这一题是讨论ggsave的本质，如果你直接在命令行里敲ggsave，他会输出ggsave的源代码，倒数第二句就是答案所在。

1
grid.draw(plot)

ggsave在绘图商调用的是grid.draw， 这是用来绘制一个grob对象。而无论是cowplot，还是meme, 它们都建立在grid系统下，也就能够用grid.draw画出来。

而如果你调用的是print(meme)，那么泛型函数会尝试调用print.grob

4. 为什么使用传统的出图方式来画meme，在循环中需要显示print(object)？而ggsave则不用？到底区别在那里？

这个问题稍微比较复杂， 我们需要先来实际的代码进行演示。

下面的for循环中，图形设备中不会出现图片，并且test.pdf打开的时候会显示图形损坏。

1
2
3
4
5
6
7
8
9
10
11
12
13
library(meme)
u <- system.file("angry8.jpg", package="meme")

for ( i in 1:10){
  meme(u, "code", "all the things!")
}

# 图片输出到新的pdf中
pdf("test.pdf")
for ( i in 1:10){
  meme(u, "code", "all the things!")
}
dev.off()

而下面的代码中，图形设备会打印图片，并且test2.pdf能出现图片。

1
2
3
4
5
6
7
8
9
10
11
for ( i in 1:10){
  p <- meme(u, "code", "all the things!")
  print(p)
}
# 图片输出到新的pdf中
pdf("test2.pdf")
for ( i in 1:10){
  p <- meme(u, "code", "all the things!")
  print(p)
}
dev.off()

为什么要在for循环里要用到print才行呢？

我们在R的控制台(console)运行meme时，实际上R会给你调用对应print函数答应。而在for循环中，它不会调用print。因此你必须要显示的调用print才行。

In a loop, automatic printing is turned off, as it is inside a function

参考: https://stackoverflow.com/questions/4716152/why-do-r-objects-not-print-in-a-function-or-a-for-loop

5. 为什么meme对象能够被ggimage和cowplot识别？

Y叔说的马甲其实就是指meme继承了recordedplot，不过现在版本的cowplot似乎搞不定

1
2
3
4
5
6
7
8
9
cowplot::plot_grid(p, p, ncol=1, labels = c("A", "B"))
# 报错如下
Error in value[[3L]](cond) : 
  invalid "recordedplot": Incompatible graphics state
In addition: Warning messages:
1: In restoreRecordedPlot(x, reloadPkgs) :
  snapshot recorded in different R version (pre 3.3.0)
2: In doTryCatch(return(expr), name, parentenv, handler) :
  snapshot recorded with different graphics engine version (pre 11 - this is version 12)

而为什么ggimage能够识别呢？ggimage是创建了一个专门的geom_image图层，为此Y叔利用ggproto，基于grid系统创造了一个GeomImage 类。这个图层就是用来绘制图片。

R

ggplot2高级:构建自己的图层

这部分内容是Extending ggplot2的学习笔记，大部分内容都是原文的简单翻译。

所有的ggplot2对象都建立自”ggproto”这套面向对象编程系统，因此想要创建出自己的一套图层，而不是简单的对已有图层进行累加，那么就需要学习”ggproto”。

创建新的stat

最简单的stat

我们会从一个最简单的stat开始: 根据已有的一组点，用一个凸壳(convex hull)包围他。

第一步，我们创建一个继承自Stat的”ggproto”对象

1
2
3
4
5
6
7
StatChull <- ggproto("StatChull", Stat,
                     compute_group = function(data, scales){
                       data[chull(data$x, data$y), , drop=FALSE]
                     },
                     required_aes = c("x","y")
)

在”ggproto”函数中，前两个是固定项，分别是类名和继承的”ggproto”类。而后续内容则是和你继承的类相关，例如compute_group()方法负责计算，required_aes则列出哪些图形属性(aesthetics)必须要存在，这两个都继承自Stat，可以用?Stat查看更多信息。。

第二步，我们开始写一个图层。由于历史设计原因，Hadley将其称作stat_()或geom_()。但实际上，Hadley认为layer_()可能更准确些，毕竟每一个图层都或多或少的有”stat”和”geom”。

所有的图层都遵循相同的格式，即你在function中声明默认参数，然后调用layer()函数，将...传递给params参数。在...的参数既可以是”geom”的参数(如果你要做一个stat封装)，或者是”stat”的参数(如果你要做geom的封装)，或者是将要设置的图形属性. layer()会小心的将不同的参数分开并确保它们存储在正确的位置:

1
2
3
4
5
6
7
8
9
10
stat_chull <- function(mapping = NULL, data = NULL, geom = "polygon",
                       position = "identity", na.rm = TRUE, show.legend = NA,
                       inherit.aes = TRUE, ...){
  layer(
    stat = StatChull, data = data, mapping = mapping, geom = geom,
    position = position, show.legend = show.legend, inherit.aes = inherit.aes,
    params = list(na.rm = na.rm, ...)
    
  )
}

(注, 在写R包的时候要注意用ggplot2::layer()或在命名空间中导入layer(), 否则会因找到函数而报错)

当我们写好了图层函数后，我们就可以尝试这个新的”stat”了

1
2
3
ggplot(mpg, aes(displ, hwy)) +
  geom_point() +
  stat_chull(fill = NA, colour= "black")

simple-stat

(后续我们会学习如何通过设置”geom”的默认值，来避免声明fill=NA)

一旦我们构建了这种基本的对象，ggplot2将会给我们带来极大的自由。举个例子，ggplot2自动保留每组中不变的图形属性，也就是说你可以分组绘制一个凸壳:

1
2
3
ggplot(mpg, aes(displ, hwy, colour = drv)) + 
  geom_point() + 
  stat_chull(fill = NA)

add group chull

我们还可以覆盖默认的图层，来以不同的形式展现凸壳:

1
2
3
ggplot(mpg, aes(displ, hwy)) +
  stat_chull(geom = "point", size = 4, colour = "red") +
  geom_point()

different chull

Stat参数

一个更加复杂的”stat”会做一些计算。我们可以通过实现一个简单版本的geom_smooth来了解。我们将会创建一个新的图层StatLm（继承自Stat）和一个的图层函数stat_lm():

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
# 基于ggproto创建StatLm
StatLm <- ggproto("StatLm", Stat,
                  required_aes = c("x", "y"),
                  compute_group = function(data, scales){
                    rng <- range(data$x, nr.rm = TRUE)
                    grid <- data.frame(x = rng)
                    mod <- lm(y ~ x, data = data)
                    grid$y <- predict(mod, newdata = grid)
                    grid
                  }
                  )
# 创建图层函数
stat_lm <- function(mapping = NULL, data = NULL, geom = "line",
                    position = "identity", na.rm = FALSE, show.legend = NA,
                    inherit.aes = TRUE, ...){
  layer(
    stat = StatLm, data = data, mapping = mapping, geom = geom,
    position = position, show.legend = show.legend, inherit.aes = inherit.aes,
    params = list(na.rm = na.rm, ...)
  )
}

调用我们写的stat_lm()图形，检查下效果

1
2
3
ggplot(mpg, aes(displ, hwy)) + 
  geom_point() + 
  stat_lm()

liner model

StatLm缺少参数不太灵活，只能做单一线性拟合。最好是允许用户能够自由修改模型公式和创建图层所需要的数据量。为了实现这一需求，我们在compute_group()增加了一些参数，代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
# 增加了参数n和formula
StatLm2 <- ggproto("StatLm2", Stat,
                  required_aes = c("x", "y"),
                  compute_group = function(data, scales, params, 
                                           n = 100, formula = y ~x){
                    
                    rng <- range(data$x, na.rm = TRUE)
                    grid <- data.frame(x = seq(rng[1], rng[2],length = n))
                    
                    mod <- lm(formula, data = data)
                    grid$y <- predict(mod, newdata = grid)
                    grid
                  })
# 固定模板
stat_lm2 <- function(mapping = NULL, data = NULL, geom = "line",
                     position = "identity", na.rm = TRUE, show.legend = NA,
                     inherit.aes = TRUE, n = 50, formula = y ~ x,
                     ...){
  layer(stat = StatLm2, data = data, mapping = mapping, geom = geom,
        position = position, show.legend = show.legend, inherit.aes = inherit.aes,
        params = list( n = n, formula = formula, na.rm = na.rm, ...))
}
# 绘图
ggplot(mpg, aes(displ, hwy)) + 
  geom_point() + 
  stat_lm() +
  stat_lm2(formula = y ~ poly(x, 10)) +
  stat_lm2(formula = y ~ poly(x, 10), geom = "point", colour = "red", n =20)

add parameter

我们并不需要显式在图层中包括新的参数，..会将这些参数放到合适的地方。但是你必须在文档中写出哪些参数是可以让用户调整的，以便用户知道他们的存在。举个一个简单的例子

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
#' @export
#' @inheritParams ggplot2::stat_identity
#' @param formula The modelling formula passed to \code{lm}. Should only 
#'   involve \code{y} and \code{x}
#' @param n Number of points used for interpolation.
stat_lm <- function(mapping = NULL, data = NULL, geom = "line",
                    position = "identity", na.rm = FALSE, show.legend = NA, 
                    inherit.aes = TRUE, n = 50, formula = y ~ x, 
                    ...) {
  layer(
    stat = StatLm, data = data, mapping = mapping, geom = geom, 
    position = position, show.legend = show.legend, inherit.aes = inherit.aes,
    params = list(n = n, formula = formula, na.rm = na.rm, ...)
  )
}

上面代码中以#' 开头内容都是roxygon语法，其中@inheritParams ggplot2::stat_identity表示在最后输出的帮助文档中会继承stat_identity的参数说明。而@export则是将函数让用户可见，否则用户无法直接调用。

挑选参数

有些时候，你会发现部分运算是针对所有数据集进行，而非每个分组。比较好的方法就是挑选明智的默认值。例如，我们需要做密度预测，我们有理由为整个图形挑选一个带宽(bandwidth)。下面的”Stat”创建了stat_density()的变体，通过选择每组最优带宽的均值作为所有分组的带宽。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
StatDensityCommon <- ggproto("StatDensityComon", Stat,
                             required_aes = "x",
                             
                             setup_params = function(data, params){
                               if (!is.null(params$bandwidth))
                                 return(params)
                               
                               xs <- split(data$x, data$group)
                               bws <- vapply(xs, bw.nrd0, numeric(1))
                               bw <- mean(bws)
                               message("Picking bandwidth of ", signif(bw,3))
                               
                               params$bandwidth <- bw
                               params
                               },
                             
                             compute_group = function(data, scales, bandwidth = 1){
                               d <- density(data$x, bw = bandwidth)
                               data.frame(x = d$x, y = d$y)
                               
                             }
                             )

stat_density_common <- function(mapping = NULL, data = NULL, geom = "line",
                                position = "identity", na.rm = FALSE,
                                show.legend = NA, inherit.aes = TRUE,
                                bandwidth = NULL, ...){
  layer(stat = StatDensityCommon, data = data, mapping = mapping,
        geom = geom, position = position, show.legend = show.legend,
        inherit.aes = inherit.aes, 
        params = list(bandwidth = bandwidth, na.rm = na.rm, ...))
}

ggplot(mpg, aes(displ, colour = drv)) +
  stat_density_common()

stat density common

作者推推荐用NULL作为默认值。如果你通过自动计算的方式挑选了重要参数，那么建议通过message()的形式告知用户(在答应浮点值参数时，用singif()可以只展示部分小数点)。

变量名和默认图形属性

这部分”stat”会阐述另外一个重要的点。当我们想要让当前”stat”对其他geoms更加有用时，我们应该返回一个变量，称之为”density”而不是”y”。之后，我们可以设置”default_aes”自动地将density映射到y, 这允许用户覆盖它从而使用不同的”geom”.

1
2
3
4
5
6
7
8
9
10
11
StatDensityCommon <- ggproto("StatDentiy2", Stat,
                             required_aes = "x",
                             default_aes = aes(y = stat(density)),
                             
                             compute_group = function(data, scales, bandwidth = 1){
                               d <- density(data$x, bw= bandwidth)
                               data.frame(x = d$x , density=d$y)
                             }
                             )
ggplot(mpg, aes(displ, drv, colour = stat(density))) +
  stat_density_common(bandwidth = 1, geom="point")

stat-area-geom

然而直接在stat中用area geom的结果可能和你想的不同。

1
2
ggplot(mpg, aes(displ, fill = drv)) + 
  stat_density_common(bandwidth = 1, geom = "area", position = "stack")

StatDensity2

密度不是一个相互累加，而是单独计算，因此预测的x没有对齐。我们可以通过在setup_params()计算数据范围的方式解决该问题

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
StatDensityCommon <- ggproto("StatDensityCommon", Stat, 
  required_aes = "x",
  default_aes = aes(y = stat(density)),

  setup_params = function(data, params) {
    min <- min(data$x) - 3 * params$bandwidth
    max <- max(data$x) + 3 * params$bandwidth
    
    list(
      bandwidth = params$bandwidth,
      min = min,
      max = max,
      na.rm = params$na.rm
    )
  },
  
  compute_group = function(data, scales, min, max, bandwidth = 1) {
    d <- density(data$x, bw = bandwidth, from = min, to = max)
    data.frame(x = d$x, density = d$y)
  }  
)

ggplot(mpg, aes(displ, fill = drv)) + 
  stat_density_common(bandwidth = 1, geom = "area", position = "stack")

stat-stack-area

使用”raster”几何形状

1
2
ggplot(mpg, aes(displ, drv, fill = stat(density))) + 
  stat_density_common(bandwidth = 1, geom = "raster")

stat-raster

练习题

  1. 拓展stat_chull，使其能够计算alpha hull, 类似于alphahull. 新的”stat”能够接受alpha做为参数
  2. 修改最终版本的StatDensityComon, 使其能够接受用户定义的min和max. 你需要同时修改layer函数和compute_group()方法
  3. 将StatLm和ggplot2::StatSmooth对比。是什么差异使得StatSmooth比StatLm更加复杂。

创建新的geom

相对于创建新的”stat”, 创建新的”geom”会将难一些，因为这需要你懂得一些grid知识。因为ggplot2基于grid，所以你得要学一些用grid绘图的知识。如果你真的打算学习如何新增一个新的”geom”，Hadley推荐你买Paul Murrell所著的R绘图系统。里面介绍所有和用”grid”绘图相关的知识。

一个简单的geom

让我们先从一个简单的案例入手，尝试实现一个类似于geom_point()的图层

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
# ggproto原型
GeomSimplePoint <- ggproto("GeomSimplePoint", Geom,
                           required_aes = c("x","y"),
                           default_aes = aes(shape = 19, size = 0.1, colour = "black"),
                           draw_key = draw_key_point,
                           
                           draw_panel = function(data, panel_params, coord){
                             coords <- coord$transform(data, panel_params)
                             grid::pointsGrob(
                               coords$x, coords$y,
                               pch = coords$shape,
                               size = unit(coords$size, "char"),
                               gp = grid::gpar(col=coords$colour)
                             )
                             
                           }
)
# 图层函数
geom_simple_point <- function(mapping = NULL, data = NULL, stat = "identity",
                              position = "identity", na.rm = FALSE, show.legend = NA,
                              inherit.aes = TRUE, ...){
  
  layer(geom = GeomSimplePoint, mapping = mapping, data = data, stat = stat,
        position = position, show.legend = show.legend, inherit.aes = inherit.aes,
        params = list(na.rm = na.rm, ...))
}                       

ggplot(mpg, aes(displ, hwy)) + 
  geom_simple_point()

geom-simple-point

上面的代码和构建新的”stat”非常的相似，我们同样需要为4块内容提供属性/方法

  • required_aes: 用户所必需的提供的美术属性
  • default_aes: 默认的图形属性值
  • draw_key: 提供在图例(legend)绘制关键信息的函数，可用?draw_key查看帮助文档
  • draw_panel: 这里就是见证奇迹的地方。该函数接受三个参数作为输入，返回一个grid的”grob”对象。它在每个面板(panel)运行一次。由于它是最复杂的内容，因此我们有必要详细地介绍它。

draw_panel有三个参数

  • data: 数据框，每一列都是一个图形属性
  • panel_params: 一个列表，里面是coord产生的每个面板的参数。你需要将其当做一个不透明的数据结构: 不要看里面的细节，只要将其传递给coord方法。
  • coord: 一个描述坐标系统的对象

你需要共同使用panel_params和coord才能对数据进行转换，即coords <- coord$transform(data, panel_params)。这会创建一个数据框，里面的位置变量会被缩放到0-1之间。得到缩放数据用于调用”grid”的grob函数。(非笛卡尔坐标系统的数据转换比较复杂，你最好是将数据转成已有ggplot2的”geom”所接受的格式，然后传递)。

分组geoms

上一步我们用到的是draw_panel，也就是为每一行元素创建一个图形元素，比如说上面的GeomSimplePoint就是每一行一个点，这是最常见的情况。当然，如果你想为每一个分组绘制一个图形元素，那么我们应该使用draw_group()。

我们用一个简化版的GeomPolygon为例讲解这个知识点:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
GeomSimplePolygon <- ggproto("GeomPolygon", Geom,
                             required_aes = c("x", "y"),
                             default_aes = aes(
                               colour = NA, fill = "grey20", size = 0.5,
                               linetype = 1, alpha = 1
                             ),
                             
                             draw_key = draw_key_polygon,
                             
                             draw_group = function(data, panel_params, coord){
                               n <- nrow(data)
                               if (n <= 2) return(grid::nullGrob())
                               
                               coords <- coord$transform(data, panel_params)
                               
                               first_row <- coords[1, , drop = FALSE]
                               
                               grid::polygonGrob(
                                 coords$x, coords$y,
                                 default.units = "native",
                                 gp = grid::gpar(
                                   col = first_row$colour,
                                   fill = scales::alpha(first_row$fill, first_row$alpha),
                                   lwd = first_row$size * .pt,
                                   lty = first_row$linetype
                                 )
                               )                                                          
                             }
                            )

geom_simple_polygon <- function(mapping = NULL, data = NULL, stat = "chull",
                                position = "identity", na.rm = FALSE, 
                                show.legend = NA, inherit.aes = TRUE, ...){
  layer(
    geom = GeomSimplePolygon, mapping = mapping, data = data, stat = stat,
    position = position, show.legend = show.legend, inherit.aes = inherit.aes,
    params = list(na.rm = na.rm, ...)
  )
  
}

ggplot(mpg, aes(displ, hwy)) +
  geom_point() + 
  geom_simple_polygon(aes(colour = class), fill = NA)

这里有几个注意点

  • 我们重写了draw_group()而不是draw_panel(), 这是因为我们希望polygon是按照绘制，而不是按行绘制。
  • 我们分组数据中不到两行，也就是没有足够的数据点去绘制polygon，因此我们返回了一个nullGrob()。你认为认为这是图形上的NULL: 这是一个grob对象，什么也不画，并且也不占任何空间
  • 在单位上，x和y都应该是native的单位。(默认pointGrob()的单位就是native，因此我这里没有做修改)。多边形线的宽度(lwd)取决于点的大小，而ggplot2计算的点大小返回的mm单位结果，因此作者将其和.pt相乘，将其调整为内部lwd接受的输入。

如果你将我们写的和实际的GeomPolygon比较，你会发现后者重写了draw_panel()，这是因为他用了一些小技巧来创建polygonGrob()从而在一次运行中得到多个polygon。这虽然更加复杂，但是在性能上更优秀。

Collective geoms

从已有的Geom中继承

有些时候，你只想对已有的图层做一些小的修改。在这种情况下，除了从Geom继承以外，你还可以从已有的子类中继承。举个例子，我们可能想要更改GeomPolygon的默认值，使其更好的在StatChull中工作:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
#GeomPolygonHollow继承自GeomPolygon
GeomPolygonHollow <- ggproto("GeomPolyHollwo", GeomPolygon,
                             default_aes = aes(colour = "black", fill = NA,
                                               size = 0.5, linetype = 1, 
                                               alpha = NA))
# layer的stat来自于创建新的stat定义的StatChull
geom_chull <- function(mapping = NULL, data = NULL,
                       position = "identity", na.rm = FALSE, show.legend = NA,
                       inheirt.aes = TRUE, ...){
  layer(stat = StatChull, geom = GeomPolygonHollow, data = data, mapping = mapping,
        position = position, show.legend = show.legend, inherit.aes = inheirt.aes,
        params = list(na.rm = na.rm, ...))
}


ggplot(mpg, aes(displ, hwy)) + 
  geom_point() + 
  geom_chull()

尽管最终的geom_chull不允许你用更改”stat”对应的”geom”, 但是在当前的情况下，凸壳最应该用的”geom”应该就是多边形。

inherit-from-existed-geom

练习题

  1. 比较GeomPoint和GeomSimplePoint
  2. 比较GeomPolygon和GeomSimplePolygon

创建你自己的主题

如果你需要创建自己的完整主题，有以下几个件事情你需要知道

  • 重写已有的元素，而不是修改他们
  • 有四个全局性元素影响几乎所有其他主题元素
  • 完整和不完整元素的比较

重写元素

默认情况下，当你新增一个主题元素，它会从一个已有主题中继承参数值。例如，如下的代码设置key颜色是红色，但它继承了已有的fill颜色。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
theme_grey()$legend.key
#> List of 5
#>  $ fill         : chr "grey95"
#>  $ colour       : chr "white"
#>  $ size         : NULL
#>  $ linetype     : NULL
#>  $ inherit.blank: logi TRUE
#>  - attr(*, "class")= chr [1:2] "element_rect" "element"

new_theme <- theme_grey() + theme(legend.key = element_rect(colour = "red"))
new_theme$legend.key
#> List of 5
#>  $ fill         : chr "grey95"
#>  $ colour       : chr "red"
#>  $ size         : NULL
#>  $ linetype     : NULL
#>  $ inherit.blank: logi FALSE
#>  - attr(*, "class")= chr [1:2] "element_rect" "element"

为了将其彻底重写，使用%+replace%而不是+:

1
2
3
4
5
6
7
8
9
new_theme <- theme_grey() %+replace% theme(legend.key = element_rect(colour = "red"))
new_theme$legend.key
#> List of 5
#>  $ fill         : NULL
#>  $ colour       : chr "red"
#>  $ size         : NULL
#>  $ linetype     : NULL
#>  $ inherit.blank: logi FALSE
#>  - attr(*, "class")= chr [1:2] "element_rect" "element"

全局元素

有四个元素会影响绘图中的全局表现

Element Theme function Description
line element_line() all line elements
rect element_rect() all rectangular elements
text element_text() all text
title element_text() all text in title elements (plot, axes & legend)

很多特殊设置继承下来的属性都可以被以上这四个属性所修改。这对于修改整体背景颜色和总体字体非常有用。

1
2
3
4
5
6
df <- data.frame(x = 1:3, y = 1:3)
base <- ggplot(df, aes(x, y)) + 
  geom_point() + 
  theme_minimal()

base

base
修改整体的字体颜色（不包括坐标）

1
base + theme(text = element_text(colour = "red"))

base-text-red

建议在创建主题的起步阶段，先从修改这些值开始。

完整和不完整的比较

你需要理解完整主题对象和不完整主题对象之间的区别。一个完整的主题对象，就是一个主题函数中设置了complete=TRUE。

以theme_grey()和theme_bw()为例，他们就是完整的主题对象。而调用theme()则会得到一个不完整的主题对象。这两个区别在于，前者是对整体的修改，而后者只是修改了部分的元素。

1
2
3
4
attr(theme_grey(), "complete")
# [1] TRUE
attr(theme(), "complete")
# [1] FALSE

如果在一个完整对象上加上一个不完整对象，那么结果是一个完整对象

1
2
3
theme_test <- theme_grey() + theme()
attr(theme_test(), "complete")
# [1] TRUE

完整主题和不完整主题在添加到ggplot对象上有一些差别

  • 在当前主题对象上增加一个不完整的主题对象，只会修改在theme()中定义的元素。
  • 而在当前主题对象上增加一个完整主题对象，则会将已有主题完全覆盖成新的主题。

参考资料

R

R语言文件读取的基础知识

在使用Python读取文件的时候，通常我需要用到open()打开一个文件，接着是readline()或readlines()读取文件内信息，最后用close()关闭文件。但是在用R语言读取数据时，我就直接用到比较高层的函数，例如read.table或者data.table::fread。R语言也有比较底层的输入输出(I/O)函数，了解这些知识，可以让我们更好的读取文件。

连接(Connection)

连接(connect)是R语言中用于多种I/O操作的基本机制。R提供了9种函数用于为不同输入建立连接，file/url/gzfile/bzfile/xzfile/unz/pipe/fifo/socketConnection，基本上从名字上就能看出这些函数的使用场景。

举个例子，如果想要爬取一个网页，我们可以用url构建目标站点的连接

1
2
3
url_con <- url(description="http://xuzhougeng.top", open="r", encoding = "UTF-8")
class(url_con)
[1] "url"        "connection"

该url_con就可以作为read.table/readLines的文件参数。

上述的description就是指用于将要被建立的连接，可以是一个文件名或URL地址，例子中用的就是一个URL。也可以是其他类型输入。

如果需要在命令行中使用R脚本，让R读取管道的输入，就是”stdin”

1
std_in <- file("stdin", "r")

如果要读取剪切板的内容，则是”clipboard”

1
clip_text <- file("clipboard", "r")

其最终返回一个”connection”对象。

“特殊”连接

这里的”特殊”连接指的是不怎么用的到建立连接的方式，可能在某些特殊的情况下使用。

场景一: 从控制台(console)里获取输入

我们可以用stdin()来从console中获取输入

1
2
3
4
5
> a <- readLines(stdin(), n=2)
hello
world
> a
[1] "hello" "world"

这里的n=2表示读取两行的输入。该方法不能用于shell命令行读取管道的输入

与stdin()对应的还有stdout()，stderr()和nulffile()，这几个命令很少使用。

场景二: 从复制的文本中获取输入

使用clipboard有些时候不清楚自己到底会从剪切板里读取什么样的输入，而且Linux命令行未必有X11支持，所以更好的方法是能够将文本粘贴到脚本中，然后进行读取。

为了解决这一问题，我们需要借助于textConnection将字符串中转成连接。

1
2
3
4
5
6
7
8
9
10
11
12
13
tmp <- "a   b   c   d
1   1   1   B
2   1   2   C
3   1   3   C
4   1   4   C
5   1   5   B
6   1   6   B
7   1   7   C
8   1   8   A
9   1   9   B
10  1   10  C
"
df <- read.table(textConnection(object = tmp),sep="\t", header = T)

和textConnection类似的，还有一个rawConnection用于构建raw对象的连接

读取/输出函数

构建的连接取决于模式可以进行读取和输出操作。R自带函数中和读取/输出相关的可以分成两类

文本模式: readLines, writeLines, cat, sink, scan, parse, read.dcf, dput, dump

二进制模式: readBin, readChar, writeBin, writeChar, load, save

连接管理

对于连接对象，除了读写以外，还有一些列的函数用于查看和管理连接。

  • open: 打开连接
  • close(): 关闭连接
  • seek(): 在连接中进行跳转
  • isOpen判断连接是否打开
  • isIncomplete判断连接的读取最后一次是否完整， 写出时是否有未输出的内容
  • flush刷新连接的输出流，保证内容都完全写出。用close()安全的关闭连接
  • showConnections: 显示目前的所有连接
  • getConnection: 获取指定连接
  • closeAllConnections: 关闭所有连接

大部分函数都不好找到使用场景，除非你要编写某些数据类型的读取函数。

生信软件工具箱

Hi-C数据可视化-组装角度

Hi-C数据可视化-组装角度

这里讨论HiC的可视化是从组装角度出发，也就是如何展示contig和contig之间的关系

Hi-C数据可视化（我所知）有下面几个

无论是何种分析工具，最核心的步骤就是提供它们需要的输入数据，也就是将fastq转成所需要的输入形式。虽然本质上都是对contact matrix进行展示，但问题在于部分工具定义了其特有的存放格式，需要进行一些格式转换。

基于上述考虑，我在数据预处理上选择使用HiC-Pro，因为HiC-Pro本身就能够将其结果转成fithic,juicebox的输入，此外HiCExplorer提供了hicCovertFormat，能够将Hic-Pro转成cool/mcool/h5/homer/ginteractions等格式。

后续假定HiC-Pro输出的文件为三个文件

  • input_20000.matrix: 原始上三角矩阵
  • input_20000_abs.bed: 坐标信息
  • input_20000_iced.matrix: ICE normalization的矩阵

HiTC

HiTC是R包，能够直接读取HiC-pro的输出结果，对应的函数是importC

1
2
3
4
5
# BiocManager::install("HiTC")
library(HiTC)
## Load Hi-C data
x <- importC("input_20000.matrix", xgi.bed = "input_20000_abs.bed")
show(x)

由于读取这一步时，它会创建NXN个 Map，如果有1,000个contig， 那么最终会得到1,000,000个map，这是一个非常可怕的计算量，R语言会因此而崩溃。因此对于contig级别的基因组而言，不应该用这个工具进行展示。

在可视化方面方面，主要是mapC，它接受一个HTClist作为输入

1
2
3
4
sset <- reduce(x, chr=c("chr5","chr6","chr7"))
x90_500 <- HTClist(mclapply(sset, binningC,
binsize=500000, bin.adjust=FALSE, method="sum", step=1))
mapC(imr90_500)

考虑到矩阵以对角线对称，实际上只需要展示上三角即可

1
2
hox <- extractRegion(x$chr6chr6, chr="chr6", from=50e6, to=58e6)
plot(hox, maxrange=50, col.pos=c("white", "orange", "red", "black"))

更多参考: https://bioconductor.org/packages/release/bioc/vignettes/HiTC/inst/doc/HiC_analysis.pdf

HiCPlotter

HiCPlotter依赖Python2.7和对应的Numpy, Scipy,Matplotlib, 由于两年没有更新，因此GitHub里提到的版本都比较古老，经测试，Numpy=1.15.1，Scipy=1.1.0， Matplotlib=2.2.3 能够正常运行。

它支持HiC-pro的输出，无需做格式转换。

举个例子

1
2
3
4
5
6
7
8
9
10
11
12
13
14
python HiCPlotter.py \
    -f input_20000_iced.matrix \
    -o Exemple \
    -wg 1
    -r 20000 \
    -tri 1 \
    -bed input_20000_abs.bed \
    -n Test -chr chrX -ptr 1
# -r: 分辨率
# -f 输入矩阵
# -tri: HiC-Pro输入要设置为1
# -bed: HiC-Pro输入要指定bed
# -wg: 全基因组的交互
# -ptr: 绘制三角还是全部

HiCPlotter只能单个染色体或者从起始染色体到指定染色体的全基因组互作，并不能绘制给定任意多个染色体的互作关系，功能比较简单。

HiCExplorer可视化

HiCExplorer开发deeptools的团队推出的HiC处理工具，包括格式转换，HiC矩阵质控和可视化。在功能上和性能上，HiCExplorer都比之前提及的工具强，但是如果只用需要的功能，它并不会很复杂。

HiCExplorer要求contact matrix以h5格式存放，因此需要对HiC-Pro的结果进行转换

1
2
3
4
5
hicConvertFormat \
	-m input_20000_iced.matrix \
	-o input_20000 \
	--bedFileHicpro input_20000_abs.bed \
	--inputFormat hicpro --outputFormat h5

输出的是input_20000.h5文件，该h5文件便是后续分析的矩阵输入文件

使用hicPlotMatrix绘制热图展示HiC数据

1
2
3
4
5
6
7
8
9
10
# 绘制单条染色体
hicPlotMatrix --matrix input_20000.h5 -out chrX_chrY \
    --region chrX 
# chrX和chrY的交互
hicPlotMatrix --matrix input_20000.h5 -out chrX_chrY \
    --region chrX \
    --region2 chrY
# 展示指定的2个染色体
hicPlotMatrix --matrix input_20000.h5 -out chrXchrY \
    --chromosomeOrder chrX chrY

更多信息参考: https://hicexplorer.readthedocs.io/en/latest/content/tools/hicPlotMatrix.html

生信软件工具箱

Canu的graph和contig生成过程

这篇文章是我发现canu输出的contig可能存在misassembly(错误组装), 为了探索这种错误是如何产生的，我尝试解决如下问题

  • Canu输出contig的基本步骤
  • contig和GFA是什么关系
  • 如何提取contig对应的read
  • 如何检查contig的graph

contig的构建过程

Canu的流程分为三个步骤，前两步是原始输入的纠错，最后一步是基于纠错后的reads来构建contigs。

步骤三的流程图

Step3

步骤三的输出文件

1
2
3
4
5
6
7
8
9
10
11
12
13
14
0-mercounts
1-overlapper
3-overlapErrorAdjustment
4-unitigger
5-consensus
prefix.ctgStore
prefix.ctgStore.coverageStat.log
prefix.ctgStore.coverageStat.stats
prefix.ovlStore
prefix.ovlStore.config
prefix.ovlStore.config.txt
prefix.ovlStore.per-read.log
prefix.ovlStore.summary
prefix.utgStore

注: prefix表示文件名前缀。

根据流程和输出文件可知这一步骤可以细分为5步

  • 构建read之间的overlap得到ovlStore
  • 分析overlap检测错误
  • 重新计算overlap的alignment
  • 使用bogart构建contig
  • 对contig进行consensus
  • 构建输出, graph, layout, sequences

这里面最核心的是bogart, 它根据overlap情况构建graph，4-unitigger有如下的GFA文件，就是运行时产生的文件。

1
2
3
4
5
6
7
prefix.best.edges.gfa
prefix.initial.assembly.gfa
prefix.final.assembly.gfa
prefix.contigs.gfa
prefix.unitigs.gfa
prefix.unitigs.aligned.gfa
prefix.contigs.aligned.gfa

bogart基于graph寻找最优路径, 最终得到unitig和contig两个结果，其中unitig更加碎一些，但更加准一些（Contigs, split at alternate paths in the graph）。

最终的contig还需要进一步的consensus，得到最终的输出

  • prefix.contigs/unitigs.fasta: Fasta序列
  • prefix.contigs/unitigs.gfa: contig.gfa已经在Canu1.9之后被移除，并非记录Fasta的构建过程
  • prefix.contigs/unitigs.layout: 上述GFA并不存放序列，实际序列在layout中
  • prefix.contigs/unitigs.layout.readToTig: contig ID和read ID的对应关系
  • prefix.contigs/unitigs.layout.tigInfo: 每个contig的信息
  • prefix.unitigs.bed: untig和contig的对应关系

Contig和GFA的关系

对于主目录输出结果中contig/untigs对应的fasta和gfa，我们需要明白其中fasta并非是由gfa生成，而是记录了contig与之前contig的overlap关系。

因此，如果用bandage可视化Canu 1.8之前的GFA文件，实际上你看到的是最终的contig之间的关系，而非read之间的关系。在Canu 1.9的更改日志中写道，

  • Output file ‘contigs.gfa’ was removed because it was misleading

我们实际需要的是canu通过解析read之间的overlap关系图来得到最终的contig的gfa文件。而Canu并没有提供符合我们需求的文件，因为Canu输出的GFA文件中都没有以P为开头的记录，即没有记录contig的路径。

假如我对其中一个contig存在怀疑，那么我应该如何检查该contig对应read的overlap关系呢？

提取contig对应的read信息

注意，不同版本的Canu输出数据库未必兼容，也就是用和建立数据库不同版本的工具会报错。

以我自己的一个数据为例，对于.gfa中的一个编号，例如tig00000007，

1
S       tig00000007     *       LN:i:62235

我们可知他的实际ID是7。它对应的read信息可以在.layout.readToTig找到。

1
2
3
$ awk '$2==7' prefix.contigs.layout.readToTig | head -n 2
632820	7	ungapped	0	37356
56109	7	ungapped	507	38683

第一列即是contig对应的readID，例如632820

利用ovStoreDump我们就可以提取和此read overlap的所有read

1
ovStoreDump -S prefix.seqStore -O unitigging/prefix.ovlStore -picture 632820

ovStoreDump可以以多种形式输出overlap的信息，例如-picture就是以ASCII展示overlap

或者根据readID用sqStoreDumpFASTQ提取实际的序列

1
2
3
sqStoreDumpFASTQ -S prefix.seqStore -o - -fasta -r 632820 -raw | seqkit seq - | head
sqStoreDumpFASTQ -S prefix.seqStore -o - -fasta -r 632820 -corrected | seqkit seq - | head
sqStoreDumpFASTQ -S prefix.seqStore -o - -fasta -r 632820 -trimmed | seqkit seq - | head

-raw/-corrected/-trimmed表示提取不同阶段的结果

输出结果里使用的实际输入序列的ID编号，可以在prefix.seqStore/readNames.txt中找到readID和原始编号的对应关系。

掌握上述的操作，就可以提取指定的contig的read，然后将read回帖到该contig上，利用IGV可视化。

可视化检查contig的GFA

提取contig对应的readID, 然后根据readID抽取gfa

1
2
3
4
5
6
contigID=$1
prefix=$2
readtotig=$prefix.contigs.layout.readToTig
gfa=unitigging/4-unitigger/$prefix.initial.assembly.gfa
awk -v id=$contigID '$2==id {printf("read%08d\n",$1)}' $readtotig > readid.txt
grep -f readid.txt $gfa > $contigID.gfa

输出的.gfa就可以用Bandage进行可视化。但实际上发现这个操作并不太可行，不如直接上一节用IGV可视化alignment来的直接。

GFA可视化

生信软件工具箱

SALSA:基于Hi-C辅助组装长读长组装结果

SALSA最初是在BMC genomics上发表，应该是当时最早提出利用Hi-C对contig进行纠错的软件，随后3D-DNA引入了这一策略。而最近升级之后的SALSA又在PLoS Computational Biology发表，可能是目前第一个用组装软件输出的GFA信息的工具。

使用Hi-C辅助组装，有以下的可能问题

  1. 原先的contig存在错误，scaffold引入了之前的错误
  2. 输入contig过短时，朝向错误引起inversion

SALSA2的优势

  1. 整合了assembly graph
  2. 能够对原始的contig进行纠错

它的工作流程如下：

workflow

SALSA2从PacBio/ONT的组装结果开始，建议提供记录模棱两可重建位置的GFA文件。Hi-C数据回贴到contig上，contig根据Hi-C的覆盖度分析潜在的组装错误，然后被打断。之后根据GFA和Hi-C信息构建hybrid scaffold graph，最后scaffold根据该graph进行重建。

SALSA提出了一种新的染色体划分方式，也就是检查每次迭代后是否存在错误的合并(mis-join)，如果发现大部分的join都是错误的，也意味着目前scaffold划分是最佳的，就可以停止了。

A mis-join detection step is performed after each iteration to check if any of the joins made during this round are incorrect. Incorrect joins are broken and the edges blacklisted during subsequent iterations. This process continues
until the majority of joins made in the prior iteration are incorrect. This provides a natural stopping condition, when accurate Hi-C links have been exhausted.

软件安装

SALSA2的依赖环境如下:

  • Python 2.7
  • Networkx version < 1.2
  • BOOST libraries

通过编译的方法进行安装

1
2
3
git clone https://github.com/marbl/SALSA.git
cd SALSA
make -j8

软件使用

软件要求以下几个输入文件:

  • 三代组装结果: contigs.fasta
  • Hi-C比对后BAM: alignment.bam
  • (可选) GFA文件: contigs_graph.gfa

第一步: 将Hi-C回贴到组装结果中，这一步主要用bowtie2/bwa比对原始的fastq文件，之后对Hi-C的比对结果进行一些预处理。目前已经有一些流程工具可以实现这一步

第二步: 将BAM转成BED，并按照Read Name进行排序。假设上一步的输出是alignment.bam

1
2
bamToBed -i alignment.bam > alignment.bed
sort -k 4 alignment.bed > tmp && mv tmp alignment.bed

如果经HiC-Pro处理的数据，那么在hic_results下有一个.allValidPairs, 可以通过awk进行转换成所需要的bed

1
awk 'BEGIN{OFS="\t"} {print $2,$3,$3+50,$1"/1","60",$4} {print $5,$6,$6+50,$1"/2","60",$7}' input.allValidPairs > alignment.bed

注意: 关于这个alignment.bed文件，经过我的测试，只需要按照Read Names排序即可，并不需要保证两行两行的配对，SALSA会预先处理。

第三步: 构建contig长度文件

1
samtools faidx contigs.fasta

第四步: 根据需求挑选命令运行. 注意-e {Your Enzyme}需要根据实际情况修改，如果是MboI，就写-e GATC.

需求1: 不需要纠错，最基本的组装

1
python run_pipeline.py -a contigs.fasta -l contigs.fasta.fai -b alignment.bed -e {Your Enzyme} -o scaffolds

需求2: 使用Hi-C对输入assembly进行纠错(加入参数-m yes)

1
python run_pipeline.py -a contigs.fasta -l contigs.fasta.fai -b alignment.bed -e {Your Enzyme} -o scaffolds -m yes

断裂位点位置存放在输出的input_breaks

注意: 关于纠错，建议阅读下assembly error detection isn’t much efficient on my data，这里面在讨论了，Hi-C信号选择对结果的影响。我也在实践中发现使用过滤后的信号会影响最终的breaks位点的数目。

需求3: 使用GFA文件辅助组装(加入了参数-g contigs_graph.gfa)

1
python run_pipeline.py -a contigs.fasta -l contigs.fasta.fai -b alignment.bed -e {Your Enzyme} -o scaffolds -m yes -g contigs_graph.gfa

需求4: 使用更加准的unitigs和unitigs GFA. 组装软件，例如Canu会输出xxx.unitigs.fasta和xxx.unitigs.gfa, 里面的序列相对于contig更短，但是错误更少，在某些情况下更适合用于scaffold。

1
python run_pipeline.py -a unitigs.fasta -l unitigs.fasta.fai -b alignment.bed -e {Your Enzyme} -o scaffolds -m yes -g unitigs_graph.gfa

最终在-o指定的目录下输出结果，其中scaffolds_FINAL.fasta和scaffolds_FINAL.agp是两个最主要的结果。

如果你因为用了unitigs导致scaffolds中gap太多，那么SALSA可以根据contigs.fasta来关闭gap

1
python stitch.py -c contigs.fasta -u unitigs.fasta -p scaffolds_iteration_x.p -o output_scaffolds.fasta

其中scaffolds_iteration_x.p里的x取决于你的迭代次数。

假如你需要在Juicebox中可视化结果，SALSA提供了conver.sh用于数据转换。你需要先从https://github.com/aidenlab/juicer/wiki/Download下载Juicer tools jarw.

1
wget https://s3.amazonaws.com/hicfiles.tc4ga.com/public/juicer/juicer_tools_1.14.08.jar

然后修改SALSA的convert.sh中的JUICER_JAR为你实际路径，最后运行

1
bash /path/to/SALSA/convert.sh SALSA_OUTDIR

输出文件是salsa_scaffolds.hic

SALSA的所有参数说明如下:

  • -a/–assembly: 初始的组装结果，header中不能有:
  • -l/–length: contigs的长度
  • -b/–bed: BAM文件转换后的bed
  • -o/–output: 输出文件夹
  • -c/–cutoff: 用于scaffold的contig最短长度
  • -g/–gfa: GFA文件
  • -e/–enzyme: 限制性内切酶
  • -i/–iter: 迭代数
  • -x/–dup: 存放重复的contig信息
  • -s/–exp: 期望的基因组大小
  • -m/–clean: 是否检查Mis-assembly

不客观的测评

SALSA这个工具无论是安装还是使用上，都非常的简洁易懂。我在使用过程中就没有遇到过安装错误，并且运行效率高（基于Python和C++）。最终的组装结果可能和文章说的那样，得到相对最优的scaffolds，无法得到染色体级别的输出结果。

Building on this work, SALSA2 does not require manual parameter tuning and is able to utilize all the contact information from the Hi-C data to generate near optimal sized scaffolds permitted by the data using a novel iterative scaffolding method

SALSA的结果的确是比较可靠的，在contig的朝向上出错更少。不过，对于mis-assembly的鉴定，目前SALSA可能有些问题，在我实践中发现有些明显的错误它没有找到。

此外，在我的几个项目实战中，我发现3D-DNA能够得到染色体水平的结果的情况，而SALSA却无法得到相同的结果。或许是输入的HiC信号不同，之后会尝试用Juicer的输出作为SALSA的输入信号，进行比较。

参考资料

其他我写过的HiC工具