細菌の16S rRNA配列や真菌のITS配列のように、配列を既知のものと照合するデータベース（greengenes, Silva, UNITEなど）が充実しているものは、確認しなくてもいいかもしれない。
しかし、例えば特定の酵素を持つ微生物種の群集構造を調べるために、その酵素をコードした遺伝子配列を使ってNGSをすることがある。
このとき、その配列を指定するためのプライマーが、ときどき対象の遺伝子ではない非特異的な配列を増やしてしまうことがある。
この非特異的なものまで含めて、対象の生物の多様性などを議論するのは見当違いだと思う。そこで、自分のデータ内にある配列が対象の生物由来かをチェックする。
1. チェックの前に
シークエンシングなどで得たOTU（ASV）のFASTA形式の配列を得る。
QIIME2でこの公式チュートリアルに沿って解析しているのであれば、rep-seqs.qzaを以下のコードで.fastaに変換できる。以下のコードを実行するとdna-sequences.fastaというファイルが生成される。これはもうFASTA形式なので、2.までスキップしていい。
qiime tools export --input-path rep-seqs.qza --output-path .R上でDADA2のチュートリアルに沿って解析しているのであれば、途中で出てくる"taxa"に配列がある。今回はtaxaをtaxonomy.txtというファイルに保存して使った。
配列をどうFASTA形式にするかという問題だが、今回はLinuxで行った。うまくやれば何でもできると思う（自分はRでもやったことがある）。
taxonomy.txtは、1行目がKingdom, Phylum, Class, ...のように分類群が書かれており、2行目以降は、1列目に配列、2列目以降に分類名が記されている。
いま欲しいのは配列だけなので、「2行目以降の1列目」が必要だ。
そしてそれぞれの配列に対して、FASTA形式になるように番号と「＞」を振っていく。
# 2行目以降の1列目を取得
tail -n +2 taxonomy.txt | cut -d ' ' -f 1 > test.txt
# 行数を数える
wc -l test.txt
# 例えばここで160行あるとわかるとする
# ASV番号と対応する配列を交互に並べていくことを160回繰り返す
# 何回やるかは配列数（上で調べた行数）による
# sequences.fastaに保存
for i in `seq 160`; do echo ">ASV"$i; sed -n ${i}p test.txt; done > sequences.fasta
# 中間生成物test.txtを消す
rm test.txt出力されるsequences.fastaはこんな感じ。
2. BLASTでチェック
"blastn"とGoogleなどで検索し、BLASTのサイトに飛ぶ。
先ほど作ったsequences.fastaをアップロードしてもいいし、内容をコピペしてもいいので、Enter Query Sequenceのところに入力し、下のBLASTのボタンを押す。
自分が160配列ぶんやったときは、数十秒で結果が出た。
NCBI上の、相同性が高いと思われる配列の候補が下に出てくる。ここで対象の生物らしきものが増えているか確認できる。
また、"Results for"のところのプルダウンで違う配列の結果も確認することができる。
しかし、何百配列もある場合、いちいちプルダウンを選択するのはめんどくさい。そこで、プルダウンの上の"Download All"で何かできないか模索した。
1配列につき100個の相同なデータが出てくるのは多すぎるので、右下の"show"を100から10に変更し、"Download All"のプルダウンから"Hit Table(csv)"を押した。今のところ、これが一番自分好みだ。
1列目から、配列名、アクセッション番号、配列一致度、…と続いていく。数値について詳しくはこちらのサイトに載っている。
アクセッション番号は、その番号でNCBI上で検索すればその配列の素性がわかるというものだ。
今回はそれを利用して、その配列に一番近い配列のタイトルを得ようと思う。
ExcelのVLOOKUPが1番最初に照合された値しか返さないことを利用し、各配列10個あるのを1番近いものだけに絞った。またわかりやすくするために各パラメータに対しヘッダーを振った。この行列をsamples.csvとして保存した。
3. Pythonでスクレイピング
探せばもっといい方法があるのかもしれないが（とくにDownload Allのところが怪しい）、このアクセッション番号を用いて、なんの配列なのかというデータを取得していきたい。
だいぶ前置きが長くなったが、やっとメインテーマだ。
※ここからはPythonです。
# 必要パッケージ類の読み込み
import numpy as np
import pandas as pd
import csv
from bs4 import BeautifulSoup
from urllib import request
import re
import requests# データ読み込みと確認
df = pd.read_csv('samples.csv')
print(df)# 全部0のdescription列を追加
df = df.assign(description = 0)
print(df)# NCBIのNucleotideの検索ページを指定しておく。
url = 'https://www.ncbi.nlm.nih.gov/nuccore/'
# 行数のぶんだけ上の検索ページでアクセッション番号を検索し、出てきたページのタイトルを取得。それを7列目（description列）に入れていくことを繰りかえす。
for i in range(0,len(df.index)):
response = requests.get(url + df.iloc[i,1])
soup = BeautifulSoup(response.text,'html.parser')
element = soup.title.string
df.iloc[i,6] = element# 保存
df.to_csv('sample_description.csv', index=False)以上の操作で各配列が何に近いのかざっくりと理解することができた（description列）。
対象の生物でない場合はそのASV（やOTU）配列を除外しその後の解析に進むといいと思う。
Comments