使用purge haplogs處理基因組雜合區域

2021-09-29 03:43:58 字數 3916 閱讀 8025

falcon和canu的組裝後會得到乙個單倍型融合的基因組，用來表示二倍體基因組。之後，falcon unzip和supernova這類軟體進一步處理其中等位基因區域，將這部分區間進行拆分。

當基因組某些區域可能有著比較高的雜合度，這會導致基因組該區域的兩個單倍型被分別組裝成primary contig， 而不是乙個為primary contig， 另乙個是associated haplotig. 如果下游分析主要關注於單倍型，這就會導致一些問題。

那麼有沒有解決方案呢？其實也很好辦，就是找到相似度很高的contig，將他們拆分。目前已經有一些軟體可以完成類似任務，如haplomerger2,redundans, 這不過這些軟體主要處理二代組裝結果。

purge_haplogs則是最近開發，用於三代組裝的基因組。它根據minimap2的比對結果，通過分析比對read的覆蓋度決定誰去誰留。該工具適用於單倍型組裝軟體，例如 canu, falcon或 falcon-unzip primary contigs, 或者是分相後的二倍體組裝(falcon-unzip primary contigs + haplotigs 。

它的工作流程如下圖所示。一般只需要兩個輸入檔案，組裝草圖(fasta格式) 和 比對的bam檔案。同時還可以提供重複序列注釋的bed檔案，輔助處理高重複的contig。

分析流程

建議: 用原來用於組裝的read進行比對。對於多個匹配的read，建議採取random best，也就是隨便找乙個。

purge_haplotigs依賴軟體比較多，手動安裝會很麻煩，但是他可以直接用bioconda裝

conda create -n purge_haplotigs_env


conda activate purge_haplotigs_env


conda install purge_haplotigs
安裝完成後需要一步測試

purge_haplotigs test
mkdir purge_haplotigs_tutorial


cd purge_haplotigs_tutorial


wget 


wget # 1.7g


wget .bai 


wget 


wget .fai
當然我們不可能直接就拿到比對好的bam檔案，我們一般是有組裝後的基因組以及用於組裝的subread，假設這兩個檔案命名為, genome.fa 和 subreads.fasta.gz.

minimap2 -ax map-pb genome.fa subreads.fasta.gz \


| samtools view -hf 256 - \


| samtools sort -@ 8 -m 1g -o aligned.bam -t tmp.ali
第一步：使用purge_haplotigs readhist從bam中統計read深度，繪製柱狀圖。

purge_haplotigs  readhist  -b aligned.bam  -g genome.fasta  -t 20


# -t 執行緒數, 不宜過高，過高反而沒有效果。
也就是下圖，你能明顯的看到圖中有兩個峰，乙個是單倍型的覆蓋度，另乙個二倍型的覆蓋度，

高雜合基因組read-depth histogram

你可能還想知道高純合基因組是什麼樣的效果，我也找了乙個純合的物種做了也做了read-depth 柱狀圖，

純合基因組read-depth histogram

之後你需要根據read-depth 柱狀圖 確定這兩個峰的位置用於下一步。下面是兩個例子。對於我們則是，20，65，190.

兩個例子

第二步: 根據read-depth資訊選擇閾值。

purge_haplotigs  contigcov  -i cns_p_ctg.aligned.sd.bam.gencov  -o coverage_stats.csv  -l 20  -m 75  -h 190
這一步生成的檔案是"coverage_stats.csv"

第三步：區分haplotigs.

purge_haplotigs purge  -g cns_p_ctg.fasta  -c coverage_stats.csv  -b cns_p_ctg.aligned.sd.bam  -t 4  -a 60
這一步會得到如下檔案

000004f,primary -> 000004f_013,haplotig


-> 000004f_017,haplotig 


-> 000004f_004,haplotig


-> 000004f_027,haplotig
在"curated.reassignments.tsv"檔案中有6列

由於我們用的是單倍型組裝primary contigs而不是二倍體組裝的parimary + haplotigs, 因此我們需要將falcon_unzip的haplotgi合併到重新分配的haplotigs中，這樣子我們依舊擁有二倍體組裝結果

cat cns_h_ctg.fasta >> curated.haplotigs.fasta
為什麼第一步的 read 覆蓋深度分析能判斷基因組是否冗餘呢？這是因為對於坍縮的單倍型，那麼含有等位的基因的read只能比對到該位置上，而如果雜合度太高被拆分成兩個不同的contig，那麼含有等位的基因的read就會分別比對到不同的read上，導致深度降低一半。下圖a就是乙個典型的包含冗餘基因組的read覆蓋度分布

read-深度分析

分析流程的第二步的任務就是人工劃分出如下圖b部分，綠色的部分是坍縮單倍型contig，藍色的部分是潛在的冗餘contig。之後，分析流程會計算這些區域中的contig的覆蓋度。 對於綠色部分中的contig，如果覆蓋度低於80%, 會進行標記用於後續分析。如果深度非常低，那麼很有可能就是組裝引入錯誤，深度非常高的部分基本就是重複序列或者是細胞器的contig，這些黃色的contig可以在後續的組裝出分開。

劃分區間

第三步就是同源序列進行識別和分配。所有標記的contig之後會用minimap2在整個組裝進行搜尋，尋找相似度較高的離散區間（如下圖c）。如果乙個contig的聯配得分大於閾值(預設70%), 那麼就會被標記為haplotigs. 如果乙個contig的最大聯配得分大於閾值(預設250%), 會被標記成重複序列，這有可能是潛在的有問題contig，或許是坍縮的contig或者低複雜度序列。

移除haplotigs

在Java SE中使用Hibernate處理資料

hibernate.dialect net.sf.hibernate.dialect.postgresqldialect hibernate.connection.driver class org.postgresql.driver hibernate.connection.url jdbc pos...

pytorch使用Dataset分批次處理資料

import torch import numpy as np from torch.utils.data import dataset from torch.utils.data import dataloader import matplotlib.pyplot as plt prepare d...

在android中使用OkHttp框架處理網路請求

okhttp網路處理框架，分成下面幾個使用過程 同步get 非同步get 不要超過1m 大小沒有限制 public void run throws exceptiond.post a form,簡單鍵值對 public void run throws exceptione.post混合引數型別 st...