[argparse]Biopythonの機能をスクリプト化して使いやすくする

動機

Biopythonには本当にいろんな機能が詰まっている。
それらを

• オプションをキーワード引数として受け付ける
• helpメッセージが出る

ようなスクリプトとして使うにはどうすればいいのか。こういうやつ。

$ ./gb2faa.py
usage: gb2faa.py [-h] -i INPUT [-o OUTPUT]

Extract protein sequences from a genbank file

optional arguments:
  -h, --help            show this help message and exit
  -i INPUT, --input INPUT
                        GenBank flat file format of the genomic sequence(s) (required)
  -o OUTPUT, --output OUTPUT
                        output fasta file (default: out.faa)

皆さんには自明なことかもしれないが、私は詰まったので一応ここに書いておきたい。

ArgumentParser (argparse)

PythonにはArgumentParser (argparse) というライブラリがデフォルトでついているので、それでコマンドライン引数を処理させる。

例: gb2faa.py

下記に示すのは、NCBIからダウンロードしたGenBank形式ファイルからタンパク質配列を抽出し、FASTA形式ファイルとして保存するスクリプト。

#!/usr/bin/env python
# coding: utf-8

import sys
import argparse
from Bio import SeqIO
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.SeqFeature import SeqFeature

def _get_args():
    parser = argparse.ArgumentParser(
        description='Extract protein sequences from a genbank file')
    parser.add_argument(
        '-i',
        '--input',
        help='GenBank flat file format of the genomic sequence(s) (required)',
        type=str,
        required=True)
    parser.add_argument(
        '-o',
        '--output',
        help='output fasta file (default: out.faa)',
        type=str,
        default="out.faa")
    if len(sys.argv) == 1:
        parser.print_help(sys.stderr)
        sys.exit(1)
    args = parser.parse_args()
    return args

def main():
    args = _get_args()
    in_gb = args.input
    out_faa = args.output
    records = SeqIO.parse(in_gb, 'genbank')
    out_records = []
    for record in records:
        for feature in record.features:
            if feature.type == 'CDS':
                if 'protein_id' in feature.qualifiers.keys():
                    record_id = feature.qualifiers['protein_id'][0]
                else: 
                    record_id = "{}_{}-{}".format(record.id,feature.location.start,feature.location.end)
                out_record = SeqRecord(
                Seq(
                        feature.qualifiers['translation'][0]),
                    id=record_id,
                    description="{} [{}]".format(
                        feature.qualifiers['product'][0],
                        record.annotations['organism']))
                out_records.append(out_record)
    with open(out_faa, "w") as output_handle:
        for out_record in out_records:
            SeqIO.write(out_record, output_handle, "fasta")

if __name__ == "__main__":
    main()

• -i (--input)引数で指定されたGenBankファイル(args.input)を読み込む
• GenBankファイル内の各DNA/RNA配列をrecordとして読み込む
• 各record からCDS (タンパク質をコードする配列) featureを読み込む
• CDS featureに付随するtranslation attributeからタンパク質配列を抜き出し、out_recordsに格納する
• out_recordsに格納されたタンパク質配列を-o (--output)引数で指定された名前のFASTA形式ファイル (args.output) に書き出す

という作業を行うスクリプトである。
太字部分がargparseの担当。

LC635856.1.gb

LOCUS       LC635856               15969 bp    DNA     circular INV 07-AUG-2021
DEFINITION  Penaeus japonicus Ginoza2017 mitochondrial DNA, complete genome.
ACCESSION   LC635856
VERSION     LC635856.1
DBLINK      BioProject: PRJDB11151
            BioSample: SAMD00276454
            Sequence Read Archive: DRR271272, DRR271273, DRR271276, DRR278743
KEYWORDS    .
SOURCE      mitochondrion Penaeus japonicus
  ORGANISM  Penaeus japonicus
            Eukaryota; Metazoa; Ecdysozoa; Arthropoda; Crustacea;
            Multicrustacea; Malacostraca; Eumalacostraca; Eucarida; Decapoda;
            Dendrobranchiata; Penaeoidea; Penaeidae; Penaeus.
REFERENCE   1
  AUTHORS   Kawato,S., Nishitsuji,K., Arimoto,A., Hisata,K., Kawamitsu,M.,
            Nozaki,R., Kondo,H., Shinzato,C., Ohira,T., Satoh,N., Shoguchi,E.
            and Hirono,I.
  TITLE     Genome and transcriptome assemblies of the kuruma shrimp,
            Marsupenaeus japonicus
  JOURNAL   G3 (Bethesda) jkab268 (2021) In press
  REMARK    DOI:10.1093/g3journal/jkab268
            Publication Status: Available-Online prior to print
REFERENCE   2  (bases 1 to 15969)
  AUTHORS   Kawato,S., Nozaki,R., Hirono,I. and Kondo,H.
  TITLE     Direct Submission
  JOURNAL   Submitted (07-JUN-2021) Contact:Hidehiro Kondo Tokyo University of
            Marine Science and Technology, Laboratory of Genome Science; 4-5-7
            Konan, Minato, Tokyo 108-8477, Japan URL
            :http://www2.kaiyodai.ac.jp/~h-kondo/genome/

$ ./gb2faa.py -i LC635856.1.gb -o LC635856.1.faa

LC635856.1.faa

>BCX71401.1 NADH dehydrogenase subunit 2 [Penaeus japonicus]
MLLSSSQILFFSTLSLGTLLSISSTSWFGAWIGLELNLLSFIPLISSKNNQYSSEAALKY
FLIQALGSSVIIMSASFMMSSNNLATLLLTVALLLKSGAAPFHFWFPSVMEGLQWPQVII
LMTIQKVAPLSLLSYLVSDHILPVIPAAIIMSALVGAIGGVNQTFLRKIVAYSSINHMAW
MMSAILISETDWLLYFSFYCLISSSVAILLSFQEAFHISHIVNHSAYSPEMKMLTFMSLL
SLGGLPPFTGFIPKWIIIQEMVSSGMFIILTVLLASALLTLFYYLRITMTALSISSPKTK
WSMKTAISSNITPSILYTNFFGLLAPSLFILIL
>BCX71402.1 cytochrome c oxidase subunit I [Penaeus japonicus]
TQRWLFSTNHKDIGTLYFIFGAWAGMVGTALSLIIRAELGQPGSLIGDDQIYNVVVTAHA
FVMIFFMVMPIMIGGFGNWLVPLMLGAPDMAFPRMNNMSFWLLPPSLTLLLSSGMVESGV
GTGWTVYPPLSASIAHAGASVDLGIFSLHLAGVSSILGAVNFMTTVINMRSTGMTMDRMP
LFVWAVFITALLLLLSLPVLAGAITMLLTDRNLNTSFFDPAGGGDPVLYQHLFWFFGHPE
VYILILPAFGMISHIISQESGKKEAFGTLGMIYAMLAIGVLGFVVWAHHMFTVGMDVDTR
AYFTSATMIIAVPTGIKIFSWLGTLHGTQLNYSPSLIWALGFVFLFTVGGLTGVVLANSS
IDIILHDTYYVVAHFHYVLSMGAVFGIFAGIAHWFPLFTGLTLNPKWLKIHFLVMFVGVN
ITFFPQHFLGLNGMPRRYSDYPDAYTAWNVVSSIGSTISLIAVLGFIIIVWEALTAARPV
LFSLFLPTSIEWQHNLPPADHSYMEIPLITN
>BCX71403.1 cytochrome c oxidase subunit II [Penaeus japonicus]
MATWGYLGLQDSASPLMEQLIFFHDHAMVVLILITTLVGYMMATLFFNTFTNRFLLEGQT
IEIVWTVLPALILIFIALPSLRLLYLLDEVNNPSVTLKTIGHQWYWSYEYSDFLQVEFDS
YMIPTNELPEDGFRLLDVDNRTILPMNTQIRVLISAADVIHSWTVPALGVKADAIPGRLN
QVSFLMNRPGLFYGQCSEICGANHSFMPIVIESVSVNSFLNWISSASEA
>BCX71404.1 ATP synthase F0 subunit 8 [Penaeus japonicus]
MPQMAPLLWLNLFLMFSATFVMFIILNYFIKVPAKIEKSPSTPQKMEMTWKW
>BCX71405.1 ATP synthase F0 subunit 6 [Penaeus japonicus]
MMTNLFSVFDPTSSVLMLPLNWLSTFLGILFLPMLYWAMPSRWSLLWTMVSTTLHKEFKT
LLGTSHLGTTLMFVSLFSLIVFNNFLGLLPYVFTSSSHLTMTLALALPLWLAFMIFGWVN
HTQHMFAHLVPQGTPGPLMPFMVLIETISNVIRPGTLAVRLAANMIAGHLLLTLLGSTGP
SLSMTLVSILIIGQILLLILEAAVAVIQSYVFAVLSTLYASEVT
...

argparseを使うことで

• helpメッセージが出る (昔書いたスクリプトが何やってんのかわからない、ということがなくなる）
• 引数の順番にとらわれなくなる（イライラしなくなる）
• 応用が広くきく（類似した別研究でも使いやすくなる）
  という利点がある。

参考サイト

https://docs.python.org/ja/3/library/argparse.html
https://atmarkit.itmedia.co.jp/ait/articles/2201/11/news031.html
https://nxdataka.netlify.app/argmemo/
https://github.com/satoshikawato/bio_small_scripts/blob/main/gb2faa.py