Bioinformatics  command チートシート [備忘録]

筆者が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
```
Primerなど、短い配列をqueryにする場合は、-task オプションを追加する。
[公式](https://www.ncbi.nlm.nih.gov/books/NBK569839/)の説明ではblastn-shortを使うことで、『BLASTN program optimized for sequences shorter than 50 bases』とある。



```
blastn -query query.fasta -task blastn-short -db database_name -evalue 10 -outfmt 6 -out result.txt
```


### blastp 
blastnをblastpに変えるだけ。
```
blastp -query query.fasta -db database_name -evalue 1e-30 -outfmt 6 -out result.txt
```

## 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などに変換して保存する必要がある。

https://support.apple.com/ja-jp/108775

## 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パッケージ](https://github.com/johnomics/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等で出力する。

### 塩基配列を用いた系統樹
塩基配列の場合も、mafft等でアライメントしたファイルを用いて、以下のようにする。

```
mafft --auto input.fasta > output.fasta

raxml-ng --msa output.fasta --all --model GTR+G+I --bs-trees 100 --threads [threads_num]
```