前言

在遺傳學(xué)中，二倍體常被視為“標(biāo)準(zhǔn)模型”，也是當(dāng)前大多數(shù)基因組組裝方法的默認(rèn)假設(shè)。然而在自然界，尤其是植物中，多倍體并不罕見(jiàn)。事實(shí)上，許多重要作物——如小麥、棉花、油菜、甘薯、馬鈴薯和甘蔗等——本身就是多倍體。相較于二倍體，多倍體通常具有更強(qiáng)的環(huán)境適應(yīng)能力和進(jìn)化潛力，這也是它們?cè)谧匀贿x擇和人工馴化過(guò)程中被反復(fù)保留下來(lái)的重要原因。 但從基因組角度看，多倍體意味著更多高度相似的染色體拷貝和更復(fù)雜的結(jié)構(gòu)關(guān)系，這也使其基因組組裝難度遠(yuǎn)高于二倍體。因此，在理解了二倍體基因組組裝的基本原理之后，今天繼續(xù)來(lái)介紹多倍體基因組的組裝的流程，以及組裝過(guò)程中存在的難點(diǎn)。

為什么多倍體基因組難拼？

要理解多倍體基因組為什么“難拼”，首先需要回到一個(gè)最根本的問(wèn)題：多倍體是如何形成的。

多倍體是如何產(chǎn)生的？

多倍體通常起源于全基因組加倍（Whole-Genome Duplication, WGD）事件，即原本的一整套染色體被整體復(fù)制了一次或多次。根據(jù)參與加倍的基因組來(lái)源不同，多倍體可以大致分為兩類：同源多倍體和異源多倍體。

1. 同源多倍體（autopolyploid） 來(lái)源于同一物種內(nèi)部的染色體加倍，例如二倍體個(gè)體發(fā)生減數(shù)分裂異常，產(chǎn)生了四倍體后代。這類多倍體中，各套染色體之間序列高度相似，幾乎沒(méi)有清晰的“亞基因組”邊界。

<center>圖 1. 同源多倍體的起源模式圖</center>

2. 異源多倍體（allopolyploid） 則來(lái)源于不同物種之間的雜交，隨后發(fā)生染色體加倍。此時(shí)，多倍體基因組中往往可以分辨出來(lái)自不同祖先物種的亞基因組（subgenome），但這些亞基因組之間仍然保持著較高的同源性。

<center>圖 2. 異源多倍體的起源模式圖</center>

難在哪里？

無(wú)論是同源多倍體還是異源多倍體，它們?cè)诨蚪M組裝中都會(huì)帶來(lái)一個(gè)共同的問(wèn)題：高度相似的多套染色體同時(shí)存在。這直接挑戰(zhàn)了基因組組裝中最核心的一條假設(shè)，即“一段序列在基因組中只有一個(gè)正確位置”。

具體來(lái)說(shuō)，多倍體基因組的組裝難點(diǎn)主要體現(xiàn)在三個(gè)方面：

1. 同源染色體難以區(qū)分

• 多套染色體在序列層面高度相似，尤其是在同源多倍體中，組裝算法很難判斷一條 read 究竟來(lái)自哪一套染色體，極易發(fā)生錯(cuò)誤合并或錯(cuò)誤拆分。

2. 重復(fù)與雜合信號(hào)高度耦合

• 在二倍體中，“重復(fù)”和“雜合”尚可區(qū)分；而在多倍體中，多拷貝同源序列、亞基因組差異和真實(shí)重復(fù)序列往往疊加在一起，使得 k-mer 分布、覆蓋深度和組裝圖結(jié)構(gòu)變得極其復(fù)雜。

這些問(wèn)題會(huì)導(dǎo)致組裝的基因組中存在：

• 兩條或多條高度相似的同源序列（allelic contigs）。這些本該分開(kāi)的序列被錯(cuò)誤地合并成了一條 contig，就稱為 坍塌（collapse）
• 亞基因組被混合
• 單倍型或亞基因組結(jié)構(gòu)被破壞
  
  

其中同源染色體被錯(cuò)誤折疊形成的坍塌區(qū)是最常出現(xiàn)，也是最難解決的問(wèn)題。

怎么判斷組裝的基因組是否存在坍塌區(qū)？

通過(guò)HIC信號(hào)圖判斷

這里引用曾老師文章中的圖來(lái)介紹塌區(qū)的形成原因和幾種不同的坍塌區(qū)情況。

1. Ideal situation

• 左右兩條是真實(shí)的同源染色體
• Ctg1 ? Ctg3 是一對(duì)等位區(qū)域
• Ctg2 ? Ctg4 是一對(duì)等位區(qū)域

最終，每一條 contig 都有一個(gè)清晰的“等位搭檔”，不會(huì)形成坍塌區(qū)。

2.Ambiguous / incorrect mapping

• Hi-C 信號(hào)開(kāi)始“變花”
• 出現(xiàn)一些額外的非對(duì)角信號(hào)
• 對(duì)角線結(jié)構(gòu)仍然存在

雖然來(lái)自于不同同源染色體的contig會(huì)出現(xiàn)信號(hào)的混亂，可能會(huì)把不同的源染色體里contig組裝到一起，但是并不會(huì)出現(xiàn)某兩個(gè)contig的合并，因此不會(huì)出現(xiàn)坍塌。

3. Collapsed contigs

• 2個(gè)或多個(gè)contig合并成了1個(gè)（坍塌）
• 坍塌的contig與兩個(gè)相鄰 contig 都出現(xiàn)強(qiáng)hic信號(hào)

來(lái)自不同同源染色體的contig因?yàn)楦叨认嗨疲M裝的過(guò)程中合并成了一條contig，直接導(dǎo)致遺傳信息的丟失。

4. Chimeric contigs

• contig 由不該連在一起的兩段contig拼成
• Hi-C 信號(hào)在 contig 內(nèi)部出現(xiàn)斷裂
• 熱圖中的信號(hào)被“切開(kāi)”，左右部分分別對(duì)應(yīng)不同 contig

來(lái)自于兩個(gè)同源染色體，高度相似的contig拼到了一起，不會(huì)出現(xiàn)遺傳信息的缺失，剪開(kāi)重新分配即可。

<center>圖 3. 坍塌區(qū)的形成原因和幾種不同的坍塌區(qū)</center>

通過(guò)測(cè)序深度判斷

在前面的內(nèi)容中，我們已經(jīng)討論了如何通過(guò) Hi-C 圖識(shí)別坍塌區(qū)。但Hi-C 僅提供了圖的空間信息，并且對(duì)測(cè)序深度的變化并不敏感。因此，我們還可以通過(guò)測(cè)序深度來(lái)輔助判斷是否存在基因組坍塌現(xiàn)象，尤其是在多倍體基因組中，這一方法尤為重要。正常情況下，每條染色體或亞基因組的測(cè)序深度應(yīng)該是相對(duì)穩(wěn)定的。例如，二倍體基因組的每條染色體應(yīng)當(dāng)有兩份拷貝，因此如果測(cè)序覆蓋度設(shè)置為 30×，則每個(gè)位置的覆蓋深度應(yīng)該大致在 30× 附近。

但當(dāng)發(fā)生基因組坍塌時(shí)，我們看到的測(cè)序深度會(huì)有以下幾種變化：

1. 同源染色體的測(cè)序深度變化

• 正常情況：同源染色體會(huì)有相似的測(cè)序深度，通常是相對(duì)均衡的。
• 坍塌情況：如果多條同源染色體的區(qū)域被錯(cuò)誤地拼接到同一條 contig 上，那么該區(qū)域的測(cè)序深度會(huì)大大增加，因?yàn)橄嗤膮^(qū)域被多次測(cè)序。

2. 區(qū)域的測(cè)序深度偏高

• 正常情況：基因組的不同區(qū)域因基因密度、GC 含量等因素可能略有起伏，但整體變化不大。
• 坍塌情況：坍塌區(qū)域的深度顯著高于其他區(qū)域，因?yàn)檫@部分區(qū)域由原本應(yīng)分開(kāi)測(cè)序的兩條同源染色體組成，導(dǎo)致該區(qū)域的重復(fù)序列被重復(fù)計(jì)數(shù)，測(cè)序深度看起來(lái)異常。

一般來(lái)講，如果是兩個(gè)同源染色體發(fā)生了坍塌，坍塌的contig或者染色體的Depth會(huì)變成常規(guī)的2倍，如果是三個(gè)染色體發(fā)生了坍塌，則是3倍，依次類推。

<center>圖 4. 坍塌區(qū)的HIC熱圖和測(cè)序深度直方圖</center>

那么如果出現(xiàn)了坍塌區(qū)，我們應(yīng)該怎么處理呢？
很簡(jiǎn)單，因?yàn)閏ontig高度相似才發(fā)生坍塌，所以直接把坍塌的contig重新復(fù)制一份，然后手動(dòng)添加到對(duì)應(yīng)的同源染色體區(qū)域中。但是這樣也不完全對(duì)，因?yàn)榫退銉蓚€(gè)contig再相似，也不可能完全一樣，但是可能是目前最好的解決辦法了。此外，如果采用上述的辦法，HIC熱圖是沒(méi)有任何改善的，還是一眼就能看出來(lái)坍塌區(qū)的存在，就算是拷貝了一份坍塌contig，重新掛載，也會(huì)因?yàn)閔ic reads 的多重比對(duì)直接導(dǎo)致兩個(gè)區(qū)域直接沒(méi)有任何信號(hào)。那么有沒(méi)有辦法既能復(fù)制contig，還能改善HIC熱圖呢？有的，有的兄弟，先賣個(gè)關(guān)子，后面會(huì)介紹！

多倍體基因組組裝實(shí)操

背景介紹可能有點(diǎn)多了，大家可能有些不耐煩了，話不多說(shuō)，直接開(kāi)始基因組的組裝。同樣還是使用hifiasm，前面介紹過(guò)了，可以直接參考過(guò)去的推文。

hifiasm -o sample.asm -t 32 \
  --h1 hic_R1.fastq.gz \
  --h2 hic_R2.fastq.gz \
  sample.hifi.fastq.gz
  
# -o 輸出文件前綴
# -t 運(yùn)行的線程數(shù)
# --h1 --h2 輸入的Hic reads

和二倍體基因組的組裝不同，因?yàn)閜_ctg.gfa是經(jīng)過(guò)一系列過(guò)濾的從而得到的高質(zhì)量的主路徑基因組，所以是有損的。因此在多倍體基因組組裝中，推薦使用無(wú)損的p_utg.gfa文件。雖然這樣損失了基因組質(zhì)量，但是可以保留更多的信息。

軟件下載

推薦使用Haphic 進(jìn)行染色體掛載，本人實(shí)測(cè)在應(yīng)對(duì)多倍體的表現(xiàn)很好。

# 下載Haphic
git clone https://github.com/zengxiaofei/HapHiC.git
# 通過(guò)conda的yml文件安裝依賴環(huán)境
conda env create -f HapHiC/conda_env/environment_py310.yml
# 啟動(dòng)Haphic環(huán)境
conda activate haphic # or: source /path/to/conda/bin/activate haphic
# 檢查軟件依賴
/path/to/HapHiC/haphic check
# 測(cè)試打印軟件幫助信息
/path/to/HapHiC/haphic -h

染色體掛載

將Hi-C測(cè)序數(shù)據(jù)比對(duì)到組裝好的基因組上，以構(gòu)建BAM文件。

# 比對(duì)
bwa index asm.fa
bwa mem -5SP -t 28 asm.fa /path/to/read1_fq.gz /path/to/read2_fq.gz | samblaster | samtools view – -@ 14 -S -h -b -F 3340 -o HiC.bam

# 通過(guò)MAPQ >1 和NM < 3篩選比對(duì)結(jié)果
/path/to/HapHiC/utils/filter_bam HiC.bam 1 –nm 3 –threads 14 | samtools view – -b -@ 14 -o HiC.filtered.bam

運(yùn)行Haphic pipline進(jìn)行染色體掛載

/path/to/HapHiC/haphic pipeline asm.fa \
HiC.filtered.bam \
nchrs \
–RE “AAGCTT” \
–threads 8

# nchr 改成預(yù)估的染色體數(shù)目
# –RE 是酶切位點(diǎn)的類型

主要的輸出結(jié)果：

• 01.cluster/corrected_asm.fa：更正后的 FASTA 格式的程序集。僅在啟用裝配體校正時(shí)生成此文件。
• 4.build/scaffolds.agp：一個(gè) SALSA 風(fēng)格的 AGP 文件，包含有關(guān) corrected_asm.fa 中所有序列的腳手架分配、排序和方向信息。如果有任何嵌合重疊群被 HapHiC 校正，損壞的重疊群將被分配新的 ID。（后面的構(gòu)建掛載圖和Juicebox可視化調(diào)整需要用到的結(jié)果）
• 04.build/scaffolds.raw.agp：一個(gè) YaHS 風(fēng)格的 AGP 文件，包含有關(guān) asm.fa 中所有序列的腳手架分配、排序和方向信息。不會(huì)為斷開(kāi)的重疊群分配新的 ID，但它們?cè)谠贾丿B群中的起始和結(jié)束坐標(biāo)將顯示在第七列和第八列中。
• 04.build/scaffolds.fa：FASTA 格式的最終腳手架。（染色體掛載后的scaffold水平的基因組）
• 04.build/juicebox.sh：用于 Juicebox 可視化和管理的 shell 腳本。（快速生成Juicebox的輸入文件腳本）

Juicebox 手動(dòng)調(diào)整

# 生成 .hic和.assembly文件
bash juicebox.sh

調(diào)整后，生成最終的基因組

/path/to/HapHiC/utils/juicer post \
-o out_JBAT \
out_JBAT.review.assembly \
out_JBAT.liftover.agp asm.fa

繪制HIC熱圖

/path/to/HapHiC/haphic plot \
out_JBAT.FINAL.agp \
HiC.filtered.bam
--min_len 10 \
--threads 20

# --min_len 繪制熱圖的sacffod最小長(zhǎng)度，可以通過(guò)調(diào)整長(zhǎng)度來(lái)設(shè)置只繪制染色體水平的sacffod
# --threads 運(yùn)行的線程數(shù)

如果沒(méi)有發(fā)現(xiàn)坍塌區(qū)，基因組的組裝到這就結(jié)束了，后面是演示如何檢測(cè)坍塌區(qū)。

檢測(cè)坍塌區(qū)

這里我們選擇通過(guò)把測(cè)序的reads比對(duì)回基因組，然后計(jì)算每條染色體的Depth。

minimap2 -t 60 -ax map-ont out_JBAT.FINAL.fa YY-2.pass.all.fq.gz|samtools view -bS -@ 60 - | samtools sort -@ 60 -o aln.sorted.bam
samtools index aln.sorted.bam

# 計(jì)算測(cè)序深度
mosdepth -t 10 \
--by 100000 \
100000_genome_depth aln.sorted.bam

生成一個(gè)后綴為“mosdepth.summary.txt”文件，有每條染色體的長(zhǎng)度和深度。以我的基因組為例，已知基因組的測(cè)序深度約為21×，可以看到Chr_3有8條同源染色體，每條染色體的深度都在21左右，這就是正常的。而Chr_1只有6條同源染色體，并且Chr_1_A、Chr_1_B和Chr_1_C的深度都明顯大于21，說(shuō)明這三條染色體肯定都存在坍塌的區(qū)域。此外，因?yàn)榘l(fā)生了坍塌，每條同染色體的長(zhǎng)度肯定也會(huì)有很大的差別。

接下來(lái)通過(guò)我自己寫(xiě)的腳本，按照滑窗繪制染色體每個(gè)區(qū)間的測(cè)序深度

# 繪制深度覆蓋圖
python plot_depth_along_chromosome.py \
-i 100000_genome_depth.regions.bed \
-c scaffold_10 \
-o scaffold_10.png

# -i mosdepth計(jì)算的測(cè)序深度
# -c 指定繪制的染色體
# -o 輸出文件

可以看到這個(gè)條染色體發(fā)生了超大面積的坍塌，大約從1000-17000kb都是坍塌區(qū)，深度在42×左右，因此可以推斷可能是2同源染色體坍塌成了一個(gè)區(qū)域。

這條染色體更特殊，一部分區(qū)域的深度達(dá)到了60×，因此是由三個(gè)同源染色體坍塌形成的。

在這里整理一下坍塌區(qū)的幾個(gè)金標(biāo)準(zhǔn)

• 測(cè)序深度是基因組覆蓋度的幾倍
• 同源染色體長(zhǎng)度相差很多
• HIC信號(hào)復(fù)雜，坍塌的contig回同時(shí)和多個(gè)contig具有很強(qiáng)的信號(hào)。

結(jié)語(yǔ)

其實(shí)本來(lái)這篇文章想把怎么處理坍塌區(qū)的部分也加上的，但是篇幅太多了，所以后面會(huì)專門寫(xiě)一篇處理多倍體組裝中坍塌區(qū)的文章。

本文由mdnice多平臺(tái)發(fā)布