写在前面
在进行染色体挂载之前一般对得到的 "primary assembly"（主要组装/初始组装结果）进行进一步的优化，以减少错误，提高基因组的质量和准确性。然后再将这些相对碎片化但polish过的高质量primary assembly进行染色体挂载以得到一个染色体水平的组装结果。如果用CLR数据组装出的primary assembly在没有进行polish的情况下，直接用HIC数据进行染色体挂载会有如下问题：（以下几点回答来自gpt3.5）
错误和噪声： CLR（Continuous Long Reads，如 PacBio 或 Nanopore 数据）通常具有较高的错误率，可能会导致组装的主要序列中存在错误和噪声。直接使用这些序列进行 Hi-C 染色体挂载可能会将这些错误传播到联系矩阵中，影响到染色体的准确性。
联系矩阵质量下降： 染色体挂载的关键是基于 Hi-C 数据生成染色体之间的联系矩阵。如果主要序列存在错误，这些错误可能会影响到联系矩阵的质量，从而降低染色体挂载的准确性。
组装断裂： 主要序列中的错误可能会导致组装断裂，使得 Hi-C 数据在某些区域无法正确匹配。这会导致染色体挂载失败或出现错误。
重复序列挑战： 高复杂性和重复性区域在 de novo 组装中是一个挑战，而 CLR 数据由于错误率较高，可能在这些区域中表现不佳。这会导致在 Hi-C 染色体挂载中难以准确匹配和定位这些区域。
因此，为了获得更准确的染色体挂载结果，建议在进行 Hi-C 染色体挂载之前，对主要序列进行纠错（polishing），以减少错误和噪声。纠错可以使用高质量的测序数据（如 illumina 数据）对主要序列进行校正，提高组装的质量和准确性。这样，在染色体挂载过程中，你会使用更可靠的主要序列和联系矩阵，从而获得更准确的染色体边界和结构信息。
另外，我们可以同时用三代和二代数据一起对我们的初始组装进行纠错。三代数据纠错主要用于解决长重复区域和复杂结构，提高基因组的连续性和准确性，有助于捕获染色体级别的信息；而二代数据纠错主要用于在主要序列中纠正小错误，提高基因组的局部准确性，从而使基因组序列精准度更接近于真实的 DNA 序列。综合而言，使用三代+二代数据进行纠错是一种综合优化的策略，有助于克服不同测序技术的优缺点，提高基因组的整体准确性和质量。这两种数据的结合能够产生更好的组装和纠错效果，为后续的基因组研究提供更可靠的基因组参考。
这里我们我们用了两种策略：
- 策略1：racon+pilon，即用racon软件进行三代数据纠错，再用pilon软件进行二代数据纠错；
- 策略2：直接使用nextpolish软件进行三代和二代数据同时纠错。
数据准备
1.primary assembly初始组装文件，这里用flye组装的初始组装做演示。
1.用于的组装的pacbio CLR三代测序文件。
2.经过fastp质控过的illumina二代测序文件。
polish策略1
racon+pilon，即用racon软件进行三代数据纠错，再用pilon软件进行二代数据纠错。
1# 软件安装
racon安装
软件主页：https://github.com/isovic/racon
## 自动安装
conda install -c bioconda racon
## 手动安装
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
## 由于是个java程序直接下载即可使用
wget https://github.com/broadinstitute/pilon/releases/download/v1.24/pilon-1.24.jar
## java程序统一使用方法：java -jar 后面接java程序
java -jar pilon-1.24.jar --help # 查看pilon用法
其他软件依赖：因为涉及到比对，所以需要安装三代和二代比对软件minimap2和bwa
minimap2主页：https://github.com/lh3/minimap2
bwa主页：https://github.com/lh3/bwa
## minimap2自动安装
conda install -c bioconda minimap2
## minimap2手动安装
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手动安装
git clone https://github.com/lh3/bwa.git
cd bwa
make
2# 运行程序
我这里我打算先用racon进行3轮三代数据纠错，再用pilon对racon的最后一轮结果再进行3轮二代数据纠错，加起来相当于6轮了。
- 用racon进行3轮三代数据纠错
## 第1轮纠错，这里详细注释了，后两轮就不展开了
/newlustre/home/jfgui/wy/anaconda3/envs/pengfang/bin/minimap2 \ # 用minimap2比对
-ax map-pb \ # 用-a参数生成sam文件而非默认的paf格式，-x参数为比对模式，这里是pacbio数据比对到参考基因组
-t 24 \ # 线程数
../../purge_dups/flye_purge/assembly.fasta \ # 参考基因组，这里第1轮就是用的flye组装出的初始组装，后面就都是前一轮的纠错结果。
../../P01TYD20308306-1_r64030_20201110_065049_1_A02.subreads.fasta \ # 用于纠错三代测序文件
| gzip -c - > minimap1.sam.gz # 压缩后输出为 minimap1.sam.gz
/newlustre/home/jfgui/wy/anaconda3/envs/pengfang/bin/racon \ # 用racon纠错
-t 24 \ # 线程数
-u \ # 纠错过和未被就错的contig都会被输出，否则只输出纠错过的contig
../../P01TYD20308306-1_r64030_20201110_065049_1_A02.subreads.fasta \ # 用于纠错三代测序文件
minimap1.sam.gz \ # 前面的比对结果
../../purge_dups/flye_purge/assembly.fasta > flye.racon1.fasta # 纠错前的基因组指向纠错后的基因组
## 第2轮纠错
/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轮纠错
/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进行3轮二代数据纠错
## 第1轮纠错，这里详细注释了，同理后两轮就不展开了
ln -s ../flye.racon3.fasta ./racon3.fasta # 将需要纠错的基因组链接到到当前工作路径
/newlustre/home/jfgui/Wangtao/software/bwa/bwa index racon3.fasta # 对基因组构建index以便后面的比对
/newlustre/home/jfgui/Wangtao/software/bwa/bwa mem \ # 用bwa的mem比对模式对二代数据进行比对
-t 24 \ #线程数
racon3.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 \ # 用samtools的sort命令对比对结果进行排序，用管道符连接
-@ 24 \ # samtools的线程数
- \ # 表示管道符前面命令的输出作为管道符后面命令的输入，这里也就是指未排序的bam文件
-o bwamem1.bam # 输出比对结果
/newlustre/home/jfgui/wy/anaconda3/envs/pengfang/bin/samtools index bwamem1.bam # 用samtools对比对文件进行index
java -Xmx128G -jar /newlustre/home/jfgui/Wangtao/software/pilon-1.24.jar \ # 用pilon纠错，-Xmx128G指的是用128G内存
--genome racon3.fasta \ # 前一轮纠错的基因组
--frags bwamem1.bam \ # 比对结果文件
--changes \ # 纠错前后变化记录文件
--diploid \ # 二倍体，不设置会默认单倍体进而识别单倍型
--threads 24 \ # 线程数
--outdir ./pilon1.out \ # 输出目录
--output racon3.pilon1 # 输出前缀
## 第2轮纠错
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轮纠错
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# 结果查看
3轮racon结果
3轮racon纠错后其结果还是比较简洁的，就3个比对文件3个纠错后基因组文件。
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跑的非常慢！不知道是不是哪里参数设置得不对，3轮pilon纠错跑了大概2.5天，结果如下：
pilon输出结果
最终纠错结果分别在pilon1.out，pilon2.out，pilon3.out中，每个目录里包含了纠错后的基因组fasta文件及纠错过程changes文件。
assembly-stats查看，基因组大小依次变小，准确性应该是上升了，等待后续组装质量评估。
$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软件进行三代和二代数据同时纠错。
1# 软件安装
nextpolish安装
软件主页：https://github.com/Nextomics/NextPolish
## 手动安装
pip install paralleltask
wget https://github.com/Nextomics/NextPolish/releases/download/v1.4.1/NextPolish.tgz
tar -vxzf NextPolish.tgz
cd NextPolish
make
## 自动安装
conda install -c bioconda nextpolish
2# 运行程序
nextPolish的使用跟前面两个软件有些不同，是要先填写一个.cfg配置文件，然后根据这个配置文件软件来决定怎么开展polish。软件官网提供了4种配置文件案例，分别为只用二代短reads纠、只用三代长reads纠、同时用短reads和长reads纠以及用短reads和hifi reads纠。我们这次选择第三种配置方式。
- 首先准备短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
- 同时准备长reads路径文件
ls ../../P01TYD20308306-1_r64030_20201110_065049_1_A02.subreads.fasta > nextpolish.lgs.fofn
- 其次配置.cfg文件
$ vi nextpolish.run.cfg
------------------------------------------------------------------------------------------------
[General]
job_type = local # 运行方式本地还是集群（local, sge, pbs）
job_prefix = nextPolish
task = best # 运行策略，2轮三代纠错+4轮二代纠错，如果后面只有三代或者二代数据，则只运行三代或者二代纠错
rewrite = yes
rerun = 3
parallel_jobs = 6 # 平行任务数
multithread_jobs = 4 # 每个任务所给的线程数
genome = ../../purge_dups/flye_purge/assembly.fasta # 初始基因组
genome_size = auto # 自动计算上面给的基因组大小
workdir = ./nextpolish_rundir # 建工作路径
polish_options = -p {multithread_jobs} # 按照上面配置的平行任务的策略跑
[sgs_option]
sgs_fofn = ./nextpolish.sgs.fofn # 前面生成的二代测序数据路径文件
sgs_options = -max_depth 100 -bwa # bwa比对，用最大100*的数据去做纠错，也可以用全部的
[lgs_option]
lgs_fofn = ./nextpolish.lgs.fofn # 前面生成的二三代测序数据路径文件
lgs_options = -min_read_len 1k -max_depth 100 #minimap2比对，同样用最大100*的数据去做纠错
lgs_minimap2_options = -x map-pb #minimap2比对方式，pacbio数据比对到基因组上。
- 最后nextPolish运行.cfg即可
/newlustre/home/jfgui/Wangtao/software/NextPolish/nextPolish ./nextpolish.run.cfg
3# 结果查看
nextpolish输出结果
最后结果在nextpolish_rundir目录里，genome.nextpolish.fasta为最终纠错结果，其中的00和01文件分别记录了2轮三代纠错结果，02、03、04、05文件夹则分别记录了2种纠错方式的4轮二代纠错。每个文件夹都有一个input.genome.fasta表示纠错前的基因组，也就前一轮纠错的结果。
最后我们查看纠错前后结果对比，基因组大小纠错后少了0.5M，准确性肯定提高了只不过这个方法看不出来😂
$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 # 纠错前
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 # 纠错后
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
