筆者がMacのターミナルでよく使用する，bioinformatics関連のツールについて，使い方を簡単にまとめる。
細かいoption等についてはほとんど紹介しないので，各ツールの公式サイトをそれぞれ参照してほしい。
(2023/12/30更新, ver.2)
seqkit
配列の取得
seqkit grep -nrp [目的のcontig名など] reference.fasta > output.fasta
相補鎖にする
seqkit seq -pr before.fasta > output.fasta
特定の染色体の，特定の位置から配列を取得
seqkit faidx reference.fasta chr1:1-10000 -o output.fasta
各染色体 or contigごとの長さを確認
seqkit fx2tab -nl data.fasta
seqkit subseq --chr Chr5 -r 38:311 omes.fasta > a00.fasta
seqkit grep -f list hoge.fasta > hoge.tig0103.fasta
local BLAST
databaseの作成
-dbtypeは塩基配列ならnucl，アミノ酸配列ならprotを指定。
makeblastdb -in reference.fasta -parse_seqids -dbtype nucl -out database_name
blastn
-outfmtは指定しなくても良い。
blastn -query query.fasta -db database_name -evalue 1e-30 -outfmt 6 -out result.txt
blastp
nucmer
nucmerで作るdot plot
nucmer --prefix=data reference.fa query.fa
delta-filter -q -r data.delta > filtered.delta
mummerplot --postscript -p mapping filtered.delta
gnuplot mapping.gp
今は.psファイルをmacのプレビューで開くことができない。Inkscapeなどのソフトで開き，pdfなどに変換して保存する必要がある。
circos
circosの前処理の部分をまとめる。
nucmer --maxmatch -c 1000 --prefix=output Reference.fasta query.fasta
delta-filter -l 4000 output.delta > output_filtered.delta
show-coords -r -c -l output_filtered.delta > output.txt
ここでoutput.txtをexcelで開き，必要なところだけを取り出し，整理する。
具体的には[S1]，[E1]，[S2]，[E2]，[TAGS]が必要となる。
以下のような順番に並び替える。
| (S1, E1の)TAGS | [S1] | [E1] | (S2, E2の)TAGS | [S2] | [E2] | color=chr_num |
|---|
tab形式ファイル(output_sort.txt)で保存し，以下のスクリプトを実行する。
cat output_sort.txt | gsed -z 's/\r/\n/g' | gsed -z 's/\t/ /g' > output.bundle.txt
output.bundle.txtをcircos plotを描画するのに用いる。
この他にkaryotypeファイルなどを用意しておき，confファイルを実行する。
minimap2
longread (Hifi)を用いたmapping
minimap2 -t [threads_num] -ax map-hifi Assembly.fasta hifi_reads.fastq.gz > output.sam
このsamファイルを用いて，unmapped readを取り出すことも可能。(ここではAssembly.fastaにmappingされなかったreadを取り出すことにする)
samtools sort -@ [threads_num] -O bam -o output.bam output.sam
samtools index output.bam
bamtools split -in output.bam -mapped
#ここで出力されたoutput.UNMAPPED.bam (unmapped read)をfastqファイルに変換可能
samtools bam2fq output.UNMAPPED.bam > unmapped_ouput.fq
#このunmapped_ouput.fqを用いてassemblyなどに活用できる。
BUSCO
以下はBUSCO v5.4.4での動作を確認した(ver. によってコマンドが異なる?)。
-lでdataset名を指定。「busco --list-datasets」で利用可能なdatasetを事前に調べておく。
-mでmodeを指定。genome，transcriptome，proteinsの3つのmodeが存在するが，いずれも他のoptionsは同じように指定する。
busco -m genome -i sequence.fasta -o output -l fungi_odb10 -c [threads_num]
Telomereの検出
tapestryパッケージを用いてtelomereを検出する。
はじめに，以下のようにパッケージをダウンロードしておく。
conda install -c bioconda tapestry
以下のようにtelomereを検出する。-rでアセンブリに用いたreadを指定。-tで既知のテロメア配列を指定。
weave -a assembly.fasta -r read.fastq.gz -d 30 -t TTAGGGTCAACA -o output -c [threads_num] -f
mosdepthによるcoverage計算
筆者はlongreadをゲノムにmappingした際のcoverageを以下のように計算した。
# longreadをゲノムにmapping
minimap2 -t [threads_num] -ax map-hifi genome.fasta hifi_reads.fastq.gz | samtools sort -@ [threads_num] -O BAM -o output.bam - && samtools index output.bam
#mosdepthによるdepthの出力
mosdepth result output.bam -Q 0 -t [threads_num]
#reslut.mosdepth.summary.txt、が各染色体ごとのdepthが出力されたファイルとなる。
mafft/RAxML/FigTreeによる系統樹
アミノ酸配列を用いた系統樹
mafft --auto input.fasta > output.fasta
raxmlHPC -f a -m PROTCATMTMAMF -s output.fasta -# 100 -p 12345 -x 12345 -n HOGE -T [threads_num]
raxmlHPC -f b -m PROTCATMTMAMF -t RAxML_bestTree.HOGE -z RAxML_bootstrap.HOGE -n BOOTSTRAP
出力されたファイルのうち，RAxML_bipartitions.BOOTSTRAPをFigtreeで開き，系統樹の形などを調節し，PDF等で出力する。
塩基配列の場合もほとんど同じだが，raxmlHPCのところで，-mオプションの種類を変える必要がある。