luria的个人博客分享 http://blog.sciencenet.cn/u/luria

博文

HiCUP的使用以及结果解读

已有 1273 次阅读 2019-1-30 14:30 |个人分类:Hi-C|系统分类:科研笔记| Hi-C, 3D genome, assembly, proximity-ligation, Hi-C辅助组装, assembly, assembly

1. HiCUP安装

HiCUP需依赖PerlBowtie/Bowtie2R以及Samtools

在官网http://www.bioinformatics.babraham.ac.uk/projects/hicup/下载最新版的HiCUP,例如目前的最新版是v0.6.1

wget http://www.bioinformatics.babraham.ac.uk/projects/download.html#hicup

# 下载到压缩包hicup_v0.6.1.tar.gz

tar -zxvf  hicup_v0.6.1.tar.gz

# HiCUPPerl写的解压即可使用,无需编译

 

HiCUP中的比对工具是Bowtie/Bowtie2,建议采用最新版的Bowtie2。和HiCUP一样,Bowtie2Linux二进制版,解压即可使用,无需编译

 

2. HiCUP的比对策略

HiCUP6Perl脚本组成,分别如下:

(1) HiCUP Digester:确定reference上的酶切位点

(2) HiCUP:为主程序,依次执行以下步骤

(3) HiCUP Truncater:在reads上寻找酶切位点,并将reads切开

(4) HiCUP Mapper:将reads比对到参考基因组上,如果输入的是PEreads,则R1R2分开单独比对到reference上。比对内部调用bowtiebowtie2比对。这一步会利用到事先建好的bowtie2 index

(5) HiCUP Filter:结合HiCUP Digester生成的酶切位点文件,过滤掉常见的Hi-C artefacts,例如Dangling Ends

(6) HiCUP Deduplicator:移除(仅保留一处最佳比对) PCR重复

 

3. HiCUP的使用

3.1 先采用bowtie2-buildreference建立索引

/path/to/bowtie2-build reference.fasta reference

 

3.2 采用hicup目录下的hicup_digesterreference上寻找酶切位点,生成酶切信息文件

/path/to/hicup_v0.6.1/hicup_digester \

     --re1 ^GATC,MboI \

     --genome reference \

     --outdir /path/to/output/dir \

     reference.fasta

 

3.3 再采用hicup主程序进行分析

有两种方式运行hicup,其一是将所有参数写到config文件中,可以先运行

hicup --example


生成样例
config文件,修改其中的参数,并运行

或者直接用命令行运行,如下:

/path/to/hicup_v0.6.1/hicup \

     --bowtie2 /path/to/bowtie2 \

     --digest Digest_reference_MboI_None_18-36-52_07-08-2018.txt \

     --format Sanger \

     --index /path/to/reference/index \

     --keep \

     --outdir /path/to/output/dir  \

     --threads 40 \

     /path/to/reads*.fastq.gz

 

4. 结果解读

最重要的两个文件是:

Ø  xx_R1_2.hicup.sam

Ø  xx_ R1_2.HiCUP_summary_report.html

前者是最终的sam文件,后者是全程汇总报告

 

每步结果具体如下:

HiCUP运行后hicup_truncater先生成xx.R1.trunc.fastqxx.R2.trunc.fastq文件,它是按照酶切位点对Reads进行切除,因此得到的reads长短不一。完成后会统计Truncated reads数,Not Truncated reads数,以两者的占比,总的reads数,平均reads长度,并记录在hicup_truncater_summary_xx.txt文件中。两个fastq对应的svg图像中仅绘制了前两项的条形图。

 

hicup_mapper调用bowtie2-align-s进行比对处理,生成xx_ R1_2.pair.sam文件,同时也统计了太短而无法比对的reads数,唯一比对的reads数,多处比对的reads数,比对不上的reads数,成对的reads数,以及以上各项的占比,记录在hicup_mapper_summary_xx.txt文件中,顺带还画了R1R2的比对结果条形图(xx_R*_trunc.fastq.mapper_barchart.svg)

 

hicup_filter对上步生成的sam文件进行过滤,分别将每种invalid ditag类型(包括same internalsame dangling endssame circularisedre_ligation及其它)过滤掉的比对结果写入到hicup_filter_ditag_rejects_xx/目录下

过滤后可用于下步分析的sam文件为xx_R1_2.filt.sam

当然,对结果ditag结果也进行了统计,包括valid pairsinvalid pairssame circularisedsame fragment dangling ends,same fragment internal, re-ligation, contiguous sequence, wrong size以及total pairs等类型,记录在hicup_filter_summary_xx.txt文件中。同时也画了Di-tag长度分布svg图和各种类型的svg piechart图。

hicup_deduplicator去重,并且计算了序列内的unique reads(<10kb>10kb分别统计)以及序列间的unique reads,并画了相应的svg

最终的结果为xx_R1_2.hicup.sam

注意,这个最终的sam文件中仅存放De-duplicationUnique Di-TagsRead Pairs,即如下图中的read pairs(强调:是read pairs!)

最后不仅给出了所有步骤reads数据变化的汇总(HiCUP_summary_report_xx.txt),同时还给给出了非常美观的html报告(xx_ R1_2.HiCUP_summary_report.html)

 

另外,需要强调一点,hicup的结果如果想接到Lachesis做组装。在hicup生成的sambam的时候一定要用sort -n,按名字排序!如果用默认的参考序列位置排序,则Lacheis中会生成Cluster信息,但是没有Order信息!

 

运行效率:

PE reads总数据量45G,参考基因组440M。运行指定40个核,估算出的CPU最大消耗35个整,最大消耗内存6.34G,运行时间约7.5h,最大消耗是hicup_mapper比对。



 




http://blog.sciencenet.cn/blog-2970729-1159899.html

上一篇:采用Stacks2分析RADseq数据[无参]
下一篇:MUMmer4的安装与使用

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备14006957 )

GMT+8, 2019-10-21 13:03

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部