DEV Community

Cover image for Show a Tanuki with Samtools!
kojix2
kojix2

Posted on

Show a Tanuki with Samtools!

↑ Tanukis (raccoon dogs) at Obihiro Zoo (Wikimedia)

    __,-─-、__
   (〆-─-ヽ)
    ( ´・ω・` )  < Show a Tanuki with Samtools!
   /  ,r‐‐‐、ヽ
    し l  x )J
   _.'、 ヽ  ノ.人
 (_((__,ノU´U. (酒)

Enter fullscreen mode Exit fullscreen mode

Introduction

For those of you who came to this page hoping for a useful and serious bioinformatics article. I'm really sorry. If you opened this page expecting some useful information, please leave quietly. Nothing like that. Read only if you don't mind wasting time on silly content. I have warned you. (From Tanu Chronicle)

In this post, I'm going to add the command tanuki to Samtools to display a ASCII art(AA) tanuki.

Statements of Need

It has long been a dream of all humanity to display a tanuki on the terminal whenever you want. I would like to offer one solution here today.

Finding ASCII art of tanuki

If you want to display a tanuki in your terminal, you must find a tanuki AA. First, get a tanuki AA from this site.

First, get the tanuki ASCII art from this page.

    __,-─-、__
   (〆-─-ヽ)
    ( ´・ω・` )
   /  ,r‐‐‐、ヽ
    し l  x )J
   _.'、 ヽ  ノ.人
 (_((__,ノU´U. (酒)

Enter fullscreen mode Exit fullscreen mode

↑ Tanuki

Tanuki

↑Shigaraki ware tanukis (raccoon dogs) (Wikimedia)

Building samtools

First, let's make sure that Samtools can be built successfully. Note that the build requires htslib and other dependencies. In my case, I put htslib directory in the same directory as samtools because I build htslib by myself for htslib binding of Crystal language which I make for my personal interest. If you put htslib here, samtools can detect it automatically. htslib requires submodules git submodule update -i --recursive. If you get stuck, please google it. It would be safer to check out a stable tag git checkout 1.16.1, and switch to a new branch git switch tanuki. It is up to you.

Then, just follow the instructions in README.md on GitHub to build it.

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
Enter fullscreen mode Exit fullscreen mode

Where is the main function?

Well, the first thing to figure out is where is the main function that is invoked when you type the samtools command.
If you open the samtools repository in VSCode and search for int main(, you will find main in the file bamtk.c.

In the middle of this file you will find help, which can be displayed by typing samtools help.

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");
Enter fullscreen mode Exit fullscreen mode

Add tanuki command help there. Of course, tanuki is an extremely important command, so many people may refer to the help file to find out how to use it. Help is important.

--- 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");
 }
Enter fullscreen mode Exit fullscreen mode

Then, scrolling further down, there is obviously code for conditional branching.

    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;
}

Enter fullscreen mode Exit fullscreen mode

Add main_tanuki here. Now samtools tanuki will call main_tanuki().

--- 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();
Enter fullscreen mode Exit fullscreen mode

Scroll up and add a prototype declaration for the main_tanuki() function.

--- 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);
Enter fullscreen mode Exit fullscreen mode

samtools can be built without it, but if it is not there, the following warning message will be displayed at compile time, so it is better to declare it properly.

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);
      |                                                    ^~~~~~~~~~~
Enter fullscreen mode Exit fullscreen mode

implement main_tanuki()

The tanuki command is very important, and we expect many users. Some of them may try to change a raccoon dog into a peacock or something. In order to prepare for such a future, the copyright information should be clearly indicated in the header.

Display a tanuki with printf.

/*  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"
"    し l  x )J\n"
"   _.'、 ヽ  ノ.人\n"
" (_((__,ノU´U. (酒)\n");

    return 0;
}

Enter fullscreen mode Exit fullscreen mode

Edit Makefile

Finally, edit the Makefile. add tanuki.o to AOBJS. Write that tanuki.o depends on tanuki.c.

--- 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
Enter fullscreen mode Exit fullscreen mode

This is the end of the preparation.

Build

Execute the build command written above.
However, if you have already done it once, you can just use make (probably).

autoheader
autoconf -Wno-syntax
./configure
make
Enter fullscreen mode Exit fullscreen mode

Run the tanuki command.

. /samtools tanuki.
Be careful not to use samtools tanuki. . /samtools tanuki to call samtools in the current directory.

Image description

The tanuki was successfully displayed.
Finally got the tanuki to show up with Samtools!

Well, what does that mean?

I don't know either. There are many mysteries in life. Why the tanuki command? I want to make it clear someday.

Conclusion

Thus, the command to display a tanuki was added to samtools, and it is now possible to display a tanuki at any time.
The sun began to rise in the east and sets in the west.
A long period of stability and peace came upon the earth.

May the fluffiness of the tanuki be with you!

Top comments (0)