Juicer实战详解

2019-12-20 13:32:15 浏览数 (1)

Juicer软件的运行是非常简单的,只需要设置几个参数就可以了,本文利用官网的小的测试测试数据集来展示该软件的基本用法。

1. 下载测试数据

从以下链接下载测试数据集

https://github.com/aidenlab/juicer/wiki/Running-Juicer-on-a-cluster

这里选用的是红框标记的小的测试数据集,如果想要体验完整的分析功能,可以option1提供的测试数据

代码语言:javascript复制
wget http://juicerawsmirror.s3.amazonaws.com/opt/juicer/work/HIC003/fastq/HIC003_S2_L001_R1_001.fastq.gz
wget http://juicerawsmirror.s3.amazonaws.com/opt/juicer/work/HIC003/fastq/HIC003_S2_L001_R2_001.fastq.gz

样本的原始序列放置在软件安装目录的work/sample/fastq目录下,sample替换成自己定义的名称。

2. 运行

这里我没有下载官方提供的参考基因组,而是采用了UCSC下载的基因组。对于自己下载的参考基因组,首先建立bwa的索引,为了方便管理,统一将基因组序列和索引文件放在软件安装目录的references文件夹下,用法如下

代码语言:javascript复制
cd references
wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/hg19.fa.gz
gunzip hg19.fa.gz
bwa index hg19.fa hg19.fa

其次建立酶切图谱,放置在restriction_sites目录下,用法如下

代码语言:javascript复制
python  misc/generate_site_positions.py HindIII hg19 references/hg19.fa

第一个参数根据实际使用的内切酶来选择,酶切图谱生成之后,可以在输出文件的基础上,生成染色体大小文件, 用法如下

代码语言:javascript复制
awk 'BEGIN{OFS="t"}{print $1, $NF}'  hg19_HindIII.txt > hg19.chrom.sizes

其实也可以从UCSC直接下载物种对应的染色质长度文件,对于其他来源的基因组文件,用上述方式更加通用。hg19.chrom.sizes文件的内容如下

代码语言:javascript复制
chr1    249250621
chr2    243199373
chr3    198022430
chr4    191154276

该文件决定了最终的Hi-C图谱包含的染色体名称,对于一些random, unplace_scaffold序列,可以直接在该文件中去除,这样在不会出现在最终结果中。 准备好样本的原始序列和参考基因组的文件之后,就可以运行juicer了。用法如下

代码语言:javascript复制
juicer.sh 
-z references/hg19.fa 
-p restriction_sites/hg19.chrom.sizes 
-y restriction_sites/hg19_HindIII.txt 
-d /home/pub/software/juicer/work/HIC003/ 
-D /home/pub/software/juicer 
-t 5

-z参数指定参考基因组fasta所在路径,在该路径下必须同时存在对应的bwa索引;-p参数指定染色体长度文件;-y指定基因组酶切图谱的路径;-d指定样本原始文件存放的路径;-D指定软件的安装路径,-t指定bwa比对使用的线程数,默认是使用全部线程。

需要注意的是, 在指定文件路径时,最好指定成绝对路径,特别是fastq文件所在路径。因为软件运行过程中会使用软链接,相对路径会出错。

软件运行完成之后,在样本对应的目录下,会生成以下目录

  1. splits
  2. aligned

splits目录下存放的是中间结果,由于hi-C数据量很大,所以会将原始序列拆分成很多份,并行运算,加快速度。默认每份包含22.5M的reads, 当然这个可以通过-C参数调整,该参数指定拆分文件的行数,默认是90000000, 注意fastq文件4行代表一条序列,所以这个参数的值必须是4的倍数。拆分后序列的R1和R2端分别通过bwa比对基因组,然后合并,筛选嵌合体序列,去重复,生成预处理后的结果文件。

aligned目录下存放的是最终结果,包含了可以导入juicebox的后缀为hic的图谱文件, inter.hicinter_30.hic, 30表示通过MAPQ > 30进行过滤之后的结果。完整流程还会进行后续处理,包括识别TAD, 染色质环等结构。其中识别染色质环的HICCUPs算法必须通过GPU加速运行才可以,所以没有安装GPU卡的普通服务器无法运行这个步骤。

从上述过程可以看到,juicer的使用确实非常简单。由于Hi-C数据的测序量非常大,以及后续分析算法的复杂度,对服务器计算资源的要求相当高,必须高性能服务器才能满足要求,而该软件所需的GPU卡成本也非常高,一块的成本在2万元左右,这些因素一定程度制约了Hi-C的普及和发展。

0 人点赞