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。之后如果用的到的话会再仔细的深入研究一下使用方法。