单倍型基因组组装工具 --- ALLHic

发表于 更新于 分类于 ， 阅读次数： 本文字数： 11k 阅读时长 ≈ 10 分钟

背景

单倍型基因组组装方法是针对于动植物基因组组装提出的更为精确的组装策略，二倍体基因组同时拥有来自父本与母本的基因型，存在复杂的等位基因变异，极大地提高了二倍型基因组的组装难度。也就是说我们如要完整组装二倍型基因组，需要将二倍型基因组进行两次单倍型基因组组装，从而精准的定位到染色体水平。目前，大多数二倍体基因组组装忽略了同源染色体之间的差异，从而将基因组组装成一个伪单倍体序列，这将会对之后基因注释产生很多不确定因素以及得到不准确的生物学解释。想要得到二倍体精确的单倍体基因组，从测序样本入手是一种常见的策略。我们可以通过不断近交以得到纯合双单倍体基因型，便可组装得到纯合基因组；又或者提供物种的配子进行测序，亦能得到准确的单倍型基因组；又或者在测序前分离单个染色体，从而实现单个单倍型基因组的测序和组装。

虽然通过各种技术想方设法从测序样本中得到单个染色体进行测序组装，能够满足部分需求；但是实验技术也存在诸多不确定的局限。如染色体分选局限于能否通过荧光强度或光散射区分特定染色体；单细胞技术 Strand-seq 局限于其数据获取困难。所以，将重心转移到二倍体基因组组装策略是一个不错的思路

认识染色体

在介绍组装工具之前，有必要先了解一下生物体内的染色体构成。就拿人染色体来讲，可以表示为 2n=2x=46 ；其中 n 代表的是单倍体染色体数目，也就是我们今天提到的组装工具 ALLHiC 的组装目标，n 是正常细胞的同源染色体的对数，则 2n 就代表着生物体体细胞的染色体数目。x 是有别于 n 的概念，代表的是染色体基数，表示一套染色体组的染色体数量，其也决定了生物体基因组的倍性。如人类基因组就是二倍体，则表示为 2x=46 ；代表人类一套完整的染色体组包含 23条 染色体；2n=46 则代表着人类配子（精子/卵子）携带的染色体数量为 23条。

染色体组：指形态、大小各不相同的一组染色体

n：生物体配子体的染色体数量 = 体细胞同源染色体对数 = 单倍体染色体数量

x：生物体一个染色体组中染色体数量 = 含有控制生物所有性状的一整套基因的染色体数量（一套染色体组中没有重复的基因，同时也没有缺失，即控制同一性状的基因只有一个，同时在一个染色体组中不存在等位基因）

这里再举个四倍体的例子：2n=4x=32 ；那我们能得到的信息如下：

1
2
3
1. 该物种的染色体总数：`32`
2. 单倍体染色体条数：`16`  --- 1n=2x=16；即来自亲本的配子有 `16` 条。
3. 一个单倍体内有两套染色体组，四倍体内有四套染色体组，各来自两个亲本的两套。

今天主要是看二倍体基因组单倍体组装流程，采用的组装工具是 ALLHiC ；二倍体基因组是由两个同源的染色体拷贝组装的，两个亲本各有一个，而我们要得到的结果便是其中一个同源染色体拷贝，即具有一套完整染色体信息的单倍体基因组；通过 ALLHic 组装，最终得到两组染色体水平的单倍型基因组。

ALLHiC 的诞生

2019年8月5日，福建农林大学基因组中心张兴坦副教授和唐海宝教授研究组在Nature Plants杂志在线发表题为“Assembly of allele-aware, chromosomal-scale autopolyploid genomes based on Hi-C data”的研究论文，开发了一种 ALLHiC 算法，成功完成了同源四倍体和同源八倍体甘蔗染色体水平组装。ALLHiC 可以从头组装获得单倍型染色体水平的多倍体基因组，分离每一个等位基因。ALLHiC 同时适用于多倍体和杂合二倍体的单倍型染色体水平组装。

作者对 ALLHiC 所有代码进行开源

Github：tangerzhang/ALLHiC: ALLHiC: phasing and scaffolding polyploid genomes based on Hi-C data (github.com)

Pipeline：Home · tangerzhang/ALLHiC Wiki (github.com)

ALLHiC 组装策略

Allelic contig table --- 等位基因重叠群表（通过参考基因组生成）

BAM --- HiC 数据 mapping 的 BAM

Prunning --- 修剪；连续的红线表示折叠区域；紫色、蓝色和绿色的连续线表示假设的多倍体基因组中的不同单倍型。粉红色虚线表示等位基因间 Hi-C 信号，而灰色虚线表示折叠区域和非折叠重叠群之间的 Hi-C 信号。

Partitioning --- 分区；彩色圆点表示在真实数据集中具有原始染色体分配的重叠群。

Rescue --- 救援（一个检查过程）；不同的颜色表示来自四种不同单倍型的重叠群，而红色实线表示折叠的重叠群。

Optimization --- 优化；不同的颜色表示四种不同的单倍型。

Building --- 构建；生成 group.asm.fasta，染色体水平的单倍体基因组

我们采用 canu 组装的伪单倍体 contigs 作为 ALLHiC 的输入。

Allhic 在很大程度上依赖于初始 contigs 组装的质量。单倍型之间序列高相似性会导致不明确的 reads 比对，可能会导致大量的 HiC reads 在 ALLHiC 管道中被删除，从而影响组装效果。

ALLHiC 组装二倍体

接下来我们就利用 ALLHiC 对杂合二倍体 Prunus armeniaca 进行组装，以获得染色体水平的单倍型基因组。

环境搭建 + 程序安装

首先，配置 ALLHiC 运行所需的依赖，我们利用 miniconda3 来配置相关环境。

  • samtools v1.9+
  • bedtools
  • matplotlib v2.0+
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
conda create -y -n allhic python=3.7 samtools bedtools matplotlib
## 下载二进制文件，解压之后可直接运行，无需编译
wget https://github.com/tanghaibao/allhic/releases/download/v0.9.13/allhic_0.9.13_Linux_64-bit.tar.gz  
## 解压
tar -zxvf allhic_0.9.13_Linux_64-bit.tar.gz
./allhic --version  ## 查考 allhic 版本号
## allhic version 0.9.13

## 文件结构
tree
.
|-- LICENSE
|-- README.md
|-- allhic
`-- allhic_0.9.13_Linux_64-bit.tar.gz

0 directories, 4 files

## 上面的方法是直接下载 ALLHiC 的二进制文件，省去了编译的过程，在唐海宝教授的 Github 上的；但是我之前是根据张兴坦教授在 Github 留的 Pipline 进行组装的，所以下面记录一下另外一种运行 ALLHiC 的方法。
git clone https://github.com/tangerzhang/ALLHiC
cd ALLHiC
chmod +x bin/*
chmod +x scripts/*  
export PATH=/your/path/to/ALLHiC/scripts/:/your/path/to/ALLHiC/bin/:$PATH
## 所需的依赖参照上面，利用 miniconda 进行配置

测序数据下载

ALLHiC 运行前的数据准备：

  • 具有良好 N50 指标的 contig 水平的目标基因组组装
  • Hi-C reads

我们通过 Canu 对 PacBio CLR 数据进行初步组装，从而得到 contig 水平的基因组，首先，我们从 ENA 数据库中下载需要的 PacBio 原始测序数据。这里利用 ascp 进行高速下载。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
ascp -v -QT -l 400m -P33001 -k1 -i \
                    ~/.aspera/connect/etc/asperaweb_id_dsa.openssh \
                    --mode recv --host fasp.sra.ebi.ac.uk --user era-fasp \
                    --file-list download.txt ./
								
## download.txt
cat download.txt
## Hi-C Data
vol1/run/ERR461/ERR4619488/4672_A_run658_ATCACGAT_S2_L005_R1_001.fastq.gz
vol1/run/ERR461/ERR4619488/4672_A_run658_ATCACGAT_S2_L005_R2_001.fastq.gz
vol1/run/ERR461/ERR4619489/4672_B_run658_GCCAATAT_S3_L005_R1_001.fastq.gz
vol1/run/ERR461/ERR4619489/4672_B_run658_GCCAATAT_S3_L005_R2_001.fastq.gz
## PacBio CLR Data
vol1/run/ERR409/ERR4093216/m54124_180918_044155.subreads.bam
vol1/run/ERR409/ERR4093217/m54124_180917_082648.subreads.bam
vol1/run/ERR409/ERR4093218/m54124_180911_161116.subreads.bam

## 因为数据库提供了多次测序数据，所以我们需要将多个合并成一个进行下面的分析
cat 4672_A_run658_ATCACGAT_S2_L005_R1_001.fastq.gz 4672_B_run658_GCCAATAT_S3_L005_R1_001.fastq.gz \
                    > hic_R1.fastq.gz
cat 4672_A_run658_ATCACGAT_S2_L005_R2_001.fastq.gz 4672_B_run658_GCCAATAT_S3_L005_R2_001.fastq.gz \
                    > hic_R2.fastq.gz
																			
## 合并 bam 文件
samtools merge pacbio.clr.bam m54124_180911_161116.subreads.bam \
                    m54124_180917_082648.subreads.bam m54124_180918_044155.subreads.bam

初步组装

下载好测序数据之后，我们就可以对 PacBio CLR 进行初步组装，如下：

我们利用 Canu 工具对 PacBio 数据进行组装，从而得到 contigs；

1
2
3
4
5
6
7
## 上面我们从 ENA 数据库中下载的是 bam 文件，所以我们需要先对其格式转换成 fasta 序列
name=pacbio.clr
samtools bam2fq ${name}.bam > ${name}.fastq
awk '{if(NR%4 == 1){print ">" substr($0, 2)}}{if(NR%4 == 2){print}}' ${name}.fastq > ${name}.fasta
canu -p preasm -d canu_result useGrid=false genomeSize=250000000 \
                    corMhapSensitivity=high corMinCoverage=0 corOutCoverage=100 \
                    correctedErrorRate=0.105 -pacbio-raw ${name}.fasta

Canu 组装很耗时间，耐心等待之后，我们便得到了初步组装的 contigs 序列。

ALLHiC 以获得染色体水平基因组组装

  1. Map Hi-C reads to draft assembly

draft.asm.fasta 是我们上面利用 Canu 组装得到的 contigs

reads_R*.fastq.gz 是 Hi-C 的双端数据

1
2
3
4
5
6
7
8
9
10
11
12
13
14
## 建立索引
bwa index -a bwtsw draft.asm.fasta  
samtools faidx draft.asm.fasta  

## 将 Hi-C reads 比对到 contigs
bwa aln -t 24 draft.asm.fasta reads_R1.fastq.gz > sample_R1.sai  
bwa aln -t 24 draft.asm.fasta reads_R2.fastq.gz > sample_R2.sai  
bwa sampe draft.asm.fasta sample_R1.sai sample_R2.sai reads_R1.fastq.gz reads_R2.fastq.gz \
                    > sample.bwa_aln.sam  
## 过滤 SAM 文件
PreprocessSAMs.pl sample.bwa_aln.sam draft.asm.fasta MBOI
(*)filterBAM_forHiC.pl sample.bwa_aln.REduced.paired_only.bam sample.clean.sam
samtools view -bt draft.asm.fasta.fai sample.clean.sam > sample.clean.bam
## *Tip: skip this step if you are using bwa mem for alignment（如果上面的比对采用的是 bwa mem ； 则跳过这个步骤）
  1. Prune

Allele.ctg.table 为等位基因重叠群表；Prune 步骤可以根据用于输入的 Allele.ctg.table 从 bam 文件中删除等位基因 contigs 之间的 Hi-C reads ，这这一步，仅保留折叠( collapsed )和分阶段重叠群( phased contigs )之间的最佳联系。

识别等位基因群，从而消除质量不好的 Hi-C 信号。

具体如何生成 Allele.ctg.table，请移步 ：ALLHiC: identify allelic contigs · tangerzhang/ALLHiC Wiki (github.com)

若构建简单二倍体基因组的 scaffolds，可跳过这一步。

1
2
3
4
5
6
7
8
9
10
11
12
Usage:  

ALLHiC_prune -i Allele.ctg.table -b sample.clean.bam -r draft.asm.fasta  

Parameters:   
          -h : help and usage  
          -i : Allele.ctg.table  ## 等位基因重叠群表，内含等位基因对信息
          -b : input bam file  
          -r : draft.asm.fasta  

Tips: Details on how to identify allelic contigs can be found in the following link: 
https://github.com/tangerzhang/ALLHiC/wiki/ALLHiC:-identify-allelic-contigs
  1. Partition

将 contigs 分配到预定数量的组中，基于分层链接算法迭代合并并包含它们之间大量链接的 contigs 。

此外，这里的 -K 值指的是目标物种的染色体数目，在最终结果中，每一个分组代表着一条染色体，所以我们最后得到的结果是染色体水平的 scaffolds

1
2
3
4
5
6
7
8
9
10
11
Usage:  

ALLHiC_partition -b prunning.bam -r draft.asm.fasta -e AAGCTT -k 16  

Parameters: 
      -h : help and usage.
      -b : prunned bam (optional, default prunning.bam)
      -r : draft.sam.fasta
      -e : enzyme_sites (HindIII: AAGCTT; MboI: GATC) ## 限制性酶切位点
      -k : number of groups (user defined K value)    ## 可以设置多个 K 值，从而比较出最佳结果
      -m : minimum number of restriction sites (default, 25)
  1. Rescue

由于 Hi-C 信号较弱，在 Partition 步骤中，可能得不到很多 contigs ；于是，这一步主要是为了恢复这些 contigs。

同时，进行进一步的排序和定向。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
Usage:  

ALLHiC_rescue -b sample.clean.bam -r draft.asm.fasta -c clusters.txt -i counts_RE.txt  

Parameters:  
      -h : help and usage.
      -b : sample.clean.bam (unpruned bam)
      -r : draft.sam.fasta
      -c : clusters.txt
      -i : counts_RE.txt
      -m : minimum single density for rescuing contigs (optional, default 0.01)  ## 这里设置保留 contigs 的最小阈值，默认是 0.01
Tips: 
    clusters.txt (-c) and counts_RE.txt (-i) can be generated by allhic extract.
    The format of clusters.txt is like below:
    #Group  nContigs    Contigs
    48g1    506         tig0000150 tig0000151 tig0000152 
    48g2    692         tig0000097 tig0000114 tig0000231 
    48g3    683         tig0000015 tig0000035 tig0000235
##  组别    该组含有的congtis个数	具体的 contigs 编号
    The first column is group name. The second is number of
    contigs anchored in the group and the third lists contigs 
    in the group.
  1. Optimize

优化每一个 Group 的排序，同时确定 contig 的正负方向。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
Usage:  

allhic extract sample.clean.bam draft.asm.fasta --RE AAGCTT  
allhic optimize group1.txt sample.clean.clm
allhic optimize group2.txt sample.clean.clm
...
allhic optimize group16.txt sample.clean.clm
## optimize 需要运行 -k 16 次；故这里我们可以添加一个循环，以便自动运行 16 次
groups=16
for i in {1..$gourps}; do allhic optimize group${i}.txt sample.clean.clm; done

Parameters:  

allhic extract  
    Input files:  
        bam file and contig-level assembly  ## bam 文件和 contig 水平的组装 fasta 文件
    Options:
        --RE value  Restriction site pattern (default: "GATC")  ## 限制性酶切位点

allhic optimize  
    Input files: 
      **counts_RE.txt - counts of restriction sites for each  ## 在每组中对酶切位点的数量统计
        contig in a clusterd group, which can be generated 	  ## 以及 contig 长度
        by allhic extract. The format is like below:
             #Contig     RECounts    Length
             tig0000001      572     157863
             tig0000002      143     33000
             tig0000003      231     60910
             tig0000004      3789    1044000
             tig0000005      646     166098
             tig0000006      67      15000
             tig0000007      1094    319000
      **clmfile -  The file records basic information of Hi-C  ## 该文件记录了两个 contigs 之间 Hi-C 链接的
      links between two contigs, including potential ordering  ## 基本信息，包含顺序和方向
      and orientation, number of supported reads and distance of ## 由 allhic extract 生成
      paired-end reads. This file can be accessed by allhic
      extract.
    Options:
           --skipGA       Skip GA step                   ## 跳过 GA 过程
           --resume       Resume from existing tour file ## 从现有的 tour 文件中恢复
           --seed value   Random seed (default: 42)      ## 随机种子？
           --npop value   Population size (default: 100)
           --ngen value   Number of generations for convergence (default: 5000)  ## 收敛代数
           --mutpb value  Mutation prob in GA (default: 0.2)  ## 遗传算法（GA）中的突变概率

  1. Build

将 tour 格式的文件转换成 fasta 序列文件和 agp 位置文件；从而得到染色体水平的 scaffolds --- groups.asm.fasta。

agp 文件记录每个 contig 的位置和方向。

1
2
3
4
5
6
7
Usage:  

ALLHiC_build draft.asm.fasta  

Input file:
    draft.asm.fasta - contig-level assmbly  ## contig 水平的组装

至此，我们已经从预组装的 contigs 得到最终染色体水平的 scaffolds ；但是单倍型体现在哪里呢？

在 Hi-C 质量高，预组装的 contigs 的 N50 指标高，且 -k 值的设置合理的情况下，Build 步骤输出的 groups.asm.fasta 文件就是我们最终得到的染色体水平的单倍型组装结果，也就是说我们能在 groups.asm.fasta 文件中获得我们需要的两个单倍型（二倍体）；只不过是需要通过与参考基因组绘制两者之间的点阵图来分析确定哪些染色体是同源的。

  1. Plot

我们将生成染色质接触矩阵来评估基因组支架

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
Usage:  

samtools faidx groups.asm.fasta                        ## 创建索引
cut -f1,2 groups.asm.fasta.fai|grep sample > chrn.list ## 取出长度信息

ALLHiC_plot sample.clean.bam groups.agp chrn.list 500k pdf

Input files:
   **bam file   - mapping bam file
   **agp file   - Placement of contigs in each Hi-C groups,
    which can be generated by ALLHiC_build
   **chrn.list  - a list of group name and length. The 
    format of this file is like below:  ## 第一列为 组名，第二列为对应的长度
        group1  125000
        group2  159000
    **bin size - heatmap bin size     ## 热图大小
    **ext      - extension of plot figures, e.g. pdf  ## 输出格式

以上 1-7 是 ALLHiC 提供的从预组装的 contigs 得到染色体水平的 scaffolds 的流程，但并不是表示其可以适用于任何基因型组装。

针对简单二倍体的 scaffolds 组装，教授也给出了更适合的流程：ALLHiC: scaffolding of a simple diploid genome · tangerzhang/ALLHiC Wiki (github.com)

同时，教授还对简单二倍体基因组 scaffolds 组装包装了一个脚本，可以更加方便的得到我们想要的单倍型基因组：ALLHiC_pip.sh

ALLHiC 作用于多倍体的单倍型组装

针对多倍体的单倍型组装，我们可能遇到 Hi-C 数据质量低的情况，从而导致最终用于评估 scaffolds 的热图效果不好，这是我们可以利用 ALLHIC_corrector 根据 Hi-C 信号检测并纠正错误连接的重叠群。它利用核心算法对 3D-DNA 进行重叠群校正，但不包括迭代 scaffolds 步骤，这为大基因组节省了大量时间。此外，ALLHiC_corrector 将 bam 文件（来自 bwa mem）作为输入，不需要 juice box mapping 。

教授在 Github issue 中有讨论到这个问题：haplotype assembly for a polyploid · Issue #86 · tangerzhang/ALLHiC (github.com)

参考链接

[1] Home · tangerzhang/ALLHiC Wiki (github.com)

[2] tanghaibao/allhic: Genome scaffolding based on HiC data in heterozygous and high ploidy genomes (github.com)

[3] Assembly of allele-aware, chromosomal-scale autopolyploid genomes based on Hi-C data | Nature Plants

[4] ALLHiC: scaffolding of a simple diploid genome · tangerzhang/ALLHiC Wiki (github.com)

[5] ALLHiC: scaffolding an auto polyploid sugarcane genome · tangerzhang/ALLHiC Wiki (github.com)

[6] Question about pruning step and input/output files · Issue #94 · tangerzhang/ALLHiC (github.com)

[7] haplotype assembly for a polyploid · Issue #86 · tangerzhang/ALLHiC (github.com)