Samtoolsでタヌキを表示する

```
　　　　__,-─-､__
　　　（〆-─-ヽ)
　　　 （ ´・ω・｀ ）  < Samtools デ タヌキ ヲ ヒョウジ スル！
　　　/ 　,r‐‐‐､ヽ
　 　 し　ｌ　 x　）J
　　　_.'､ ヽ　 ノ.人
　(_((__,ノＵ´U. (酒)

```

## はじめに

> 役に立つバイオインフォマティックスのまじめな記事を期待してこのページを開いた皆さん。本当に申し訳ない。なにか有益な情報を期待してこのページを開いた方は静かに帰ってほしい。そういったものは一切ない。ふざけた内容で時間を失っても構わない人だけ読んでください。ちゃんと注意したからな。（タヌ記より）

この記事では、Samtoolsに `tanuki` というコマンドを追加して、タヌキのAAを表示することをやっていきます。



## モチベーション

そんなものはない。タヌキ表示したいって言ってんだろ！！

## タヌキのAAを探す

タヌキをターミナルに表示する以上、まずはタヌキのAAを探さなければならない。まずはこちらのサイトからタヌキのAAを入手する。

https://jihu.blog.fc2.com/blog-entry-487.html

```
　　　　__,-─-､__
　　　（〆-─-ヽ)
　　　 （ ´・ω・｀ ）
　　　/ 　,r‐‐‐､ヽ
　 　 し　ｌ　 x　）J
　　　_.'､ ヽ　 ノ.人
　(_((__,ノＵ´U. (酒)

```

↑タヌキ。

## Samtoolsのビルド

https://github.com/samtools/samtools

まずはSamtoolsが正常にビルドできることを確認しておこう。ビルドには htslib などの依存が必要になるので注意しよう。自分の場合は[趣味で作ってるCrystal言語のhtslibバインディング](https://github.com/bio-cr/hts.cr)のために、htslibを自分でビルドしたりするので、Samtoolsと同じ階層にhtslibディレクトリを置いている。ここに置くと、自動でsamtoolsはhtslibを検知してくれたりする。htslibのビルドは `git submodule update -i --recursive` などが必要。まあいろいろ引っかかると思うので詰まったらその都度ググってください（投げやり）。`git checkout 1.16.1` など安定版のタグに移動して `git switch tanuki` などとブランチを切っておくと安全だろう。そのあたりは好みで。

あとはGitHubのREADME.mdの指示に忠実にビルドすればよい。

```sh
autoheader            # Build config.h.in (this may generate a warning about
                      # AC_CONFIG_SUBDIRS - please ignore it).
autoconf -Wno-syntax  # Generate the configure script
./configure           # Needed for choosing optional functionality
make
```

## メイン関数はいずこ？

さて、最初に考えなければならないのは、`samtools` コマンドを入力したときのメイン関数はどこにあるのだろうか？ということである。
VSCodeでsamtoolsのリポジトリを開いて `int main(` などで検索。すると `bamtk.c` というファイルに `main` があることがわかる。

このファイルの途中に `samtools help` とやると表示されるヘルプっぽいものがある。

```c
static void usage(FILE *fp)
{
    /* Please improve the grouping */

    fprintf(fp,
"\n"
"Program: samtools (Tools for alignments in the SAM format)\n"
"Version: %s (using htslib %s)\n\n", samtools_version(), hts_version());
    fprintf(fp,
"Usage:   samtools <command> [options]\n"
"\n"
"Commands:\n"
"  -- Indexing\n"
"     dict           create a sequence dictionary file\n"
"     faidx          index/extract FASTA\n"
"     fqidx          index/extract FASTQ\n"
"     index          index alignment\n"
"\n"
"  -- Editing\n"
"     calmd          recalculate MD/NM tags and '=' bases\n"
"     fixmate        fix mate information\n"
"     reheader       replace BAM header\n"
"     targetcut      cut fosmid regions (for fosmid pool only)\n"
"     addreplacerg   adds or replaces RG tags\n"
"     markdup        mark duplicates\n"
"     ampliconclip   clip oligos from the end of reads\n"
"\n"
"  -- File operations\n"
"     collate        shuffle and group alignments by name\n"
"     cat            concatenate BAMs\n"
"     consensus      produce a consensus Pileup/FASTA/FASTQ\n"
"     merge          merge sorted alignments\n"
"     mpileup        multi-way pileup\n"
"     sort           sort alignment file\n"
"     split          splits a file by read group\n"
"     quickcheck     quickly check if SAM/BAM/CRAM file appears intact\n"
"     fastq          converts a BAM to a FASTQ\n"
"     fasta          converts a BAM to a FASTA\n"
"     import         Converts FASTA or FASTQ files to SAM/BAM/CRAM\n"
"     reference      Generates a reference from aligned data\n"
"\n"
"  -- Statistics\n"
"     bedcov         read depth per BED region\n"
"     coverage       alignment depth and percent coverage\n"
"     depth          compute the depth\n"
"     flagstat       simple stats\n"
"     idxstats       BAM index stats\n"
"     phase          phase heterozygotes\n"
"     stats          generate stats (former bamcheck)\n"
"     ampliconstats  generate amplicon specific stats\n"
"\n"
"  -- Viewing\n"
"     flags          explain BAM flags\n"
"     head           header viewer\n"
"     tview          text alignment viewer\n"
"     view           SAM<->BAM<->CRAM conversion\n"
"     depad          convert padded BAM to unpadded BAM\n"
"     samples        list the samples in a set of SAM/BAM/CRAM files\n"
"\n"
"  -- Misc\n"
"     help [cmd]     display this help message or help for [cmd]\n"
"     version        detailed version information\n"
"\n");
```

そこに `tanuki` コマンドのヘルプ追記する。あたりまえだが、tanuki は極めて重要なコマンドなので、ヘルプを参照して使い方を調べる人も多いかも知れぬ。ヘルプは重要である。

```diff
--- a/bamtk.c
+++ b/bamtk.c
@@ -203,6 +203,9 @@ static void usage(FILE *fp)
 "  -- Misc\n"
 "     help [cmd]     display this help message or help for [cmd]\n"
 "     version        detailed version information\n"
+"\n"
+"  -- Tanuki\n"
+"     tanuki         show a tanuki\n"
 "\n");
 }
```

それから、さらに下の方にスクロールしていくと、明らかに条件分岐をしているコードがある。


```c
    int ret = 0;
    if (strcmp(argv[1], "view") == 0)           ret = main_samview(argc-1, argv+1);
    else if (strcmp(argv[1], "import") == 0)    ret = main_import(argc-1, argv+1);
    else if (strcmp(argv[1], "mpileup") == 0)   ret = bam_mpileup(argc-1, argv+1);
    else if (strcmp(argv[1], "merge") == 0)     ret = bam_merge(argc-1, argv+1);
    else if (strcmp(argv[1], "sort") == 0)      ret = bam_sort(argc-1, argv+1);
    else if (strcmp(argv[1], "index") == 0)     ret = bam_index(argc-1, argv+1);
    else if (strcmp(argv[1], "idxstat") == 0 ||
             strcmp(argv[1], "idxstats") == 0)  ret = bam_idxstats(argc-1, argv+1);
    else if (strcmp(argv[1], "faidx") == 0)     ret = faidx_main(argc-1, argv+1);
    else if (strcmp(argv[1], "fqidx") == 0)     ret = fqidx_main(argc-1, argv+1);
    else if (strcmp(argv[1], "dict") == 0)      ret = dict_main(argc-1, argv+1);
    else if (strcmp(argv[1], "head") == 0)      ret = main_head(argc-1, argv+1);
    else if (strcmp(argv[1], "fixmate") == 0)   ret = bam_mating(argc-1, argv+1);
    else if (strcmp(argv[1], "rmdup") == 0)     ret = bam_rmdup(argc-1, argv+1);
    else if (strcmp(argv[1], "markdup") == 0)   ret = bam_markdup(argc-1, argv+1);
    else if (strcmp(argv[1], "ampliconclip") == 0) ret = amplicon_clip_main(argc-1, argv+1);
    else if (strcmp(argv[1], "flagstat") == 0 ||
             strcmp(argv[1], "flagstats") == 0) ret = bam_flagstat(argc-1, argv+1);
    else if (strcmp(argv[1], "calmd") == 0)     ret = bam_fillmd(argc-1, argv+1);
    else if (strcmp(argv[1], "fillmd") == 0)    ret = bam_fillmd(argc-1, argv+1);
    else if (strcmp(argv[1], "reheader") == 0)  ret = main_reheader(argc-1, argv+1);
    else if (strcmp(argv[1], "cat") == 0)       ret = main_cat(argc-1, argv+1);
    else if (strcmp(argv[1], "targetcut") == 0) ret = main_cut_target(argc-1, argv+1);
    else if (strcmp(argv[1], "phase") == 0)     ret = main_phase(argc-1, argv+1);
    else if (strcmp(argv[1], "depth") == 0)     ret = main_depth(argc-1, argv+1);
    else if (strcmp(argv[1], "coverage") == 0)  ret = main_coverage(argc-1, argv+1);
    else if (strcmp(argv[1], "bam2fq") == 0 ||
             strcmp(argv[1], "fastq") == 0 ||
             strcmp(argv[1], "fasta") == 0)     ret = main_bam2fq(argc-1, argv+1);
    else if (strcmp(argv[1], "pad2unpad") == 0) ret = main_pad2unpad(argc-1, argv+1);
    else if (strcmp(argv[1], "depad") == 0)     ret = main_pad2unpad(argc-1, argv+1);
    else if (strcmp(argv[1], "bedcov") == 0)    ret = main_bedcov(argc-1, argv+1);
    else if (strcmp(argv[1], "bamshuf") == 0)   ret = main_bamshuf(argc-1, argv+1);
    else if (strcmp(argv[1], "collate") == 0)   ret = main_bamshuf(argc-1, argv+1);
    else if (strcmp(argv[1], "stat") == 0 ||
             strcmp(argv[1], "stats") == 0)     ret = main_stats(argc-1, argv+1);
    else if (strcmp(argv[1], "flag") == 0 ||
             strcmp(argv[1], "flags") == 0)     ret = main_flags(argc-1, argv+1);
    else if (strcmp(argv[1], "split") == 0)     ret = main_split(argc-1, argv+1);
    else if (strcmp(argv[1], "quickcheck") == 0)  ret = main_quickcheck(argc-1, argv+1);
    else if (strcmp(argv[1], "addreplacerg") == 0) ret = main_addreplacerg(argc-1, argv+1);
    else if (strcmp(argv[1], "pileup") == 0) {
        fprintf(stderr, "[main] The `pileup' command has been removed. Please use `mpileup' instead.\n");
        return 1;
    }
    else if (strcmp(argv[1], "tview") == 0)   ret = bam_tview_main(argc-1, argv+1);
    else if (strcmp(argv[1], "ampliconstats") == 0)     ret = main_ampliconstats(argc-1, argv+1);
    else if (strcmp(argv[1], "samples") == 0)     ret = main_samples(argc-1, argv+1);
    else if (strcmp(argv[1], "consensus") == 0) ret = main_consensus(argc-1, argv+1);
    else if (strcmp(argv[1], "reference") == 0) ret = main_reference(argc-1, argv+1);
    else if (strcmp(argv[1], "version") == 0 || \
             strcmp(argv[1], "--version") == 0)
        long_version();
    else if (strcmp(argv[1], "--version-only") == 0) {
        printf("%s+htslib-%s\n", samtools_version(), hts_version());
    }
    else {
        fprintf(stderr, "[main] unrecognized command '%s'\n", argv[1]);
        return 1;
    }
    return ret;
}

```

ここに `main_tanuki` を追加する。これで `samtools tanuki` を実行すると `main_tanuki()` が呼び出されるようになる。

```diff
--- a/bamtk.c
+++ b/bamtk.c
@@ -287,6 +290,7 @@ int main(int argc, char *argv[])
     else if (strcmp(argv[1], "samples") == 0)     ret = main_samples(argc-1, argv+1);
     else if (strcmp(argv[1], "consensus") == 0) ret = main_consensus(argc-1, argv+1);
     else if (strcmp(argv[1], "reference") == 0) ret = main_reference(argc-1, argv+1);
+    else if (strcmp(argv[1], "tanuki") == 0) ret = main_tanuki();
     else if (strcmp(argv[1], "version") == 0 || \
              strcmp(argv[1], "--version") == 0)
         long_version();
```

上にスクロールして `main_tanuki()` 関数のプロトタイプ宣言を追加しておく。

```diff
--- a/bamtk.c
+++ b/bamtk.c
@@ -72,6 +72,7 @@ int main_import(int argc, char *argv[]);
 int main_samples(int argc, char *argv[]);
 int main_consensus(int argc, char *argv[]);
 int main_reference(int argc, char *argv[]);
+int main_tanuki(void);
```

これはなくてもビルドできるが、ないとコンパイル時に次のような警告メッセージが出るので、ちゃんと宣言しておく。

```c
bamtk.c: In function ‘main’:
bamtk.c:293:52: warning: implicit declaration of function ‘main_tanuki’ [-Wimplicit-function-declaration]
  293 |     else if (strcmp(argv[1], "tanuki") == 0) ret = main_tanuki(void);
      |                                                    ^~~~~~~~~~~
```

## main_tanuki() を実装する

`tanuki` コマンドはきわめて重要なので、多数の利用者が見込まれる。中にはタヌキをクジャクなどに改造しようとする者も現れるかもしれない。そのような将来を見据えて、きちんとヘッダーに著作権情報を明記しておく。

```c
/*  tanuki.c -- show you a tanuki.

    Copyright (C) 2022 Tanuki Institute.

    Author: kojix2

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
DEALINGS IN THE SOFTWARE.  */

#include <config.h>

#include <stdlib.h>
#include <stdio.h>

int main_tanuki(void)
{
    printf(
"　　　　__,-─-､__\n"
"　　　（〆-─-ヽ)\n"
"　　　 （ ´・ω・｀ ）\n"
"　　　/ 　,r‐‐‐､ヽ\n"
"　 　 し　ｌ　 x　）J\n"
"　　　_.'､ ヽ　 ノ.人\n"
"　(_((__,ノＵ´U. (酒)\n");

    return 0;
}

```

## Makefileを編集する

最後にMakefileを編集する。AOBJSに `tanuki.o` を追加する。`tanuki.o` が `tanuki.c` に依存する旨を記載する。

```diff
--- a/Makefile
+++ b/Makefile
@@ -44,7 +44,7 @@ AOBJS=      bam_aux.o bam_index.o bam_plcmd.o sam_view.o bam_fastq.o \
             bam_tview.o bam_tview_curses.o bam_tview_html.o bam_lpileup.o \
             bam_quickcheck.o bam_addrprg.o bam_markdup.o tmp_file.o \
             bam_ampliconclip.o amplicon_stats.o bam_import.o bam_samples.o \
-            bam_consensus.o consensus_pileup.o reference.o
+            bam_consensus.o consensus_pileup.o reference.o tanuki.o
 LZ4OBJS  =  $(LZ4DIR)/lz4.o
 
 prefix      = /usr/local
@@ -206,6 +206,7 @@ sam_opts.o: sam_opts.c config.h $(sam_opts_h)
 sam_utils.o: sam_utils.c config.h $(samtools_h)
 sam_view.o: sam_view.c config.h $(htslib_sam_h) $(htslib_faidx_h) $(htslib_khash_h) $(htslib_kstring_h) $(htslib_thread_pool_h) $(htslib_hts_expr_h) $(samtools_h) $(sam_opts_h) $(bam_h) $(bedidx_h)
 sample.o: sample.c config.h $(sample_h) $(htslib_khash_h)
+tanuki.o: tanuki.c
 stats_isize.o: stats_isize.c config.h $(stats_isize_h) $(htslib_khash_h)
 stats.o: stats.c config.h $(htslib_faidx_h) $(htslib_sam_h) $(htslib_hts_h) $(htslib_hts_defs_h) $(samtools_h) $(htslib_khash_h) $(htslib_kstring_h) $(stats_isize_h) $(sam_opts_h) $(bedidx_h)
 amplicon_stats.o: amplicon_stats.c config.h $(htslib_sam_h) $(htslib_khash_h) $(samtools_h) $(sam_opts_h) bam_ampliconclip.h
```

以上で準備は終了である。

## ビルドする

上に書いたビルドコマンドを実行する。
ただしすでに一回実行していたら`make`だけでもいけると思う（多分）。

```
autoheader
autoconf -Wno-syntax
./configure
make
```

## tanuki コマンドを実行する

`./samtools tanuki`
※ `samtools tanuki` としないように注意しよう。`./samtools tanuki` とすることで現在のディレクトリのsamtoolsを呼ぶことができる。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/42e8191c-7ded-b4d0-1ac9-1efc8496e9e4.png)



無事にタヌキが表示された。Samtoolsでタヌキが表示された！

はて、それに何の意味があるのだろう？？？

私にもわからない。人生にはわからない謎がたくさんある。なぜタヌキコマンドなのか。いつかきっと明らかにしたい。

## おわりに

このようにして、Samtools にタヌキを表示するコマンドが追加され、いつでもタヌキを表示ことができるようになった。
太陽は東から登り、西へ沈むようになった。
地球に長い安定と平和が訪れた。

タヌ！