小马的生信笔记

引言

记录一次使用quarTeT+TGS-GapCloser+Gapsuite补基因组gap的过程，纯记录文章，只有流程的描述，不会有任何原理的讲解。

安装软件

				
					# quarTeT
conda create -n gapfill --channel conda-forge --channel bioconda python=3.11.4 minimap2=2.26 mummer4=4.0.0rc1 trf=4.09.1 cd-hit=4.8.1 blast=2.14.0 tidk=0.2.31 r=4.3 r-rideogram=0.2.2 r-ggplot2=3.4.4 gnuplot=5.4 unimap=0.1
git clone https://github.com/aaranyue/quarTeT
# mm2plus
mamba install bioconda::mm2plus
# TGS-GapCloser
mamba install bioconda::tgsgapcloser
# Gapsuite
mamba install bioconda::jellyfish
mamba install bioconda::seqkit
git clone https://github.com/panlab-bioinfo/RAfilter.git
cd RAfilter/src
make

填补Gap

hifi数据的bam文件转fa

				
					pbindex ccs.bam
bam2fasta -o out ccs.bam

ont的测序数据fq转为fa

				
					seqkit fq2fa -j 50 pass.raw.fq > pass.raw.fa

quarTeT初步填补gap

				
					python ~/Software/quarTeT/quartet.py GapFiller \
-d all_ont_fill.genome.filled.fasta \
-g out.fasta hifiasm.hic.haps.p_ctg.gfa.fa pass.raw.fa \
-t 60 -p all_ont_hifi_fill

TGS-GapCloser再次填补

				
					tgsgapcloser --scaff all_ont_hifi_fill.genome.filled.chr10a.fasta \
--reads pass.raw.fa \
--output tgsgapcloser_ont --ne
				
			
				
					tgsgapcloser --scaff all_ont_hifi_fill.genome.filled.chr10a.fasta \
--reads  out.fasta \
--output tgsgapcloser_ont --ne

这一步结束，只有1个gap没有被补上，随后使用Gapsuite手动填补。Gapsuite产生的文件很大，因此我们只提取存在gap的染色体跑Gapsuite。

				
					seqkit grep -p  "chr10a" all_ont_hifi_fill.genome.filled.fasta > all_ont_hifi_fill.genome.filled.chr10a.fasta

Gapsuite手动补gap

				
					bash pipeline.sh -r ont -t 60 all_ont_hifi_fill.genome.filled.chr10a.fasta pass.raw.fa

发表评论

基因组专栏