ローカルアライメントで16SrRNA 塩基配列をトリミングする

# はじめに
　ロングリードで得られたシークエンスデータを持ちいて `Taxonomy 推定` を行うときに、全長16S を用いるべきか、V1-V5 でトリミングして用いるべきか、迷ったのでトリミングツールを実装しました。V5 領域のユニバーサルプライマーの配列に対して、`re` モジュールを用いた正規表現による相同配列検索、または`Bio.pairwise2` を用いたローカルアライメントによる相同配列検索により、V5 領域末端の検出を目指します。`fasta` データは[NCBIからダウンロードしたもの](https://qiita.com/jtnutrient/items/eca85d687787aa1a2600)を使用します。
# 実装
```py
# Forward primer (515F): GTGYCAGCMGCCGCGGTAA
# Reverse primer (806R): GGACTACNVGGGTWTCTAAT
# Reverse Primer (926R): CCGYCAATTYMTTTRAGTTT

import numpy as np
from Bio import Entrez, SeqIO, Align, pairwise2
from Bio.Seq import Seq
from Bio.pairwise2 import format_alignment
import re
import matplotlib.pyplot as plt

# Write IUPAC code base sequence in regular expressions
# https://www.ddbj.nig.ac.jp/ddbj/code.html
def expand_iupac(dna_seq):
    iupac_dict = {
        "A": "A",
        "C": "C",
        "G": "G",
        "T": "T",
        "R": "[AG]",  # A or G
        "Y": "[CT]",  # C or T
        "S": "[GC]",  # G or C
        "W": "[AT]",  # A or T
        "K": "[GT]",  # G or T
        "M": "[AC]",  # A or C
        "B": "[CGT]", # C or G or T
        "D": "[AGT]", # A or G or T
        "H": "[ACT]", # A or C or T
        "V": "[ACG]", # A or C or G
        "N": "[ACGT]" # any base
    }
    return ''.join(iupac_dict[nuc] for nuc in dna_seq.upper())

# primer, 926R
primer_seq_r = "CCGYCAATTYMTTTRAGTTT"
primer_rev_r = "".join(list(reversed(primer_seq_r)))
expanded_primer_r = expand_iupac(primer_rev_r)
```
```py
# primer matching using re module
path = r"file_path"
record = SeqIO.parse(path, 'fasta')
end_re = []
count = 0
for i, content in enumerate(record):
    seq = content.seq
    complement_seq = seq.complement()
    matches = re.finditer(expanded_primer_r, str(complement_seq))
    count += 1
    for match in matches:
        end_re.append(match.span()[1])
        
    
plt.hist(end_re, bins = 140)
plt.title(f"V5 reverse primer matching, re module")
plt.ylabel("count")
plt.xlabel("end base of seq which match primer seq")
plt.show()
```
<img src = "https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3230781/b349d827-fe21-2817-53fc-b78fcfa10f70.png
" width ="500">

　ローカルアライメントを用いた方法では、計算時間の短縮のため、アライメントする配列に対して600番目の塩基から1200番目の塩基を切り出し、プライマーとの比較を行います。
```py
# primer matching using pairwise2.localms module
path = r"file_path"
record = SeqIO.parse(path, 'fasta')
end = []
match = 2
mis_match =  -1
gap_init = -1
gap_extend = -0.5
trim_init = 600
trim_term = 1200
list_alignments = []


for i, content in enumerate(record):
    seq = content.seq
    complement_seq = seq.complement()
    # match; +2, mis_match: -1, gap_initiation: -1, gap_extention: -0.5 localms
    alignments = pairwise2.align.localms(primer_rev_r, str(complement_seq)[trim_init:trim_term],match, mis_match, gap_init, gap_extend) #, match, mis_match, gap_init, gap_extend

    # best alignment score
    best_alignment = max(alignments, key=lambda x: x.score)
    end.append(best_alignment.end)
    list_alignments.append(alignments)

plt.hist((np.array(end)+trim_init).flatten(), bins=140)
plt.title("V5 reverse primer,  localxx alignment")
plt.ylabel("count")
plt.xlabel("end base of seq which match primer seq")
plt.show()
```
　正規表現との結果を比較します。
```py
plt.hist(end_re, bins=140, label="re module")
plt.hist((np.array(end)+trim_init).flatten(), bins=140, label="localms")
plt.legend()
plt.ylabel("count")
plt.xlabel("end base of seq which match primer seq")
plt.show()
```
<img src = "https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3230781/db47e25c-904a-6c65-25c0-31b61aa51fc8.png" width ="500">

# おわりに
　実務では、実際に読まれたシークエンス配列に対して、同様にV5領域を検出できるのかを確認します。どなたかのお役に立てば幸いです。