NC结果复现 | 可视化R包endplot开发测试

endplot算是自己对END-seq数据分析工作的总结和理解。END-seq并不像single-cell广为人知，网上可用的资源满天飞。学习新技术阅读一手的文献必不可少，理解其精华并消化，如果能够联系作者答疑解惑就更好了。
作为应用者最主要的目的自然是分析自有数据，谁愿意把时间花在深入了解技术细节上，通常也就是跟着教程走一遭，跑一遍自己的数据。
然后，END-seq数据的下游分析并没有完善的流程，每次分析时都挺折腾，将重复且复杂的事情简单化本身也挺有意义。故而，将分析数据的代码整理成endplot包，以便后续使用。

schematic of END-seq break pattern

endplot主要用于END-seq数据的下游可视化，即用bam和peak作为输入文件。数据处理过程参考文章<ATM and PRDM9 regulate SPO11-bound recombination intermediates during meiosis>。

library(endplot)

bam <- c(wt='wt_sort.bam', atm='atm_sort.bam')
peak <- c(wt='wt_peaks.bed', atm='atm_peaks.bed')
obj <- create_obj(bam, peak, cpu=2)
obj
An object of class endobj:
8421 peaks across 2 samples
flank size of DBS: 6000

输入为两个命名向量，分别为bam和peak文件，元素名作为样本名。读取多个样本时指定cpu参数，可以多样本并行读取加快速度。数据读取后会生成一个endobj对象，后续所有绘图都基于此对象。

str(obj)
Formal class 'endobj' [package "endplot"] with 5 slots
  ..@ mat_sense    : num [1:8421, 1:12001] 0.0373 0 0 0 0 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:8421] "chr1_4604329_-:wt" "chr1_4917816_+:wt" "chr1_5021486_-:wt" "chr1_5151136_+:wt" ...
  .. .. ..$ : NULL
  ..@ mat_antisense: num [1:8421, 1:12001] 0 0 0 0 0 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:8421] "chr1_4604329_-:wt" "chr1_4917816_+:wt" "chr1_5021486_-:wt" "chr1_5151136_+:wt" ...
  .. .. ..$ : NULL
  ..@ peak         :'data.frame':       8421 obs. of  7 variables:
  .. ..$ chr   : chr [1:8421] "chr1" "chr1" "chr1" "chr1" ...
  .. ..$ start : int [1:8421] 4609192 4924684 5026434 5158108 5524611 9942224 10090345 10161820 10453548 10615561 ...
  .. ..$ end   : int [1:8421] 4609192 4924684 5026434 5158108 5524611 9942224 10090345 10161820 10453548 10615561 ...
  .. ..$ peak  : chr [1:8421] "chr1_4603201_4617457" "chr1_4916882_4930750" "chr1_5020423_5034550" "chr1_5150395_5163878" ...
  .. ..$ count : int [1:8421] 7 10 18 5 30 7 43 14 7 24 ...
  .. ..$ strand: chr [1:8421] "-" "+" "-" "+" ...
  .. ..$ sample: chr [1:8421] "wt" "wt" "wt" "wt" ...
  ..@ flank        : num 6000
  ..@ rpm_sfactor  : Named num [1:2] 26.8 48
  .. ..- attr(*, "names")= chr [1:2] "wt" "atm"

为了查看和操控数据便捷性，实现了一些常用的方法，例如，下面的查看metadata和取子集：

head(obj)
                   chr   start     end                 peak count strand sample
chr1_4604329_-:wt chr1 4609192 4609192 chr1_4603201_4617457     7      -     wt
chr1_4917816_+:wt chr1 4924684 4924684 chr1_4916882_4930750    10      +     wt
chr1_5021486_-:wt chr1 5026434 5026434 chr1_5020423_5034550    18      -     wt
chr1_5151136_+:wt chr1 5158108 5158108 chr1_5150395_5163878     5      +     wt
chr1_5519544_-:wt chr1 5524611 5524611 chr1_5518208_5532881    30      -     wt
chr1_9936865_-:wt chr1 9942224 9942224 chr1_9935931_9949799     7      -     wt

atm <- obj[obj$sample == 'atm',]
atm
An object of class endobj:
4261 peaks across 1 samples
flank size of DBS: 6000

为了验证数据处理过程是否正确，下载了文献里面的两个数据：WT和ATM。关于ATM的功能这里就不在介绍，感兴趣可以阅读文章。这里主要来重现文章的结果，验证数据处理过程是否正确。

Fig.4B

p <-  plot_longResection(obj)
p

因为没有下载SSDS数据，所以无法重现文章的结果。这里展示一下可以利用函数很方便地绘制long-resection、short-resection的这类图。

Fig.4H复现

p <- plot_frequency(obj, cols=c('blue','green',"orange2","purple4"))

文章结果：

从结果对比来看，整体pattern一致，而且atm的 minimum resection和maximum resection均值也比较接近。

Fig.5A复现

p <- plot_curve(obj, flank=5000, cols=c('red', 'black'))
p

文章结果：

从上面两个结果的对比可知，虽然细节上稍微有一些差异，但整体特征基本一致。

Fig.5B复现

plot_heat(obj, sname='wt', binsize=20)

文章结果：

绘制其他图形使用的是ggplot2，而热图使用的pheatmap。从上面的结果可知，复现与文章的结果整体特征一致。

Fig.8C复现

文章结果：

从上面的结果来看，正负链峰值的间距完全一致，曲线的pattern也几乎一致。

NC结果复现 | 可视化R包endplot开发测试