コンティグからラフにテロメア配列を探したい

背景

HiFi リードをhifiasmでアセンブルすると、運がいいとT2Tで繋がったコンティグも出力されるようになった。

目的

実際にどのコンティグがT2Tなのか調べたい。

方法

テロメア配列が密集している領域をコンティグから探してくる。
###環境
samtools, bedtools がインストールされていればOK.
seqkitはお好みで。

Step.0 コンティグを長さ順にソート(やってもやらなくてもよい)

ついでに100bpごとにfastaを区切っている

seqkit sort -l -r hoge.p_ctg.fa | \
seqkit seq -w 100 \
>sorted_hoge.bp.p_ctg.fa

Step.1 samtoolsでindex作ってchrom.sizes作る

samtools faidx sorted_hoge.bp.p_ctg.fa
cut -f 1,2 sorted_hoge.bp.p_ctg.fa.fai \
>sorted_hoge.bp.p_ctg.chrom.sizes

Step.2　bedtools makewindowsで任意のサイズのbinに区切ったbedファイルを作る

今回は10kbごとに区切った

bedtools makewindows -g sorted_hoge.bp.p_ctg.chrom.sizes -w 10000 >sorted_hoge.bp.p_ctg_10kb.bed

Step.3　bedtools getfastaで任意のサイズのbinに区切ったbedファイルに配列情報を付け足す

bedtools getfasta -fi sorted_hoge.bp.p_ctg.fa \
-bed sorted_hoge.bp.p_ctg_10kb.bed -tab \
>sorted_hoge.bp.p_ctg_10kb_seq.bed

Step.4　シェルでテロメア配列が任意のbinの中に何個存在するか数え上げてwigファイルとして出力

cat sorted_hoge.bp.p_ctg_10kb_seq.bed | \
awk 'BEGIN {OFS="\t"}{print $1,$2,$3, gsub(/AAACCCT|AGGGTTT/,"",$4)}' \
> sorted_hoge.bp.p_ctg_10kb_seq_telomere.wig

結果

IGVで眺めてみた

コンティグの両末端にテロメア配列っぽいのピークがみえるのでT2Tっぽいが、中央部にもテロメア配列がかたまっていたりするので、よくわからん。

考察

なんかもっとスマートな方法がある気がするので誰か教えて！