通过第一讲:三维基因组学习笔记,我们了解了3D基因组研究范围,然后根据我在生信技能树发布的生信工程师标准提炼出基础技能,也就是第二讲:生信基础技能 。最后提炼出了数据分析流程,并且安装好了对应的软件,也就是第3讲:流程及软件 。不过中间我还插播了一个文献解读 。
现在才正式开始数据处理实战,其中实战的测试数据,参考基因组以及对应的软件安装都是在第3讲:流程及软件 。看懂了这些准备工作,现在就可以跟我一起来一步步走通Hic数据分析流程啦,首先回忆一下准备工作咯。 conda create -n hic python=2 bowtie2 source activate hic conda search hiclab conda install -y sra-tools samtools conda install -y bedtools trim-galore source activate hic conda install -y numpy scipy matplotlib h5py cython numexpr statsmodels scikit-learn pandas cd ~ ln -s /data/training/zengjianming/ data mkdir -p /home/zengjianming/data/biosoft/hiclib cd /home/zengjianming/data/biosoft/hiclib pip install https:///mirnylab/mirnylib/get/tip.tar.gz pip install https:///mirnylab/hiclib/get/tip.tar.gz ## 17.7MB , 44kB/s mkdir -p /home/zengjianming/data/biosoft/hicpro cd /home/zengjianming/data/biosoft/hicpro git clone https://github.com/nservant/HiC-Pro.git ## 46.95 MiB | 134.00 KiB/s, done cd HiC-Pro/ mkdir /home/zengjianming/data/biosoft/hicpro/bin mkdir -p /data/training/zengjianming/biosoft/hicpro/bin/HiC-Pro_2.10.0 mkdir -p /home/zengjianming/bin
再次了解 Hic-pro 软件Hic-pro安装比较麻烦,我就不重复写了,大家仔细查看第3讲:流程及软件 。 很明显,前面教程我们讲到我们的Hic-pro安装在: source activate hic /home/zengjianming/biosoft/hicpro/bin/HiC-Pro_2.10.0/bin/HiC-Pro -h
但是使用的学员服务器没有配置好,导致我的home目录下面空间不足,所以我挪动了data盘下面的空间,这样就有28T的空间啦。 所以软件挪了地方: /data/training/zengjianming/biosoft/hicpro/bin/HiC-Pro_2.10.0/HiC-Pro/bin/HiC-Pro --help
然后我们这一讲的重点是Hic-pro软件的使用。 如果仅仅是看软件的帮助文档,会发现其实非常容易使用; usage : HiC-Pro -i INPUT -o OUTPUT -c CONFIG [-s ANALYSIS_STEP] [-p] [-h] [-v] -i|--input INPUT : input data folder; Must contains a folder per sample with input files -o|--output OUTPUT : output folder -c|--conf CONFIG : configuration file for Hi-C processing
重点就是3个参数而已,而最重要的当然是配置文件的准备咯。作者给出的示例文件不是一般的复杂; ## 记住 /data/training/zengjianming/ 目录就是 /home/zengjianming/data 目录。 cat /data/training/zengjianming/biosoft/hicpro/bin/HiC-Pro_2.10.0/HiC-Pro/config-hicpro.txt mkdir -p ~/data/project/hic/ cd ~/data/project/hic/
参考基因组mkdir -p ~/data/project/hic/ref cd ~/data/project/hic/ref wget ftp://ftp.ensemblgenomes.org/pub/bacteria/release-40/fasta/bacteria_20_collection/caulobacter_crescentus_na1000/dna/Caulobacter_crescentus_na1000.ASM2200v1.dna.toplevel.fa.gz gunzip Caulobacter_crescentus_na1000.ASM2200v1.dna.toplevel.fa.gz source activate hic bowtie2-build Caulobacter_crescentus_na1000.ASM2200v1.dna.toplevel.fa bacteria cat >genome.size Chromosome 4016942
测试数据来自于Tung B. K. Le et al. Science 2013 :https://www.ncbi.nlm./sra/?term=srr824846 mkdir -p ~/data/project/hic/fq/s1/ cd ~/data/project/hic/fq/s1/ 858M Jul 3 16:21 SRR824846_Q20L10_1.fastq.gz 857M Jul 3 16:22 SRR824846_Q20L10_2.fastq.gz
运行Hic-pro软件上面说过使用这个软件需要修改一个配置文件, 简单的包括: BOWTIE2_IDX_PATH = # bowtie2建立的索引所在的路径,记住绝对路径 REFERENCE_GENOME = # bowtie2建立的索引 GENOME_SIZE = # 一个文件记录着参考基因组中每条序列的大小,格式为chr01 10000000
复杂的是GENOME_FRAGMENT = ,也就是 HiC消化片段位点文件,这个文件在每个HIC流程中都需要生成,一般HiC建库的酶为hindiii(A^AGCTT ) 或者dnpiii 生成命令为 /PATH/HiC-Pro-master/bin/utils/digest_genome.py -r hindiii -o Refgenome.fasta 在软件例子里面是: GENOME_FRAGMENT = HindIII_resfrag_hg19.bed LIGATION_SITE = AAGCTAGCTT
但是我们这里的测试数据是 新月柄杆菌 Caulobacter crescentus , 文章里面 Le TB et al., 'High-resolution mapping of the spatial organization of a bacterial chromosome.', Science, 2013 Oct 24;342(6159):731-4 提到的是限制性酶切位点NcoI 也是 用的内切酶,其它还有 BamHI、EcoRI、HindIII、NdeI、XhoI, 参考:https://www.biomart.cn/experiment/430/457/741/43956.htm 那么我需要针对特定酶和特定的参考基因组序列来设置好GENOME_FRAGMENT参数所需文件。 代码如下: mkdir -p ~/data/project/hic/digest cd ~/data/project/hic/digest bin=/data/training/zengjianming/biosoft/hicpro/bin/HiC-Pro_2.10.0/HiC-Pro/bin/utils/digest_genome.py $bin -r C^CATGG -o bacteria.bed ../ref/Caulobacter_crescentus_na1000.ASM2200v1.dna.toplevel.fa
得到的 bacteria.bed 文件如下; Chromosome 0 3188 HIC_Chromosome_1 0 + Chromosome 3188 3551 HIC_Chromosome_2 0 + Chromosome 3551 3737 HIC_Chromosome_3 0 + Chromosome 3737 4500 HIC_Chromosome_4 0 + Chromosome 4500 5578 HIC_Chromosome_5 0 + Chromosome 5578 10834 HIC_Chromosome_6 0 + Chromosome 10834 10840 HIC_Chromosome_7 0 + Chromosome 10840 19183 HIC_Chromosome_8 0 + Chromosome 19183 19381 HIC_Chromosome_9 0 + Chromosome 19381 21967 HIC_Chromosome_10 0 +
然后修改配置文件后即可成功运行hicpro软件啦。 cd ~/data/project/hic/ source activate hic cp /data/training/zengjianming/biosoft/hicpro/bin/HiC-Pro_2.10.0/HiC-Pro/config-hicpro.txt ./ bin=/data/training/zengjianming/biosoft/hicpro/bin/HiC-Pro_2.10.0/HiC-Pro/bin/HiC-Pro nohup $bin -i fq -o out -c config-hicpro.txt 1> run.log 2>&1 &
这个配置文件问题简单,重要的地方如下: BOWTIE2_IDX_PATH = /home/zengjianming/data/project/hic/ref/ REFERENCE_GENOME = bacteria GENOME_SIZE = /home/zengjianming/data/project/hic/ref/genome.size PAIR1_EXT = SRR824846_Q20L10_1 PAIR2_EXT = SRR824846_Q20L10_2 GENOME_FRAGMENT = /home/zengjianming/data/project/hic/digest/bacteria.bed LIGATION_SITE = CCATGG
但是 我遇到的报错居然是在 input的fastq文件的folder这个参数,太尴尬:https://groups.google.com/forum/#!msg/hic-pro/PrX00iA0N1U/apzlw8ArBwAJ 也就是说每个样本的fq文件都需要在自己独立的文件夹下面。 因为软件包含的步骤较多,所以会比较耗时! 输出结果解读所有的输出,都在自己运行软件的时候指定的out文件夹下面。 共有bowtie_results和hic_results两个结果文件夹。 其中bowtie_results文件夹下共有三个文件夹:bwt2、bwt2_global和bwt2_local,分别是序列比对结果、染色体间关联比对结果和染色体内部关联比对结果。其中bwt2文件夹下有一些数据统计结果的输出文件,如mpairstat文件(如下)。 Total_pairs_processed 23546961 100.0 Unmapped_pairs 22197 0.094 Low_qual_pairs 0 0.0 Unique_paired_alignments 17744193 75.357 Multiple_pairs_alignments 462162 1.963 Pairs_with_singleton 5318409 22.586 Low_qual_singleton 0 0.0 Unique_singleton_alignments 0 0.0 Multiple_singleton_alignments 0 0.0 Reported_pairs 17744193 75.357
hic_results文件夹下共有data、matrix以及pic三个文件夹。 重点必须是hic_results文件,如下: hic_results/ |-- [ 24] data | `-- [4.0K] s1 | |-- [504M] _bacteria.bwt2pairs.DEPairs | |-- [ 93K] _bacteria.bwt2pairs.DumpPairs | |-- [ 83M] _bacteria.bwt2pairs.REPairs | |-- [ 303] _bacteria.bwt2pairs.RSstat | |-- [ 73M] _bacteria.bwt2pairs.SCPairs | |-- [ 0] _bacteria.bwt2pairs.SinglePairs | |-- [1.2G] _bacteria.bwt2pairs.validPairs | |-- [1.2G] s1_allValidPairs | |-- [ 150] s1_allValidPairs.mergestat | `-- [ 486] s1.mRSstat |-- [ 24] matrix | `-- [ 41] s1 | |-- [ 99] iced | | |-- [ 85] 1000000 | | | |-- [ 275] s1_1000000_iced.matrix | | | `-- [ 125] s1_1000000_iced.matrix.biases | | |-- [ 83] 150000 | | | |-- [6.6K] s1_150000_iced.matrix | | | `-- [ 675] s1_150000_iced.matrix.biases | | |-- [ 81] 20000 | | | |-- [345K] s1_20000_iced.matrix | | | `-- [4.9K] s1_20000_iced.matrix.biases | | |-- [ 81] 40000 | | | |-- [ 86K] s1_40000_iced.matrix | | | `-- [2.4K] s1_40000_iced.matrix.biases | | `-- [ 83] 500000 | | |-- [ 792] s1_500000_iced.matrix | | `-- [ 225] s1_500000_iced.matrix.biases | `-- [ 99] raw | |-- [ 99] 1000000 | | |-- [ 139] s1_1000000_abs.bed | | |-- [ 165] s1_1000000.matrix | | `-- [ 18] s1_1000000_ord.bed -> s1_1000000_abs.bed | |-- [ 96] 150000 | | |-- [ 783] s1_150000_abs.bed | | |-- [4.1K] s1_150000.matrix | | `-- [ 17] s1_150000_ord.bed -> s1_150000_abs.bed | |-- [ 93] 20000 | | |-- [5.9K] s1_20000_abs.bed | | |-- [216K] s1_20000.matrix | | `-- [ 16] s1_20000_ord.bed -> s1_20000_abs.bed | |-- [ 93] 40000 | | |-- [2.9K] s1_40000_abs.bed | | |-- [ 53K] s1_40000.matrix | | `-- [ 16] s1_40000_ord.bed -> s1_40000_abs.bed | `-- [ 96] 500000 | |-- [ 253] s1_500000_abs.bed | |-- [ 469] s1_500000.matrix | `-- [ 17] s1_500000_ord.bed -> s1_500000_abs.bed `-- [ 24] pic `-- [ 188] s1 |-- [4.9K] plotHiCContactRanges_s1.pdf |-- [5.2K] plotHiCFragment_s1.pdf |-- [6.3K] plotHiCFragmentSize_s1.pdf |-- [5.0K] plotMappingPairing_s1.pdf `-- [5.1K] plotMapping_s1.pdf
18 directories, 40 files
因为这个菌参考基因组很小,才4M左右,所以取500K,就只有8个bin了,其实意义不大。 随便秀一个pdf图吧。
似乎是还停留在对比对结果的统计层面上!
后续再更新结果的生物学意义解读,以及后续的compartment以及TAD和loop的分析,敬请期待哦! 里面提到的各种文献,综述,测试数据,软件,待这个三维基因组更新完毕会统一打包发放哦,请大家不要着急。
|