Hi-C分析练习（从fastq文件到contact矩阵）
最近纽约复工，因为每天要做实验（还得全程戴口罩，非常闷得慌），比不了之前每天都在家可以为所欲为的自学生信，所以更新也少了很多~
Hi-C的技术相关背景学习笔记我在之前有写过（Hi-C技术的初步了解）。本篇笔记按照哈佛医学院贴在网上的教程走一遍流程，虽然这个教程是2018年的，但我觉得仍然很有参考价值。具体可参考网址：https://github.com/hms-dbmi/hic-data-analysis-bootcamp
（一）下载练习文件
根据教程，你需要先下载练习用的fastq文件和相关基因组索引文件：
# download input fastq files下载fastq练习文件
$ wget https://s3.amazonaws.com/4dn-dcic-public/hic-data-analysis-bootcamp/input_R1.fastq.gz
$ wget https://s3.amazonaws.com/4dn-dcic-public/hic-data-analysis-bootcamp/input_R2.fastq.gz
$ ll
-rw------- 1 fy04 fy04 27701813 May 16 2018 input_R1.fastq.gz
-rw------- 1 fy04 fy04 27626086 May 16 2018 input_R2.fastq.gz
#解压fastq文件
$ gunzip input_R1.fastq.gz
$ gunzip input_R2.fastq.gz
$ ll
-rw------- 1 fy04 fy04 98637520 May 16 2018 input_R1.fastq
-rw------- 1 fy04 fy04 91381760 May 16 2018 input_R2.fastq
# download bwa genome index下载bwa的基因组索引文件，你也可以自己构建
$ wget https://s3.amazonaws.com/4dn-dcic-public/hic-data-analysis-bootcamp/hg38.bwaIndex.tgz
$ tar -xzf hg38.bwaIndex.tgz
$ rm hg38.bwaIndex.tgz
$ ll
-rw------- 1 fy04 fy04 17478 May 11 2018 hg38.fasta.amb
-rw------- 1 fy04 fy04 27715 Dec 8 2016 hg38.fasta.ann
-rw------- 1 fy04 fy04 3099922624 Dec 8 2016 hg38.fasta.bwt
-rw------- 1 fy04 fy04 774980637 Dec 8 2016 hg38.fasta.pac
-rw------- 1 fy04 fy04 1549961320 Dec 8 2016 hg38.fasta.sa
# download chromsizes
$ wget https://s3.amazonaws.com/4dn-dcic-public/hic-data-analysis-bootcamp/hg38.mainonly.chrom.size
这时你应该下好的文件有：
然后可以简单的查看一下我们下载的文件：
#查看fastq文件
$ head ./input_R1.fastq
@1
GATCACACCACTGCACTCCAGCCTGGGCGACAGAGCAAGACTCTATCTCAAAAAAAAAAAAAAAAATCAAGAAATAATGGCTGAAAATTTCCAAAATTTG
+1
AABBBFFFB@BFGGGFFGGGFGGHHHCGCGGGGHFBGCFFGFFFHHHHHGHGBEE0EEE/EEEEEE03?F33BF4BEB22?2/B2BFD22FFDFGHFCG2
@2
CATCAGAAGGATGGATCCCTTGTGCTACCCATTTTATAGTCACAAACACGTTCCTCACTCCCTGACCCCTGGCAACCACTAATGGGTTTTTCATCTCTGT
+2
CCCCCFFFFFFFGGGGGGGGGGHHHHHHHHHHHHHHHHHHHGHHHHHHHGHHHHHHHGHHHHHHHHHHGGGHHHHHHGGHHEGHHGFFGGGHHHHHHHHE
@3
CTTTTTTTTGTAAAATGAGAAAGCTGGATGAAAGGATCTCTAAGGTCCCTTCACCTTTGAAATCTGATAACTCTAAGGTTTAATTTCCTCATCTGTCAAA
#查看基因组大小文件
$ head ./hg38.mainonly.chrom.size
chr1 248956422
chr2 242193529
chr3 198295559
chr4 190214555
chr5 181538259
chr6 170805979
chr7 159345973
chr8 145138636
chr9 138394717
chr10 133797422
（二）bwa比对
上面的文件都下载好后，就可以按照官方给出的流程PPT走分析流程了：
https://hms-dbmi.github.io/hic-data-analysis-bootcamp/#1
第一步当然就是比对了。在Hi-C的分析流程里，比对使用的软件是bwa。这不同于以往使用的bowtie2、hisat2和STAR（比对软件STAR的使用）。bwa软件的官方网站请见：here。
#这里用bwa比对后，直接使用管道符号把sam文件转成bam文件
$ bwa mem -SP5M -t8 hg38.fasta input_R1.fastq input_R2.fastq | samtools view -bhS - > output.bam
#-t8: use 8 cores
#-SP5M : Hi-C-specific options
#查看bam文件
$ samtools view output.bam | less -S
1 97 chr7 99922629 48 49M1I50M = 76330531 -23592008 GATCACACCACTGC
1 145 chr7 76330531 60 92M = 99922629 23592008 AGTACAGATGGTCTAACTCAGT
2 113 chrY 56859590 26 100M chrUn_GL000216v2 129167 0 ACAGAGATGAAAAACCCATTAG
2 177 chrUn_GL000216v2 129167 0 92M chrY 56859590 0 TTATCCATATATTCACCTGTCT
3 81 chr10 110772888 60 100M = 110908302 135316 TTTGACAGATGAGGAAATTAAACCTTAGAG
3 161 chr10 110908302 60 92M = 110772888 -135316 NGACATCTTTTTCTCCTATCACAAGGATAT
4 97 chr18 65855270 60 83M1I16M = 63677027 -2178190 TTTTACAGAAATCA
（三）过滤
把上面得到的bam文件进行过滤。过滤这一步官网使用的软件是pairsamtools（这个软件现已更名为pairtools）。
pairtools软件官网：https://readthedocs.org/projects/pairtools/downloads/pdf/latest/
这里主要是要生成.pairs文件。这里涉及到Hi-C的技术原理，请参看我之前的文章（Hi-C技术的初步了解），先简单的了解Hi-C的样品取样流程，这里更容易理解，文字可能不好表达，请意会：
而pairtools在这一步是干嘛的呢？见下图：
安装pairtools：
$ conda install -c conda-forge -c bioconda pairtools
或者
$ pip install pairtools
接下来就是过滤步骤了：
#读取bam文件，并生成Hi-C pairs
$ samtools view -h output.bam | pairtools parse -c hg38.mainonly.chrom.size -o parsed.pairsam
#对.pairsam文件进行sort
$ pairtools sort --nproc 8 -o sorted.pairsam parsed.pairsam
#查看输出文件
$ less -S parsed.pairsam
这时你查看的文件应该是长这样：
（四）去重
这一步仍然需要用pairtools这个软件。
$ pairtools dedup --mark-dups -o deduped.pairsam sorted.pairsam
（五）Select
根据pair类型筛选pairs。
$ pairtools select \
'(pair_type == "UU") or (pair_type == "UR") or (pair_type == "RU")' \
-o filtered.pairsam deduped.pairsam
（六）拆分文件
把上面得到的.pairsam文件拆分成.pairs文件和.sam文件。
$ pairtools split --output-pairs output.pairs filtered.pairsam
#查看output.pairs文件
$ less -S output.pairs
（七）对pairs文件建立索引以及查找区域内的pairs
这一步需要用到的软件是：pairix
pairix软件官方网址：https://github.com/4dn-dcic/pairix
首先安装pairix:
$ conda install -c bioconda pairix
#把上一步得到的output.pairs文件进行压缩，得到gz文件
$ bgzip output.pairs
#建立索引
$ pairix -f output.pairs.gz
#可以试着查询一下染色体某一区域内的pairs
$ pairix output.pairs.gz 'chr2:1-60000000|chr4:5000000-6000000'
156768 chr2 16158981 chr4 5970814 - + UU
448703 chr2 31363676 chr4 5157768 + + UU
199145 chr2 44090655 chr4 5811114 - - UU
#上面的各列代表: readID chr1 pos1 chr2 pos2 strand1 strand2 pair_type
（八）Binning! (生成contact矩阵)
这一步需要用的软件是：cooler。软件官网：https://github.com/mirnylab/cooler。
首先安装cooler：
$ conda install -c conda-forge -c bioconda cooler
输入文件需要: pairs文件, chromsize文件
输出文件 : cool文件
bin大小 : 比如5000 (高分辨率), 500000 (低分辨率)
$ cooler cload pairix hg38.mainonly.chrom.size:50000 output.pairs.gz output.cool
INFO:cooler.cli.cload:Using 8 cores
INFO:cooler.create:Creating cooler at "output.cool::/"
INFO:cooler.create:Writing chroms
INFO:cooler.create:Writing bins
......
INFO:cooler.create:Writing indexes
INFO:cooler.create:Writing info
INFO:cooler.create:Done
利用cooler读取cool文件：
$ cooler dump -t pixels --header --join -r chr19 output.cool
chrom1 start1 end1 chrom2 start2 end2 count
chr19 250000 300000 chr19 250000 300000 4
chr19 250000 300000 chr19 300000 350000 1
chr19 250000 300000 chr19 350000 400000 1
chr19 250000 300000 chr19 5300000 5350000 1
chr19 250000 300000 chr19 6900000 6950000 1
......＃这里你会发现最后一列是count，并且都是整数
#上面各列代表：chrom1 start1 end1 chrom2 start2 end2 value
（九）标准化
我们在做RNA-seq和CHIP-seq以及其他各种高通量测序实验分析中，都有一步是标准化。Hi-C也一样。需要把上面我们得到的这个output.cool文件进行标准化。
#标准化
$ cooler balance output.cool
# 查看输出文件(w/ -b for normalized matrix)
$ cooler dump -t pixels --header --join -r chr19 -b output.cool
chrom1 start1 end1 chrom2 start2 end2 count balanced
chr19 250000 300000 chr19 250000 300000 4
chr19 250000 300000 chr19 300000 350000 1
chr19 250000 300000 chr19 350000 400000 1
......
chr19 900000 950000 chr19 1300000 1350000 1 0.323781
chr19 900000 950000 chr19 2600000 2650000 1 0.249956
......
#标准化后会多出一列
（十）Aggregation! (For HiGlass view)
这一步是为了进行后续使用HiGlass进行可视化的。这一步会生成一个多分辨率的cool文件。
$ cooler zoomify output.cool
后续的可视化可以按照教程https://docs.higlass.io/higlass_docker.html操作了。哈佛的这个Hi-C流程的可视化使用的是Higlass，也有其他的可视化软件。比如Juicebox和WashU Epigenome Browser。之后如果用的到的话会再仔细的深入研究一下使用方法。
