1.输入文件的准备
创建一个文件夹，以物种名命名，比如:
mkdir spalax_galili
在spalax_galili文件夹中再创建2个文件夹，分别命名为fastq和reference，如下图
在fastq文件夹中软连接HiC测序数据，hic测序数据命名为“ 物种名.HiC_R1.fastq.gz / 物种名.HiC_R2.fastq.gz “。可以压缩，也可以不压缩，但是要严格按照这个格式命名，如下图：
#对参考基因组建立索引
bwa index spalax_galili.fasta #参考基因组的bwa索引
samtools faidx spalax_galili.fasta
cut -f1,2 spalax_galili.fasta.fai > spalax_galili.length #参考基因组的染色体长度文件
#参考基因组的酶切文件
python /data/01/user164/software/juicer/misc/generate_site_positions.py MboI spalax_galili spalax_galili.fasta
##会生成一个文件：spalax_galili_MboI.txt，这个文件后续会用到
##106服务器的路径是:/home/106public/software/juicer/misc/generate_site_positions.py
##python generate_site_positions.py 测序用的内切酶名 物种名 参考基因组
在reference文件夹中存放参考基因组的bwa索引和染色体长度文件，如图所示：
2.运行juicer
最后确定所有的输入文件已经准备好。Hi-C双端测序数据、参考基因组以及bwa索引、参考基因组的染色体长度文件、参考基因组的酶切文件，如下图：
在以物种名命名的文件夹，和fastq/reference同级的目录写shell脚本
vim spalax_galili.juicer.sh
/data/01/user164/software/juicer/scripts/juicer.sh -z 参考基因组 -p 染色体长度文件 -y 参考基因组酶切文件 -d 以物种名命名的文件夹 -D /data/01/user164/software/juicer -g 物种名 -s 内切酶名 -t 线程数
#-D参数是软件路径, 不变即可
#在超算上测试耗费资源, 内存在70-80G左右。
#CPU 10-40个(10-40个运行的时间差不多相同，都是2-3天)。
#106服务器的路径:
#/home/106public/software/juicer/scripts/juicer.sh
#-D 设置为: /home/106public/software/juicer
#比如（请注意绝对路径的问题）：
/home/kuangzhr17/xlliang/software/juicer/scripts/juicer.sh -z ./spalax_galili.fasta -p ./spalax_galili.length -y ./spalax_galili_MboI.txt -d /home/kuangzhr17/ Pan.3DGenom/spalax_galili -D /home/106public/software/juicer -g spalax_galili -s MboI -t 40
运行完毕后，会生成aligned的文件夹
里面的inter.hic/inter_30.hic，即最后的输出文件, 其余文件和文件夹都可以删除
改名，改为 物种名.hic/物种名.30.hic，比如spalax_galili.hic/spalax_galili.30.hic
后续区室的鉴定/TAD的鉴定等都要基于此文件
如果有使用Centromics鉴定着丝粒序列的需求, 可以保留merged_nodups.txt
https://github.com/zhangrengang/Centromics
- AB区室、TAD的鉴定
使用软件https://h1d.readthedocs.io/en/latest/
spalax_galili的染色体号是这么命名的：Chr1/Chr2/Chr3
以Chr1为例, 其余染色体批量生成命令行即可
3.1 互作矩阵的提取 (VC标准化)
h1d basic dump spalax_galili.30.hic 25000 Chr1 --datatype rawhic -o spalax_galili.VC –gt spalax_galili.length –normalize VC
#h1d basic dump juicer输出文件 分辨率 染色体号 --datatype rawhic(照抄) -o 物种名 -gt 染色体长度 –normalize 标准化方式-VC/KR
#会输出一个叫spalax_galili.VC的文件夹，里面有个叫25000(分辨率)的文件夹，再进去就是Chr1互作矩阵。不同的染色体不同的分辨率都指定同一个文件夹即 -o 参数都设置相同的路径。
#最后的互作矩阵是：observed.VC.chrChr1.matrix.gz
#这里是25k分辨率的互作矩阵，还需要提取5k，10k，50k分辨率的互作矩阵
3.2 AB区室的鉴定（以50k分辨率为例）
#基因密度的计算 – 需要基因注释文件
bedtools makewindows -g spalax_galili.length -w 50000 > 50k.windows
grep -w "gene" input.gff | awk '{print $1"\t"$4"\t"$5}' > gene.pos
bedtools intersect -a 50k.windows -b gene.pos -c > spalax_galili.50k.gene.density
h1d one PC1 observed.VC.chrChr1.matrix.gz 50000 Chr1 -o spalax_galili.Chr1.PC1 --gt spalax_galili.length -p spalax_galili.50k.gene.density
# spalax_galili.Chr1.PC1即最后的AB区室文件，以基因密度更高的作为A区室
3.3 TAD的鉴定（以25k分辨率为例）
h1d call stripeTAD observed.VC.chrChr1.matrix.gz 25000 Chr1 -o spalax_galili.Chr1.TAD --gt spalax_galili.length
4.这个软件还能算什么，还能鉴定什么
