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

はじめに

　ロングリードで得られたシークエンスデータを持ちいて Taxonomy 推定 を行うときに、全長16S を用いるべきか、V1-V5 でトリミングして用いるべきか、迷ったのでトリミングツールを実装しました。V5 領域のユニバーサルプライマーの配列に対して、re モジュールを用いた正規表現による相同配列検索、またはBio.pairwise2 を用いたローカルアライメントによる相同配列検索により、V5 領域末端の検出を目指します。fasta データはNCBIからダウンロードしたものを使用します。

実装

# 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)

# 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()

　ローカルアライメントを用いた方法では、計算時間の短縮のため、アライメントする配列に対して600番目の塩基から1200番目の塩基を切り出し、プライマーとの比較を行います。

# 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()

　正規表現との結果を比較します。

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()

おわりに

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