为什么要Polish?
通过软件初步拼接后的基因组序列还需要执行polish(也就是校正)环节.
由于测序错误和拼接错误存在,测序中的reads在中的剪辑并非完全可靠.
可在此使用测序reads数据,将reads比对到组装后的基因组序列中,通过海量的reads识别及校正基因组序列中的错误位点.
这个过程俗称基因组Polish.polish后,基因组的可靠性得到了显著提升.
软件
根据用来polish的数据的不同可以分为,三代polish和二代polish:
Pilon: 主要针对使用二代测序数据对基因组序列进行polish
Racon: 主要针对使用三代测序数据对基因组序列进行polishPacBio、Nanopore
NextPolish: PacBio、Nanopore、Illumina都可以使用,即可使用单一的二代或三代测序数据,也可实现二+三代的混合polish
#一般进行多轮polish,例如:Pilon三轮+Racon三轮
一.pilon–对二代测序数据polish
mkdir pilon_1
cd pilon_1
1.链接测序组装后的基因组文件
ln -s 组装序列文件 contig.fa #为了方便流程化运行,将链接文件命名为contig.fa或scaffolds.fa;
#组装序列文件就是需要抛光polish的基因组文件
2.利用bwa和samtools软件比对并排序基因组文件
bwa适用于短reads序列比对
bwa index contig.fa #创建基因组文件的index
索引用于bwa比对
bwa mem -t 12 contig.fa $fq1 $fq2 | samtools sort - -o reads.bam #$fq1和$fq2是illumina二代测序得到的短reads,将其比对到组装基因组上并排序.
samtools index reads.bam #创建对比结果的index文件
3.pilon纠错
java -Xmx20G -jar
pilion.jar路径–genome ./contig.fa –frags reads.bam –fix all –changes –vcf –diploid –threads 12 –outdir ./ –output contig_pilon
sed -i ‘s/_//‘ contig_pilon
常用选项参数说明：
输入参数
–genome 提供输入参考基因组
–frags 表示输入 < 1kb 的文库BAM
–jumps 输入 > 1kb 的文库BAM
–unpaired 输入非配对的BAM。
输出参数
–output表示输出的前缀
–outdir表示输出文件夹
–changes 列出发生变化的部分,以FASTA形式保存
–vcf 以VCF形式保存。
运行参数
–fix 声明对参考基因组做哪方面的改进，包括“snps”,”indels”,”gaps”,”local”, 默认是”all”
–threads : 设置线程数
可能是因为java写的程序的缘故，pilon的运行速度比较慢。最终会生成*
contig_pillon.changes：展示哪些位点经过了纠错；
contig_pillon.fasta ：最终纠错得到的结果；
contig_pillon.vcf ：vcf格式展示纠错的位点。
二代polish三代组装的基因组，反复几轮循环后，changes文件一般都能到0。
二代polish二代组装的基因组，除非特别小的基因组，否则可能永远达不到0。这种情况下，
可以先来几个循环，从已得到的所有changes文件中找一个最小的，其对应的fasta文件可作为最终结果。
第二轮polish
mkdir pilon_2
cd pilon_2
bwa index ../pilon_1/contig_pilon.fasta
bwa mem -t 12 ../pilon_1/contig_pilon.fasta $fq1 $fq2 | samtools sort - -o reads.bam
samtools index reads.bam
java -Xmx20G -jarpilion.jar路径–genome ../pilon_1/contig_pilon.fasta –frags reads.bam –fix all –changes –vcf –diploid –threads 12 –outdir ./ –output contig_pilon
第三轮polish
mkdir polish_3
cd polish_3
bwa index ../polish_2/contig_pilon.fasta
bwa mem -t 12 ../polish_2/contig_pilon.fasta $fq1 $fq2 |samtools sort - -o reads.bam
samtools index reads.bam
java -Xmx20G -jarpilion.jar路径–genome ../pilon_2/contig_pilon.fasta –frags reads.bam –fix all –changes –vcf –diploid –threads 12 –outdir ./ –output contig_pilon
二.Racon–对三代测序数据polish
Racon 可以使用illumina、pacbio和nanopore数据进行polish的软件,一般我们使用三代数据作为输入.
cd $workdir
mkdir polish_racon
cd polish_racon
ln -s ../polish_pilon/pilon_3/contig_pilon.fasta contig.fa #将contig.fa链接需要polish的序列文件
第一轮纠错polish
mkdir racon_1
cd racon_1
minimap2 -a -x map-ont -t 12 ../contig.fa $ont | gzip - >racon_1.sam.gz #将测序的nanopore数据比对到基因组序列上,并重定向输出用gzip压缩成.sam.gz文件
racon -t 12 $ont racon_1.sam.gz ../contig.fa > racon_1.fasta #根据比对结果,利用nanopore数据对基因组序列进行纠错polish
第二轮纠错(重复第一轮操作,将原基因组序列文件替换为第一轮纠错后的基因组序列文件)
mkdir racon_2
cd racon_2
minimap2 -a -x map-ont -t 12 ../racon_1/racon_1.fasta $ont | gzip - > racon_2.sam.gz
racon -t 12 $ont racon_2.sam.gz racon_1.fasta > racon_2.fasta
第三轮纠错
mkdir racon_3
cd racon_3
minimap2 -a -x map-ont -t 12 ../racon_2/racon_2.fasta $ont | gzip - > racon_3.sam.gz
racon -t 12 $ont racon_3.sam.gz ../racon_2/racon_2.fasta > racon_3.fasta
三.nextPolish
1.使用nextPolish需要准备
一个基因组文件contig.fa
ln -s 需要Polish的基因组文件路径 contig.fa
一个配置文件run.cfg
#参数说明：https://nextpolish.readthedocs.io/en/latest/OPTION.html
echo ‘
[General]
job_type = local ## 运行环境，可选local、sge、pbs(可以使用任务系统,但是我目前还不会-.-)；
job_prefix = nextPolish #任务名称，随便起个；
task = best ## task need to run [all, default, best, 1, 2, 5, 12, 1212…], 1, 2 are different algorithm modules for short reads, while 5 is the algorithm module for long reads, all=[5]1234, default=[5]12, best=[55]1212. (default: best)一个数字就代表进行一次相应模式的Polish
rewrite = yes
rerun = 3
parallel_jobs = 6 #并行运行的任务数；
multithread_jobs = 5 #每个任务占用的线程数；
genome = contig.fa ## 基因组文件
genome_size = auto #基因组大小，设置为auto就可以了；
workdir = ./output #指定输出工作路径；
polish_options = -p {multithread_jobs}
#短reads数据
[sgs_option]
sgs_fofn = ./sgs.fofn #需要将短reads数据文件路径输入到sga.fofn中
sgs_options = -max_depth 100 -bwa
长reads数据
[lgs_option]
lgs_fofn = ./lgs.fofn #需要将长reads数据文件路径输入到sga.fofn中
lgs_options = -min_read_len 1k -max_depth 100
lgs_minimap2_options = -x map-ont -t {multithread_jobs} ##其他数据修改minimap2 -x参 map-pb/map-hifi/map-ont/map-iclr - CLR/HiFi/Nanopore/ICLR vs reference mapping
#hifi数据
[hifi_option]
hifi_fofn = ./hifi.fofn
hifi_options = -min_read_len 1k -max_depth 100
hifi_minimap2_options = -x map-hifi‘ > run.cfg
一些数据文件sgs.fofn
echo “$datadir/illumina_1.fastq.gz” >./sgs.fofn
echo “$datadir/illumina_2.fastq.gz” >>./sgs.fofn
lgs.fofn
echo “$datadir/nanopore.fastq.gz” >lgs.fofn
hifi.fofn
echo “$datadir/pacbio_hifi.fastq.gz” >hifi.fofn
2.开始运行nextPolish
nextPolish run.cfg
