[bioinfo]nanopore readの精度向上のためのbasecalling法 (250129 ver.)
- Obtain pod5 data by Minknow(Minknowでpod5データを取得する)
- Basecall using latest version of dorado and the model with –emit-moves(doradoで–emit-movesを付けてベースコールする)
- Identify duplex reads by duplex_tools(duplex_toolsでduplex readを判別する)
- Basecall duplex reads by using dorado(doradoでduplex basecallを行う)
- Split fused template and complement reads by duplex_tools’s split-pair(連続してポアに入ったduplex readを分割する)
- Basecall duplex reads by using dorado(splitしたreadについてdoradoでduplex basecallを行う)
- Convert unmapped bam to fastq of simplex reads(unmapped bamからfastqに変換)
- Purge duplex reads and low quality reads from simplex fastq data and concatenate(simplexデータからクオリティが低いデータとduplex readを排除して結合する)
1. Obtain pod5 data by Minknow(Minknowでpod5データを取得する)
普通にシーケンスランを行う、ベースコールはしなくても良いがデータの概要はした方が正確にわかる。
2. Basecall using latest version of dorado and the model with –emit-moves(doradoで–emit-movesを付けてベースコールする)
duplex_toolsの解説ではhacでベースコールするとあるが、結局最新版のモデルでベースコールした方がsimplexリードの正確性が高いため、GPUの余裕があればsupでベースコールする”sup”と書くだけでpod5データから適したモデルを自動的にダウンロードしてベースコールしてくれる。
dorado basecaller sup -r ./ --emit-moves > unmapped _reads_with_moves.bam
3. Identify duplex reads by duplex_tools(duplex_toolsでduplex readを判別する)
duplex_toolsはpython 3.10下以外ではエラーが生じる為、condaでpython 3.10環境を作った上で、pipでインストールする。
#create python 3.10 env
conda create -n python310 python=3.10
conda activate python310
#install duplex_tools
python -m venv venv --prompt duplex
. venv/bin/activate
pip install duplex_tools
#judge duplex reads using duplex_tools
duplex_tools pair --output_dir pairs_from_bam unmapped_reads_with_moves.bam
4. Basecall duplex reads by using dorado(doradoでduplex basecallを行う)
#basecall all accessible pod5 data under sequence data folder
dorado duplex sup -r ./ --pairs pairs_from_bam/pair_ids_filtered.txt > duplex_orig.sam
5. Split fused template and complement reads by duplex_tools’s split-pair(連続してポアに入ったduplex readを分割する)
duplex_tools split_pairs unmapped_reads_with_moves.bam ./ pod5_splitduplex/
#concat extracted split ids
cat pod5_splitduplex/*_pair_ids.txt > split_duplex_pair_ids.txt
6. Basecall duplex reads by using dorado(splitしたreadについてdoradoでduplex basecallを行う)
今度はsplitしたreadのベースコールを行う。
dorado duplex sup pod_splitduplex/ --pairs split_duplex_pair_ids.txt > duplex_splitduplex.sam
7. Convert unmapped bam to fastq of simplex reads(unmapped bamからfastqに変換)
gatkを使ってsamあるいはbamからfastqに変換する。
gatk SamToFastq -I unmapped_reads_with_moves.bam -FASTQ unmapped_reads_with_moves.fastq
gatk SamToFastq -I duplex_orig.sam -FASTQ duplex_orig.fastq
gatk SamToFastq -I duplex_splitduplex.sam -FASTQ duplex_splitduplex.fastq
8. Purge duplex reads and low quality reads from simplex fastq data and concatenate(simplexデータからクオリティが低いデータとduplex readを排除して結合する)
NanoFiltとseqkit grepを使ってクオリティが低いリードと重複したリードを排除する。クオリティが低いリードが入るのはすべてのpod5も入れてベースコールしているため。
#convert space of list file to change row for seqkit
sed 's/ /\n/g' pairs_from_bam/pair_ids_filtered.txt > pairs_from_bam/pair_ids_filtered+.txt
sed 's/ /\n/g' split_duplex_pair_ids.txt > split_duplex_pair_ids+.txt
#concatenate duplex read list
cat pairs_from_bam/pair_ids_filtered+.txt split_duplex_pair_ids+.txt > duplex_list.txt
#purge duplex reads from simplex basecall data
seqkit grep -f duplex_list.txt -v unmapped_reads_with_moves.fastq -o unmapped_reads_with_moves_only_simplex.fastq
#concatenate after quality filtration
cat unmapped_reads_with_moves_only_simplex.fastq duplex_orig.fastq duplex_splitduplex.fastq | NanoFilt -q 10 > simplex+duplex.fastq