简介
在本文中，我们将使用 Caki2 细胞系的 Hi-C 数据来说明利用 Hi-C 数据发现 SVs 的过程。
数据集
Caki2 是一种人类肾癌细胞系。我们将从 ENCODE（https://www.encodeproject.org/，登录号：ENCSR401TBQ）下载公开可用的 Caki2 Hi-C 数据集。
数据下载
首先，我们将建立一个工作目录，然后把数据下载到该目录中，例如，在你的 home 目录下创建一个名为 “project” 的文件夹。
mkdir project && cd project
wget -O Caki2.R1.fastq.gz https://www.encodeproject.org/
files/ENCFF190NLS/@@download/ENCFF190NLS.fastq.gz
wget -O Caki2.R2.fastq.gz https://www.encodeproject.org/
files/ENCFF914EQR/@@download/ENCFF914EQR.fastq.gz
计算方法
已经提出了若干种计算方法，利用 Hi-C 作为工具来检测结构变异和拷贝数变异（CNVs）。
在这些方法中，Hi-C breakfinder 能够以高达 1 kb 的分辨率检测出最多种类型的 SVs。在本文中，我们将专注于使用 Hi-C breakfinder 发现 SVs。
Hi-C 分析流程
Hi-C 数据分析通常包括以下步骤：
- 将 Hi-C reads 映射到参考基因组；
- 生成并归一化 Hi-C 矩阵；
- 检测染色质结构的层级，包括 A/B compartments、Topologically Associating Domains（TADs）、sub-TADs 或 loop domains。
已有多项计算方法被开发出来，用以简化上述流程。
Hi-C breakfinder
Hi-C breakfinder 唯一需要用户提供的输入是 reads 的比对文件。因此，人们可以使用由各种 Hi-C 分析流程生成的比对文件（通常是 “bam” 文件）。
根据我们的观察，不同流程产生的比对文件之间的差异对 Hi-C breakfinder 的结果没有影响。在接下来的内容中，我们将展示一条内部流程，从 Hi-C reads 的比对开始，到利用 Hi-C breakfinder 发现 SVs，并最终在 Hi-C 图上可视化 SVs。
Pipeline
我们首先使用一条改编自 4D Nucleome（https://data.4dnucleome.org/resources/data-analysis/hi_c-processing-pipeline）的流程来处理原始 Hi-C reads。
此步骤将生成一个比对文件（.bam），它将成为下一步 SV 检测的输入；同时还会生成一个 cooler 格式的矩阵文件，我们可借此在 SV 位点绘制 Hi-C 热图。
实战
我们使用 BWA aligner 将 Caki2 原始 reads 比对到人类基因组参考序列 GRCh38。
用户可以根据自身需求选择参考版本。首先，我们需要下载参考数据并为参考序列建立索引。
mkdir -p project/ref; cd project/ref
wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/405/
GCA_000001405.15_GRCh38/seqs_for_alignment_pipelines.ucs-
c_ids/GCA_000001405.15_GRCh38_no_alt_plus_hs38d1_analysis_-
set.fna.gz
gunzip GCA_000001405.15_GRCh38_no_alt_plus_hs38d1_analysis_-
set.fna.gz
bwa index -a bwtsw GCA_000001405.15_GRCh38_no_alt_plus_hs38-
d1_analysis_set.fna
接下来，我们将 Caki2 的paired-end reads 比对到已建立索引的参考基因组。
请注意，旧版流程可能要求用户分别比对每一对 reads，因为在比对过程中难以准确估计插入片段大小。在新版（> 0.7.15）的 BWA 中已无此必要，可以一次性用一条命令把一对 reads 的两个 mate 同时比对，而无需强制它们来自同一基因组片段。
cd project
REF= ref/GCA_000001405.15_GRCh38_no_alt_plus_hs38d1_analy-
sis_set.fna
bwa mem -SP5M -t 8 $REF Caki2.R1.fastq.gz Caki2.R2.fastq.gz |
samtools view -bS > Caki2.bam
接下来，我们要用 Pairtools 过滤掉未比对、单侧比对、重复比对以及多重比对的 reads。最小比对质量由 --min-mapq 设定。在 Juicer 中常用的两个值是 1 和 30，这里我们把 --min-mapq 设为 1。
Pairtools 会把比对结果转换成接触对（contact pairs）。若一条 read 满足 --min-mapq 定义的“唯一比对”，则其配对类型标记为 “U”；若来自嵌合比对但仍是有效的 Hi-C 接触，则标记为 “R”。我们将保留标记为 “UU”、“UR” 或 “RU” 的配对。
CHROM_SIZES=/path/to/GRCh38.chrom.sizes
samtools view Caki2.bam | \
pairtools parse -c $CHROM_SIZES --min-mapq 1 --assembly hg38 |
pairtools sort | \
pairtools select ’(pair_type == "UU" or pair_type == "UR" or
pair_type == "RU")’ | \
pairtools dedup --max-mismatch 3 --output \
>( pairtools split --output-pairs Caki2.nodups.pairs.gz --
output-sam Caki2.nodups.bam )
由 Pairtools 生成的 “Caki2.nodups.bam” 将作为 Hi-C breakfinder 的输入。
如果已知限制酶信息，可以进一步过滤掉落在同一限制片段内的 reads 或自环（self-circles）。为此，我们需要先生成一个限制位点文件，然后用该文件对配对 reads 进行过滤。这里我们以 “MboI” 限制酶为例。
cooler digest -o restrict_sites.MboI.hg38.bed $CHROM_SIZE $REF
MboI
pairtools select ’len(chrom1) < 6 and len(chrom2) < 6 and
chrom1 != "chrM" and chrom2 != "chrM"’ Caki2.nodups.pairs.gz |
\ paritools restrict -f restrict_sites.MboI.hg38.bed | \
awk ’$0 ~ /^#/{print} $0 !~ /^#/ && $9 != $12{print}’ | bgzip
-c > Caki2.nodups.valid.pairs.gz
未完待续，欢迎关注！
