De novo組裝#05 | 基因組糾錯(cuò)polish（racon+pilon, nextpolish）

寫在前面

在進(jìn)行染色體掛載之前一般對(duì)得到的 "primary assembly"（主要組裝/初始組裝結(jié)果）進(jìn)行進(jìn)一步的優(yōu)化，以減少錯(cuò)誤，提高基因組的質(zhì)量和準(zhǔn)確性。然后再將這些相對(duì)碎片化但polish過的高質(zhì)量primary assembly進(jìn)行染色體掛載以得到一個(gè)染色體水平的組裝結(jié)果。如果用CLR數(shù)據(jù)組裝出的primary assembly在沒有進(jìn)行polish的情況下，直接用HIC數(shù)據(jù)進(jìn)行染色體掛載會(huì)有如下問題：（以下幾點(diǎn)回答來自gpt3.5）

  1. 錯(cuò)誤和噪聲： CLR（Continuous Long Reads，如 PacBio 或 Nanopore 數(shù)據(jù)）通常具有較高的錯(cuò)誤率，可能會(huì)導(dǎo)致組裝的主要序列中存在錯(cuò)誤和噪聲。直接使用這些序列進(jìn)行 Hi-C 染色體掛載可能會(huì)將這些錯(cuò)誤傳播到聯(lián)系矩陣中，影響到染色體的準(zhǔn)確性。

  2. 聯(lián)系矩陣質(zhì)量下降： 染色體掛載的關(guān)鍵是基于 Hi-C 數(shù)據(jù)生成染色體之間的聯(lián)系矩陣。如果主要序列存在錯(cuò)誤，這些錯(cuò)誤可能會(huì)影響到聯(lián)系矩陣的質(zhì)量，從而降低染色體掛載的準(zhǔn)確性。

  3. 組裝斷裂： 主要序列中的錯(cuò)誤可能會(huì)導(dǎo)致組裝斷裂，使得 Hi-C 數(shù)據(jù)在某些區(qū)域無法正確匹配。這會(huì)導(dǎo)致染色體掛載失敗或出現(xiàn)錯(cuò)誤。

  4. 重復(fù)序列挑戰(zhàn)： 高復(fù)雜性和重復(fù)性區(qū)域在 de novo 組裝中是一個(gè)挑戰(zhàn)，而 CLR 數(shù)據(jù)由于錯(cuò)誤率較高，可能在這些區(qū)域中表現(xiàn)不佳。這會(huì)導(dǎo)致在 Hi-C 染色體掛載中難以準(zhǔn)確匹配和定位這些區(qū)域。

因此，為了獲得更準(zhǔn)確的染色體掛載結(jié)果，建議在進(jìn)行 Hi-C 染色體掛載之前，對(duì)主要序列進(jìn)行糾錯(cuò)（polishing），以減少錯(cuò)誤和噪聲。糾錯(cuò)可以使用高質(zhì)量的測(cè)序數(shù)據(jù)（如 illumina 數(shù)據(jù)）對(duì)主要序列進(jìn)行校正，提高組裝的質(zhì)量和準(zhǔn)確性。這樣，在染色體掛載過程中，你會(huì)使用更可靠的主要序列和聯(lián)系矩陣，從而獲得更準(zhǔn)確的染色體邊界和結(jié)構(gòu)信息。

另外，我們可以同時(shí)用三代和二代數(shù)據(jù)一起對(duì)我們的初始組裝進(jìn)行糾錯(cuò)。三代數(shù)據(jù)糾錯(cuò)主要用于解決長重復(fù)區(qū)域和復(fù)雜結(jié)構(gòu)，提高基因組的連續(xù)性和準(zhǔn)確性，有助于捕獲染色體級(jí)別的信息；而二代數(shù)據(jù)糾錯(cuò)主要用于在主要序列中糾正小錯(cuò)誤，提高基因組的局部準(zhǔn)確性，從而使基因組序列精準(zhǔn)度更接近于真實(shí)的 DNA 序列。綜合而言，使用三代+二代數(shù)據(jù)進(jìn)行糾錯(cuò)是一種綜合優(yōu)化的策略，有助于克服不同測(cè)序技術(shù)的優(yōu)缺點(diǎn)，提高基因組的整體準(zhǔn)確性和質(zhì)量。這兩種數(shù)據(jù)的結(jié)合能夠產(chǎn)生更好的組裝和糾錯(cuò)效果，為后續(xù)的基因組研究提供更可靠的基因組參考。

這里我們我們用了兩種策略：

  • 策略1：racon+pilon，即用racon軟件進(jìn)行三代數(shù)據(jù)糾錯(cuò)，再用pilon軟件進(jìn)行二代數(shù)據(jù)糾錯(cuò)；
  • 策略2：直接使用nextpolish軟件進(jìn)行三代和二代數(shù)據(jù)同時(shí)糾錯(cuò)。

數(shù)據(jù)準(zhǔn)備

1.primary assembly初始組裝文件，這里用flye組裝的初始組裝做演示。
1.用于的組裝的pacbio CLR三代測(cè)序文件。
2.經(jīng)過fastp質(zhì)控過的illumina二代測(cè)序文件。

polish策略1

racon+pilon，即用racon軟件進(jìn)行三代數(shù)據(jù)糾錯(cuò)，再用pilon軟件進(jìn)行二代數(shù)據(jù)糾錯(cuò)。

1# 軟件安裝
racon安裝
軟件主頁：https://github.com/isovic/racon

## 自動(dòng)安裝
conda install -c bioconda racon

## 手動(dòng)安裝
git clone --recursive https://github.com/lbcb-sci/racon.git racon
cd racon
mkdir build
cd build
cmake -DCMAKE_BUILD_TYPE=Release ..
make

pilon安裝
軟件主頁：https://github.com/broadinstitute/pilon

## 由于是個(gè)java程序直接下載即可使用
wget https://github.com/broadinstitute/pilon/releases/download/v1.24/pilon-1.24.jar 
## java程序統(tǒng)一使用方法：java -jar 后面接java程序
java -jar pilon-1.24.jar --help  # 查看pilon用法

其他軟件依賴：因?yàn)樯婕暗奖葘?duì)，所以需要安裝三代和二代比對(duì)軟件minimap2和bwa
minimap2主頁：https://github.com/lh3/minimap2
bwa主頁：https://github.com/lh3/bwa

## minimap2自動(dòng)安裝
conda install -c bioconda minimap2
## minimap2手動(dòng)安裝
wget https://github.com/lh3/minimap2/releases/download/v2.26/minimap2-2.26_x64-linux.tar.bz2
tar -jxvf minimap2-2.26_x64-linux.tar.bz2
./minimap2-2.26_x64-linux/minimap2
---------------------------------------------------------------------------------------------------------------------------------
## bwa手動(dòng)安裝
git clone https://github.com/lh3/bwa.git
cd bwa
make

2# 運(yùn)行程序
我這里我打算先用racon進(jìn)行3輪三代數(shù)據(jù)糾錯(cuò)，再用pilon對(duì)racon的最后一輪結(jié)果再進(jìn)行3輪二代數(shù)據(jù)糾錯(cuò)，加起來相當(dāng)于6輪了。

  • 用racon進(jìn)行3輪三代數(shù)據(jù)糾錯(cuò)
## 第1輪糾錯(cuò)，這里詳細(xì)注釋了，后兩輪就不展開了
/newlustre/home/jfgui/wy/anaconda3/envs/pengfang/bin/minimap2 \  # 用minimap2比對(duì)
      -ax map-pb \  # 用-a參數(shù)生成sam文件而非默認(rèn)的paf格式，-x參數(shù)為比對(duì)模式，這里是pacbio數(shù)據(jù)比對(duì)到參考基因組
      -t 24 \ # 線程數(shù)
      ../../purge_dups/flye_purge/assembly.fasta \ #  參考基因組，這里第1輪就是用的flye組裝出的初始組裝，后面就都是前一輪的糾錯(cuò)結(jié)果。
      ../../P01TYD20308306-1_r64030_20201110_065049_1_A02.subreads.fasta \ # 用于糾錯(cuò)三代測(cè)序文件 
      | gzip -c - > minimap1.sam.gz  # 壓縮后輸出為 minimap1.sam.gz
/newlustre/home/jfgui/wy/anaconda3/envs/pengfang/bin/racon \ # 用racon糾錯(cuò)
      -t 24 \ # 線程數(shù)
      -u \ # 糾錯(cuò)過和未被就錯(cuò)的contig都會(huì)被輸出，否則只輸出糾錯(cuò)過的contig
      ../../P01TYD20308306-1_r64030_20201110_065049_1_A02.subreads.fasta \ # 用于糾錯(cuò)三代測(cè)序文件
      minimap1.sam.gz \ # 前面的比對(duì)結(jié)果
      ../../purge_dups/flye_purge/assembly.fasta > flye.racon1.fasta # 糾錯(cuò)前的基因組指向糾錯(cuò)后的基因組

## 第2輪糾錯(cuò)
/newlustre/home/jfgui/wy/anaconda3/envs/pengfang/bin/minimap2 -ax map-pb -t 24 flye.racon1.fasta ../../P01TYD20308306-1_r64030_20201110_065049_1_A02.subreads.fasta | gzip -c - > minimap2.sam.gz
/newlustre/home/jfgui/wy/anaconda3/envs/pengfang/bin/racon -t 24 -u ../../P01TYD20308306-1_r64030_20201110_065049_1_A02.subreads.fasta minimap2.sam.gz flye.racon1.fasta > flye.racon2.fasta

## 第3輪糾錯(cuò)
/newlustre/home/jfgui/wy/anaconda3/envs/pengfang/bin/minimap2 -ax map-pb -t 24 flye.racon2.fasta ../../P01TYD20308306-1_r64030_20201110_065049_1_A02.subreads.fasta | gzip -c - > minimap3.sam.gz
/newlustre/home/jfgui/wy/anaconda3/envs/pengfang/bin/racon -t 24 -u ../../P01TYD20308306-1_r64030_20201110_065049_1_A02.subreads.fasta minimap3.sam.gz flye.racon2.fasta > flye.racon3.fasta
  • 再用pilon進(jìn)行3輪二代數(shù)據(jù)糾錯(cuò)
## 第1輪糾錯(cuò)，這里詳細(xì)注釋了，同理后兩輪就不展開了
ln -s ../flye.racon3.fasta  ./racon3.fasta  # 將需要糾錯(cuò)的基因組鏈接到到當(dāng)前工作路徑
/newlustre/home/jfgui/Wangtao/software/bwa/bwa index racon3.fasta    # 對(duì)基因組構(gòu)建index以便后面的比對(duì)
/newlustre/home/jfgui/Wangtao/software/bwa/bwa mem \ # 用bwa的mem比對(duì)模式對(duì)二代數(shù)據(jù)進(jìn)行比對(duì)
    -t 24 \ #線程數(shù)
    racon3.fasta \ # 對(duì)比用的參考基因組，即前一輪的糾錯(cuò)結(jié)果
    ../../../../../XX_survey/fastp.out/liver.clean.CLEAN.R1.fastq.gz ../../../../../XX_survey/fastp.out/liver.clean.CLEAN.R2.fastq.gz \ # 用于糾錯(cuò)的二代數(shù)據(jù)
    | /newlustre/home/jfgui/wy/anaconda3/envs/pengfang/bin/samtools sort \ # 用samtools的sort命令對(duì)比對(duì)結(jié)果進(jìn)行排序，用管道符連接 
           -@ 24 \ # samtools的線程數(shù)
           -  \ # 表示管道符前面命令的輸出作為管道符后面命令的輸入，這里也就是指未排序的bam文件
           -o bwamem1.bam # 輸出比對(duì)結(jié)果
/newlustre/home/jfgui/wy/anaconda3/envs/pengfang/bin/samtools index bwamem1.bam # 用samtools對(duì)比對(duì)文件進(jìn)行index
java -Xmx128G -jar /newlustre/home/jfgui/Wangtao/software/pilon-1.24.jar \  # 用pilon糾錯(cuò)，-Xmx128G指的是用128G內(nèi)存
    --genome racon3.fasta \  # 前一輪糾錯(cuò)的基因組
    --frags bwamem1.bam \   # 比對(duì)結(jié)果文件
    --changes \ # 糾錯(cuò)前后變化記錄文件
    --diploid \ # 二倍體，不設(shè)置會(huì)默認(rèn)單倍體進(jìn)而識(shí)別單倍型
    --threads 24 \   # 線程數(shù)
    --outdir ./pilon1.out \  # 輸出目錄 
    --output racon3.pilon1   # 輸出前綴

## 第2輪糾錯(cuò)
ln -s ./pilon1.out/racon3.pilon1.fasta  ./racon3.pilon1.fasta
/newlustre/home/jfgui/Wangtao/software/bwa/bwa index racon3.pilon1.fasta
/newlustre/home/jfgui/Wangtao/software/bwa/bwa mem -t 24 racon3.pilon1.fasta ../../../../../XX_survey/fastp.out/liver.clean.CLEAN.R1.fastq.gz ../../../../../XX_survey/fastp.out/liver.clean.CLEAN.R2.fastq.gz | /newlustre/home/jfgui/wy/anaconda3/envs/pengfang/bin/samtools sort -@ 24 - -o bwamem2.bam
/newlustre/home/jfgui/wy/anaconda3/envs/pengfang/bin/samtools index bwamem2.bam
java -Xmx128G -jar /newlustre/home/jfgui/Wangtao/software/pilon-1.24.jar --genome racon3.pilon1.fasta --frags bwamem2.bam --changes  --diploid --threads 24 --outdir ./pilon2.out --output racon3.pilon2

## 第3輪糾錯(cuò)
ln -s ./pilon2.out/racon3.pilon2.fasta  ./racon3.pilon2.fasta
/newlustre/home/jfgui/Wangtao/software/bwa/bwa index racon3.pilon2.fasta
/newlustre/home/jfgui/Wangtao/software/bwa/bwa mem -t 24 racon3.pilon2.fasta ../../../../../XX_survey/fastp.out/liver.clean.CLEAN.R1.fastq.gz ../../../../../XX_survey/fastp.out/liver.clean.CLEAN.R2.fastq.gz | /newlustre/home/jfgui/wy/anaconda3/envs/pengfang/bin/samtools sort -@ 24 - -o bwamem3.bam
/newlustre/home/jfgui/wy/anaconda3/envs/pengfang/bin/samtools index bwamem3.bam
java -Xmx128G -jar /newlustre/home/jfgui/Wangtao/software/pilon-1.24.jar --genome racon3.pilon2.fasta --frags bwamem3.bam --changes  --diploid --threads 24 --outdir ./pilon3.out --output racon3.pilon3

3# 結(jié)果查看

3輪racon結(jié)果

3輪racon糾錯(cuò)后其結(jié)果還是比較簡潔的，就3個(gè)比對(duì)文件3個(gè)糾錯(cuò)后基因組文件。
assembly-stats查看，基因組大小是依次變大了。

$assembly-stats flye.racon1.fasta flye.racon2.fasta flye.racon3.fasta
stats for flye.racon1.fasta
sum = 774215185, n = 680, ave = 1138551.74, largest = 25240035
N50 = 9242029, n = 27
N60 = 7581968, n = 36
N70 = 6421899, n = 47
N80 = 5066544, n = 61
N90 = 2268818, n = 81
N100 = 517, n = 680
N_count = 0
Gaps = 0
-------------------------------------------------------------------------------
stats for flye.racon2.fasta
sum = 775085635, n = 680, ave = 1139831.82, largest = 25259389
N50 = 9253497, n = 27
N60 = 7586889, n = 36
N70 = 6427024, n = 47
N80 = 5073214, n = 61
N90 = 2269709, n = 81
N100 = 517, n = 680
N_count = 0
Gaps = 0
-------------------------------------------------------------------------------
stats for flye.racon3.fasta
sum = 775642180, n = 680, ave = 1140650.26, largest = 25269729
N50 = 9260756, n = 27
N60 = 7588618, n = 36
N70 = 6429082, n = 47
N80 = 5078123, n = 61
N90 = 2270025, n = 81
N100 = 517, n = 680
N_count = 0
Gaps = 0

pilon跑的非常慢！不知道是不是哪里參數(shù)設(shè)置得不對(duì)，3輪pilon糾錯(cuò)跑了大概2.5天，結(jié)果如下：


pilon輸出結(jié)果


最終糾錯(cuò)結(jié)果分別在pilon1.out，pilon2.out，pilon3.out中，每個(gè)目錄里包含了糾錯(cuò)后的基因組fasta文件及糾錯(cuò)過程changes文件。
assembly-stats查看，基因組大小依次變小，準(zhǔn)確性應(yīng)該是上升了，等待后續(xù)組裝質(zhì)量評(píng)估。

$assembly-stats pilon1.out/racon3.pilon1.fasta pilon2.out/racon3.pilon2.fasta pilon3.out/racon3.pilon3.fasta
stats for pilon1.out/racon3.pilon1.fasta
sum = 774270266, n = 680, ave = 1138632.74, largest = 25232310
N50 = 9243114, n = 27
N60 = 7577994, n = 36
N70 = 6420177, n = 47
N80 = 5068385, n = 61
N90 = 2266107, n = 81
N100 = 517, n = 680
N_count = 0
Gaps = 0
-------------------------------------------------------------------------------
stats for pilon2.out/racon3.pilon2.fasta
sum = 774210928, n = 680, ave = 1138545.48, largest = 25230882
N50 = 9242790, n = 27
N60 = 7578144, n = 36
N70 = 6419986, n = 47
N80 = 5068374, n = 61
N90 = 2266189, n = 81
N100 = 517, n = 680
N_count = 0
Gaps = 0
-------------------------------------------------------------------------------
stats for pilon3.out/racon3.pilon3.fasta
sum = 774169945, n = 680, ave = 1138485.21, largest = 25229904
N50 = 9242646, n = 27
N60 = 7577699, n = 36
N70 = 6420023, n = 47
N80 = 5068254, n = 61
N90 = 2266164, n = 81
N100 = 517, n = 680
N_count = 0
Gaps = 0

polish策略2

直接使用nextpolish軟件進(jìn)行三代和二代數(shù)據(jù)同時(shí)糾錯(cuò)。
1# 軟件安裝
nextpolish安裝
軟件主頁：https://github.com/Nextomics/NextPolish

## 手動(dòng)安裝
pip install paralleltask
wget https://github.com/Nextomics/NextPolish/releases/download/v1.4.1/NextPolish.tgz
tar -vxzf NextPolish.tgz
cd NextPolish
make

## 自動(dòng)安裝
conda install -c bioconda nextpolish

2# 運(yùn)行程序
nextPolish的使用跟前面兩個(gè)軟件有些不同，是要先填寫一個(gè).cfg配置文件，然后根據(jù)這個(gè)配置文件軟件來決定怎么開展polish。軟件官網(wǎng)提供了4種配置文件案例，分別為只用二代短reads糾、只用三代長reads糾、同時(shí)用短reads和長reads糾以及用短reads和hifi reads糾。我們這次選擇第三種配置方式。

  • 首先準(zhǔn)備短reads路徑文件
ls ../../../../XX_survey/fastp.out/liver.clean.CLEAN.R1.fastq.gz  ../../../../XX_survey/fastp.out/liver.clean.CLEAN.R2.fastq.gz > nextpolish.sgs.fofn
  • 同時(shí)準(zhǔn)備長reads路徑文件
ls ../../P01TYD20308306-1_r64030_20201110_065049_1_A02.subreads.fasta > nextpolish.lgs.fofn
  • 其次配置.cfg文件
$ vi nextpolish.run.cfg
------------------------------------------------------------------------------------------------
[General]
job_type = local     # 運(yùn)行方式本地還是集群（local, sge, pbs）
job_prefix = nextPolish
task = best     # 運(yùn)行策略，2輪三代糾錯(cuò)+4輪二代糾錯(cuò)，如果后面只有三代或者二代數(shù)據(jù)，則只運(yùn)行三代或者二代糾錯(cuò)
rewrite = yes
rerun = 3
parallel_jobs = 6    # 平行任務(wù)數(shù)
multithread_jobs = 4    # 每個(gè)任務(wù)所給的線程數(shù)
genome = ../../purge_dups/flye_purge/assembly.fasta    # 初始基因組
genome_size = auto    # 自動(dòng)計(jì)算上面給的基因組大小
workdir = ./nextpolish_rundir    # 建工作路徑
polish_options = -p {multithread_jobs}     # 按照上面配置的平行任務(wù)的策略跑

[sgs_option]
sgs_fofn = ./nextpolish.sgs.fofn    # 前面生成的二代測(cè)序數(shù)據(jù)路徑文件
sgs_options = -max_depth 100 -bwa     # bwa比對(duì)，用最大100*的數(shù)據(jù)去做糾錯(cuò)，也可以用全部的

[lgs_option]
lgs_fofn = ./nextpolish.lgs.fofn     # 前面生成的二三代測(cè)序數(shù)據(jù)路徑文件
lgs_options = -min_read_len 1k -max_depth 100     #minimap2比對(duì)，同樣用最大100*的數(shù)據(jù)去做糾錯(cuò)
lgs_minimap2_options = -x map-pb      #minimap2比對(duì)方式，pacbio數(shù)據(jù)比對(duì)到基因組上。
  • 最后nextPolish運(yùn)行.cfg即可
/newlustre/home/jfgui/Wangtao/software/NextPolish/nextPolish  ./nextpolish.run.cfg

3# 結(jié)果查看

nextpolish輸出結(jié)果

最后結(jié)果在nextpolish_rundir目錄里，genome.nextpolish.fasta為最終糾錯(cuò)結(jié)果，其中的00和01文件分別記錄了2輪三代糾錯(cuò)結(jié)果，02、03、04、05文件夾則分別記錄了2種糾錯(cuò)方式的4輪二代糾錯(cuò)。每個(gè)文件夾都有一個(gè)input.genome.fasta表示糾錯(cuò)前的基因組，也就前一輪糾錯(cuò)的結(jié)果。

最后我們查看糾錯(cuò)前后結(jié)果對(duì)比，基因組大小糾錯(cuò)后少了0.5M，準(zhǔn)確性肯定提高了只不過這個(gè)方法看不出來??

$assembly-stats  nextpolish_rundir/00.lgs_polish/input.genome.fasta  nextpolish_rundir/genome.nextpolish.fasta
stats for nextpolish_rundir/00.lgs_polish/input.genome.fasta # 糾錯(cuò)前
sum = 771989320, n = 680, ave = 1135278.41, largest = 25186124
N50 = 9211426, n = 27
N60 = 7566198, n = 36
N70 = 6405554, n = 47
N80 = 5050486, n = 61
N90 = 2266136, n = 81
N100 = 493, n = 680
N_count = 0
Gaps = 0
-------------------------------------------------------------------------------
stats for nextpolish_rundir/genome.nextpolish.fasta  # 糾錯(cuò)后
sum = 771468932, n = 680, ave = 1134513.14, largest = 25167817 
N50 = 9204720, n = 27
N60 = 7567313, n = 36
N70 = 6408529, n = 47
N80 = 5048648, n = 61
N90 = 2265774, n = 81
N100 = 493, n = 680
N_count = 0
Gaps = 0

其他相關(guān)好文推薦：racon+pilon進(jìn)行三代+二代數(shù)據(jù)糾錯(cuò) - 簡書 (jianshu.com)

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡書系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容