【preseq】 评估文库复杂性

问题：

前面看到这个图，描述随着测序深度逐渐增加，Repetition rate 逐渐增加。觉得自己进行library size 抽样，还是比较麻烦，samtools flagstat 计算distinct reads 数目也花时间。
看到提到使用preseq 来评估文库复杂性，于是来尝试一下。

https://www.nature.com/articles/ncomms7033

简单介绍图片内容：

• 左侧图a 展示ULI-NChIP 实验方法与NChIP / Low -input ChIP-seq 差别（标黄色的框框，是不同的step）
• 右侧是使用ULI-NChIP对三种组蛋白进行文库复杂性评估，黑色线条是对照组，可看出H3k4me3 实验相比H3K27me3 等等，出现的重复序列更多一些（看红色线条和黑色线条高度差），可能因为同一个样本，捕获的H3K4me3 DNA 量比H3K27me3 要低一些，测序之前需要多进行几轮PCR引起的。

---

安装

• 存在C++ /R 语言版本，用conda 可以安装。
  https://bioconda.github.io/recipes/preseq/README.html
  https://github.com/smithlabcode/preseq

conda install preseq

主要有两个子命令。lc_extrap 和 c_curve，支持的输入文件包括bam 和bed 文件

• 输入文件需要排序
  Input files can be either in BED or BAM file format. The file should be sorted by chromosome, start position, strand position, and finally strand if in BED format. If the file is in BAM format, then the file should be sorted using bamtools or samtools sort.

• 帮助： 常用的参数：
  -o；设置输出目录；;
  -s：library size 步长；
  -B ：默认是bed 文件;
  -P : 默认是单端；

$ preseq c_curve 
Usage: c_curve [OPTIONS] <sorted-bed-file>

Options:
  -o, -output   yield output file (default: stdout) 
  -s, -step     step size in extrapolations (default: 1e+06) 
  -v, -verbose  print more information 
  -P, -pe       input is paired end read file 
  -H, -hist     input is a text file containing the observed histogram 
  -V, -vals     input is a text file containing only the observed counts 
  -B, -bam      input is in BAM format 
  -l, -seg_len  maximum segment length when merging paired end bam reads 
                (default: 5000) 
  -r, -seed     seed for random number generator 

Help options:
  -?, -help     print this help message 
      -about    print about message

测试：

• 计算文库大小与distinct read 关系

bamToBed -i trim_25.uniq.bam >./try/trim_25.uniq.bed
preseq c_curve -s 1e+03 -o estimates.txt trim_25.uniq.bed

estimates.txt.png

• 画图

library(tidyverse)
data <- read_delim("estimates.txt",delim="\t")
ggplot(data,aes(x=total_reads,y=distinct_reads)) + geom_line()

image.png

---

思考：

• lc_extrap 功能和c_curve 应该是类似的，测试时候提示使用defect mode 不清楚ERROR是什么原因,加上了-D参数，运行正常。

(R_env) [17:34:44] kcao@localhost:~/hexinyi/kcao/try
$ preseq lc_extrap  -o yield_estimates.txt trim_25.uniq.bed 
ERROR:  too many defects in the approximation, consider running in defect mode
(R_env) [17:35:25] kcao@localhost:~/hexinyi/kcao/try
$ preseq lc_extrap  -D -o yield_estimates.txt trim_25.uniq.bed 

yield_estimates.txt.png

欢迎评论交流~😀