はじめての環境DNA
研究室から現場まで、実践で使える知識を!
データの解析方法

Qiime2とBLAST+を用いた環境DNAのデータ解析

QIIME2はアンプリコンシーケンスのデータ解析を中心に現在もサポートが続いているデータ解析パイプラインです。またQIIME2と聞くと菌叢解析を想像しますが専用品ではないのでMiseqなどのショートリードシーケンサーから得られたデータであれば好きな解析が可能です。

今回は手持ちのデータがなくても解析できるように、先行研究で配列決定された魚類の環境DNAシーケンスデータをSRA toolkitで取得後、QIIME2でデータ解析をしたのちにBlast+で種のアサイメントを実施します。

SRA toolkitを用いた配列取得とQIIME2での解析、Blast+での種のアサイメントは独立なので、例えば配列取得がうまく行かなくてもfastq.gzさえ手元にあればQIIME2以降の解析は可能です。

動作環境

  • Ubuntu 22.04 LTS
    (多分WSL, WSL2 + Ubuntuでも同様に動作すると思います)
  • python 3.8
  • SRA toolkit : github
    • fasterq-dump version 3.0.2
    • prefetch version 3.0.2
  • BLAST+ version 2.13: github

SRA toolkitは以前の記事からversionが上がっている点に注意してください。最新のものに置き換えて使用されるPCにインストールください。

またDesktop上にqiime2というフォルダを作成して、そこを中心にファイル生成を行います。

デスクトップ上で、ターミナルを開いてqiime2フォルダや、その他の必要フォルダの作成をしておきましょう。

Miniconda3を使ってQiime2をインストール

パッケージマネージャーとしてcondaを提供するMiniconda3を使ってQIIME2をインストールします(QIIME2の推奨法)。Minicondaのインストールは以下の記事を参考に自分の環境に合ったものを選択してください。

QIIME2のインストールはQIIME2のチュートリアル(https://docs.qiime2.org/2022.11/install/native/)に沿って進めます。2022年12月末現在のQIIME2のバージョンは2022.11です。

Minicondaのアップデート

とりあえずcondaのアップデートを行います。ターミナルを開いて以下のコマンドでアップデート。

# All requested packages already installed.と表示されていたらアップデートは終わってます。

mambaのインストール

チュートリアルにはありませんが、パッケージマネージャーとして優秀なmambaをインストールしておきます。

wgetのインストール

ダウンローダーとしてwgetを使用するので、mambaでwgetをインストールしておきます。

qiime2 core 2022.11 distributionをcondaの仮想環境にインストール

以前にQiime2をインストールしたことがある、もしくはインストールを試みたことがある場合も含め、各リリース専用の仮想環境を構築して運用することを推奨します。

以前のバージョンは2022.8でしたが2022.11になっていくつか機能追加がされているようです(次のリリースは2023.2)。

では、仮想環境としてqiime2-2022.11を作成していきましょう。まず、wgetでパッケージバージョンなどが記載されたymlファイルをダウンロードします。

lessなどで中身を見ると、qiime2-2022.11-py38-linux-conda.ymlにはQIIME2での解析に必要な基本的なライブラリがバージョンとともに記載されています。

ymlファイルの情報を使って、2022.11バージョンの仮想環境を作成します。

特にエラーなく以下のような指示が出たらOK。

では、指示に従って仮想環境に入っておきましょう。

仮想環境を出るときは、mamba deactivate、仮想環境がどんな名前だったか確認したい時はmamba info -eで確認できます。最後に使用したymlファイルはいらないので削除。

ここれでQIIME2のインストールは完了。

後で使用するライブラリをインストールしておく
  • parallel : 並列処理 HP
  • seqkit : 配列データ処理 github
  • pigz : 並列解凍・圧縮 HP
  • seqfu : QIIME2のmanifestファイル作成 github

次はSRA toolkitで先行研究で配列決定されたデータをダウンロードします。

SRA toolkitを使ってハイスループットシーケンスデータを取得する

今回は2022年に宮先生が市民科学×環境DNAというテーマで実施された、MiFishプライマーによる魚類群集構造解析論文のデータを使って進めていきたいと思います。

沿岸海域における地域の生物多様性を評価するための魚の eDNA メタバーコーディングにおける市民科学の使用: パイロット研究
The use of citizen science in fish eDNA metabarcoding for evaluating regional biodiversity in a coastal marine region: A pilot study
DOI : 10.3897/mbmg.6.80444

文中には以下のようにあるので、NCBIのSRAデータベースより登録されているデータのアクセッション番号(SRR〜)を確認します。

All raw DNA sequence data and associated information were deposited in the DDBJ/EMBL/GenBank database and are available under accession number DRA012840.

確認したところ、DRA012840には6つのデータがあってDRR321580~DRR321585が割り振られているのでその情報を使って配列データを取得します。

prefetch によるSRAデータのダウンロード

prefetchはSRAデータをダウンロードしてくるコマンドです。一度.sraファイルをダウンロードしておくとfastqへの変換に失敗してしまったりした際にやり直しが楽です。

DRR321580~DRR321585の配列データを一括でダウンロードします。

実行するとoutputフォルダが作成され、その中にDRR~.sraという形式で6データがダウンロードされます。次はfastqerq-dumpを使ってダウンロードした.sraファイルからFastqファイルを生成します。

fasterq-dumpによるsraファイルからのfastq生成

先程outputフォルダに生成した.sraファイルをもとにfastqファイルを00fastqgzフォルダに生成します。

オプションを2つ使用しています。-pは進捗状況を表示、-eは複数のコアを使用しての並列処理オプションです。

_1.fastqと_2.fastq はイルミナのHTSから出力された状態だとR1, R2にあたるファイルです。FowardプライマーとReverseプライマーで読まれた配列になります。

  • DRR32158X_1.fastq
  • DRR32158X_2.fastq

.fastqファイルはpigzで圧縮した後に名前を_1 →_R1、_2 → _R2に変換しておきます。

00fastqgzフォルダの中身がDRR32158X_R1.fastqとDRR32158X_R2.fastqとなっていればOKです。

QIIME2によるデータ処理

ここからは、SRA-toolkitでダウンロードした配列データを用いてQIIME2によるデータ解析処理を進めていきます。ざっくりいうと以下のような感じです。

  1. FastqファイルをQIIME2へインポート
    • seqfu
    • qiime tools import
  2. Demultiplex (シーケンスをサンプルに配分)
  3. 配列上の非生物学的部分(ex. primer)の除去
    • qiime cutadapt trim-paired
  4. Quality scoreの可視化
    • qiime demux summarize
  5. クオリティフィルタリング、デノイジング、キメラ除去
    • qiime dada2 denoise-paired
    • qiime feature-table filter-samples
    • qiime feature-table filter-seqs
  6. データのエクスポート
    • qiime tools export
    • biom convert

Manifestファイルの作成とQiime2へのデータのインポート

QIIME2へfastqファイルをインポートする際、ファイル情報をまとめたmanifestファイルが必要になります。

MGIやPacBioなどのシーケンサーから出力されたfastqファイルはどうか分かりませんが、最近Illunaのシーケンサーから取得したようなfastqファイルであれば、–typeにはSampleData[PairedEndSequencesWithQuality]を–input-formatにはPairedEndFastqManifestPhred33V2を指定すれば良いかと思います。

Read typePhred score–type–input-format備考
Single-endPHRED 33SampleData[SequencesWithQuality]SingleEndFastqManifestPhred33V2
Paired-endPHRED 33SampleData[PairedEndSequencesWithQuality]PairedEndFastqManifestPhred33V2
Single-endPHRED 64SampleData[SequencesWithQuality]SingleEndFastqManifestPhred64V2
Paired-endPHRED 64SampleData[PairedEndSequencesWithQuality]PairedEndFastqManifestPhred64V2
Paired-endPHRED 64SampleData[PairedEndSequencesWithQuality]CasavaOneEightSingleLanePerSampleDirFmtfastqファイルの形式が.+_.+_L[0-9][0-9][0-9]_R[12]_001\.fastq\.gzである必要あり
(URL : https://docs.qiime2.org/2022.11/tutorials/importing/#fastq-manifest-formats)

指定できる値は以下のコマンドで確認できます。

Manifestファイルはawkコマンドを使用して作成することもできますがseqfuを使うと簡単です。両者はManifestファイルの形式が少し異なる点に注意です。

開発者githubはこちら。オプションは以下の通り。

Manifest.tsvは以下のようにな形式になっているはずです。

sample-idforward-absolute-filepathreverse-absolute-filepath
DRR321580$HOME/Desktop/qiime2/00fastqgz/DRR321580_R1.fastq.gz$HOME/Desktop/qiime2/00fastqgz/DRR321580_R2.fastq.gz
DRR321581$HOME/Desktop/qiime2/00fastqgz/DRR321581_R1.fastq.gz$HOME/Desktop/qiime2/00fastqgz/DRR321581_R2.fastq.gz
DRR321582$HOME/Desktop/qiime2/00fastqgz/DRR321582_R1.fastq.gz$HOME/Desktop/qiime2/00fastqgz/DRR321582_R2.fastq.gz
DRR321583$HOME/Desktop/qiime2/00fastqgz/DRR321583_R1.fastq.gz$HOME/Desktop/qiime2/00fastqgz/DRR321583_R2.fastq.gz
DRR321584$HOME/Desktop/qiime2/00fastqgz/DRR321584_R1.fastq.gz$HOME/Desktop/qiime2/00fastqgz/DRR321584_R2.fastq.gz
DRR321585$HOME/Desktop/qiime2/00fastqgz/DRR321585_R1.fastq.gz$HOME/Desktop/qiime2/00fastqgz/DRR321585_R2.fastq.gz
$HOMEの部分はDesktopまでのfull pathが表記される

ちなみに、先程のURL先のページを参考にawkコマンドでManifest.csvが作成できます。

上記コマンドで作成したManifest.csvは以下のようにな形式になっているはずです。

sample-idabsolute-filepathdirection
DRR321580$HOME/Desktop/qiime2/00fastqgz/DRR321580_R1.fastq.gzforward
DRR321580$HOME/Desktop/qiime2/00fastqgz/DRR321580_R2.fastq.gzreverse
DRR321581$HOME/Desktop/qiime2/00fastqgz/DRR321581_R1.fastq.gzforward
DRR321581$HOME/Desktop/qiime2/00fastqgz/DRR321581_R2.fastq.gzreverse
DRR321582$HOME/Desktop/qiime2/00fastqgz/DRR321582_R1.fastq.gzforward
DRR321582$HOME/Desktop/qiime2/00fastqgz/DRR321582_R2.fastq.gzreverse
DRR321583$HOME/Desktop/qiime2/00fastqgz/DRR321583_R1.fastq.gzforward
DRR321583$HOME/Desktop/qiime2/00fastqgz/DRR321583_R2.fastq.gzreverse
DRR321584$HOME/Desktop/qiime2/00fastqgz/DRR321584_R1.fastq.gzforward
DRR321584$HOME/Desktop/qiime2/00fastqgz/DRR321584_R2.fastq.gzreverse
DRR321585$HOME/Desktop/qiime2/00fastqgz/DRR321585_R1.fastq.gzforward
DRR321585$HOME/Desktop/qiime2/00fastqgz/DRR321585_R2.fastq.gzreverse
$HOMEの部分はDesktopまでのfull pathが表記される

seqfu, awkともに今回のデータはPaired-end sequenceデータなのでforwardとreverseの絶対パスが表記されています。

ManifestファイルができたらQIIME2に配列データをインポートします。seqfuで作成したtsv形式のManifestファイルを使用する場合、–input-formatの語尾にV2を、awkコマンドベースで作成したようなcsv形式のManifestファイルの場合は–input-formatの語尾にV2をつけなくて大丈夫です。

Plug in : qiime tools import

と出ていればOKです。

Cutadaptによるプライマー除去とデータの視覚化

以降は、一般的にDADA2によるデノイジングに移ることが多いです。

DADA2は5’末端からプライマー部分の配列が固定長である場合、5’末端から指定した配列長を除去する–p-trun-left-f, –p-trun-left-rオプションでプライマー配列の除去が可能です。ただ、今回は複数プライマー(MiFish-U, U2, Ev2)を使用しているということもあり、cutadaptによってプライマー配列を除去するステップを追加します。

この論文で使用しているプライマーは以下の通り

プライマー名配列(5’→3′)配列長
MiFish-U-F (硬骨魚類)NNNNNNGTCGGTAAAACTCGTGCCAGC27
MiFish-U-R (硬骨魚類)NNNNNNCATAGTGGGGTATCTAATCCCAGTTTG33
MiFish-Ev2-F (軟骨魚類)NNNNNNRGTTGGTAAATCTCGTGCCAGC28
MiFish-Ev2-R (軟骨魚類)NNNNNNGCATAGTGGGGTATCTAATCCTAGTTTG34
MiFish-U2-F (アナハゼ類)NNNNNNGCCGGTAAAACTCGTGCCAGC27
MiFish-U2-R (アナハゼ類)NNNNNNCATAGGAGGGTGTCTAATCCCCGTTTG33
MiFish-Unity-F (コンセンサス配列)NNNNNNGYYGGTAAAWCTCGTGCCAGC27
MiFish-Unity-R (コンセンサス配列)NNNNNNCATAGKRGGGTRTCTAATCCYMGTTTG33

今回は複数プライマーの配列をもとに変異のある部分は縮合塩基に置き換えたコンセンサス配列を指定します。(実際にプライマー配列を指定する際、NNNNNNを入れなかたのは生データをseqfu viewで確認した際、N部分の塩基数がバラバラであったためです:備考の章を参照)

Plug in : qiime cutadapt trim-paired

と出ていればOKです。

使用しているオプションのうち–p-discard-untrimmedは、プライマー配列を認識できなかったリードは廃棄するオプションです。ライブラリ調整や配列決定時のエラーが発生していそうなものは廃棄したいので使用しました。

Quality scoreの可視化

QIIME2は分析に使用するqzaファイルの他に、qiime demux summarizeでqzaファイルを変換した可視化用ファイルqzvが利用できます。生成されたqzvファイルはqiime2 viewで確認できます。

importした配列データ(00demux.qza)とCutadaptで処理した後(01primer-trimmed-demux.qza)のクオリティプロットをqzvファイルを生成して比較してみます。

Plug in : qiime demux summarize

qiime2 viewを開いて、

Drag and drop or click hereの部分に00demux.qzv, 01primer-trimmed-demux.qzvを入れてみると、Bar plotが表示されます。上段のInteractive Quality PlotタブをクリックするとR1, R2のQuality Plotが出てきます。

00demux.qzvのクオリティプロット
01primer-trimmed-demux.qzvのクオリティプロット

プライマー部位が除去されて151bpより短くなっているのが分かります。

DADA2によるデノイジング、エラー配列の修正

DADA2はクオリティスコアをエラーモデルに使用しているので、低品質な配列を切り取っておくとRare variantな配列の感度が向上するようです。

そのため、R1, R2ともに後半の配列(3’側)を–p-trunc-len-f, –p-trunc-len-rオプションで切り捨てます。また、Paired-endシーケンスは、Forwardプライマー側(R1)とReverseプライマー側(R2)つなげる必要があります。

最終的に、以下の図でいうところのRead1, Read2がInsertとかぶる部分の配列を残して結合したいので、ForwardとReverseプライマー配列が除去される(ている)ことと、種判別領域は可変長であることが多いことを念頭に置き、平均産物長をベースに–p-trunc-len-f, –p-trunc-len-rの値を決定します。

まとめると考慮すべき事項は以下の通り。

  • Forward primer length 最大 27bp (NNNNNNを含む)
  • Reverse primer length 最大 33bp (NNNNNNを含む)
  • Overlap length >= 20bp (最低12bpまで可)
  • R1よりR2のほうがクオリティが早く落ちやすいので、–p-trunc-len-f >= –p-trunc-len-rとしておく
  • MiFishの種判別領域の平均産物長は172bp、軟骨魚類やヤツメ類などは180bp程度ある(経験則)

となると最低でも約180bp程度は残るほうが良さそうです。

R1はFプライマーを除くと最低で124bp、R2はRプライマーを除くと最低でも118bpは残る予定です。20bpのOverlap regionを加味しても222bpは残るので少し削って結合しても種判別領域は生成できそうです。

項目配列長残り (おおよそ)
サイクル数151 + 151 = 302bp302bp
Fプライマー27bp275bp
Rプライマー33bp242bp
Overlap20bp222bp

–p-trunc-len-f >= –p-trunc-len-rというルールは守りつつ222bp以上を目指すと、–p-trunc-len-fは115bp, –p-trunc-len-rも115bpとしてDADA2の解析を進めても問題なさそうです。

その他にも色々とオプションはありますが、指定の値以外はデフォルトで実施。

Plug in : qiime dada2 denoise-paired

リードフィルタリング

データを取得元のMiya et al.2020でもやっているので、Singleton, Doubleton, Tripletonを除くために4リード未満の配列を削除します。

Plug in : qiime feature-table filter-samples

リードフィルタリングしたテーブルデータをもとにフィルタリングされて生き残った配列データを出力します。

Plug in : qiime feature-table filter-seqs

データのエクスポート

配列データのエクスポート

配列データはdna-sequences.fastaとして出力指定したフォルダに出力されます。

Plug in : qiime tools export

Summary tableのエクスポート

Table形式のデータはBIOM形式のファイルとしてエクスポートされます。BIOM形式のデータはbiom convertによって.tsv形式のデータへ変換することが出来ます。

Plug in : qiime tools export, biom convert

リードフィルタリングする前のテーブルデータをエクスポート

リードフィルタリングした後のテーブルファイルをエクスポート

左 : フィルタリング前のテーブルデータの最終行 右 : 4リードでフィルタリング後のテーブルデータの最終行

フィルタリングを行う前後で比較すると行数が減っていることが分かります。またこのテーブルデータを参考に出力した、Representative_seq.fastaの配列数もフィルタリング後のテーブルデータのASV数と同数になっているはずです。

ここまでこれば、得られた配列データに対して相同性検索をかけて種を推定することが出来ます。次はBlast+によるBLASTnを実施します。

BLASTnによる塩基配列の相同性検索

まず、魚類の12s rRNAの配列データを用意します。手元に配列データがない場合、2023/1/4にダウンロードしたものがこちらにあります(PMiFish用のデータベース形式になっています)。

データベース用のfastaファイルを02refフォルダに入れて、名前をMiFishDB.fastaに変更しておきます。コマンドで指定するのでそこを変えるなら他の名前でもOKです。

BLAST+のダウンロード

Blast+は配列の相同性検索を行うパッケージです。DNA配列に対してはBLASTnを実行します。mambaを使ってインストールします。

完了したら、確認。

一番下に DESCRIPTION Nucleotide-Nucleotide BLAST 2.13.0+ と出ていたらよいです。インストールできない場合は、必要なcondaのchannelが登録されていないかもしれないので、確認・追加して再度実行してみてください。

Local Blast (BLASTn) の実行

まずfasta形式のデータベースファイルからBLASTn用のデータベースファイル群を作成します。

makeblastdbのオプションは以下の通り。

makeblastdbの実行。

02refフォルダの中には以下の8ファイルが生成されていればOKです。

・MiFishDB.ndb
・MiFishDB.nhr
・MiFishDB.nin
・MiFishDB.njs
・MiFishDB.not
・MiFishDB.nsq
・MiFishDB.ntf
・MiFishDB.nto

次にBlastnを実行します。

オプションは以下の通り。

実行するコードは以下の通り。

色々と指定しているオプションは以下の通りです(NCBI)。

オプション説明
query相同性検索したい配列データ
db検索対象となるデータベース
evalueE-value
lengthquery配列に対するデータベース配列のアライメント長の閾値
num_threadsblast検索時に使用するthread数
perc_identity配列相同性の閾値
qcov_hsp_percquery配列に対してデータベース配列が重複している割合の閾値
bitscoreBit score
max_target_seqs結果を出力する最大配列数
outfmt6を指定するとタブ区切り形式で結果を出力(7はコメントありタブ区切り形式)

出力形式を指定するoutfmtに使用した値に関する情報は以下のとおりです。

オプション説明
qseqidクエリ配列の名前
sseqidデータベース配列の名前
lengthアライメント配列長
mismatcahデータベース配列に対するクエリ配列のミスマッチ数
gapopenデータベース配列に対するクエリ配列のギャップ数
qstartクエリ配列のアライメント開始位置
qendクエリ配列のアライメント終了位置
sstartデータベース配列のアライメント開始位置
sendデータベース配列のアライメント終了位置
evalueE-value
bitscoreBit score
qseqクエリ配列

出力された結果ファイルはタブ区切りテキストなのでエクセルなどを使って確認することが出来ます。

blastres.tsv(ヘッダーは後で付与するので無し)

参照したデータベースの型が悪く、sseqid(2列目)の形式が微妙なので整形します。

出力結果の1行目にoutmftで指定した情報を基本としたヘッダーを作成

blastres_sep_hd.tsv

これで、クエリ配列に対してBlast検索で最大100本までの検索結果の情報を得ることが出来ました。ただ、このままでは同じ種に関する情報が重複していて非常に扱いづらいです。

クエリ配列に対して、Scientific name, Length, Mismatch, Gapopen, E-value, Bitscoreが重複している場合、重複を除去することでスッキリします(データベースがキュレーションされている前提)。

また、出力されているBlastの結果はBit scoreに基づいてSortされており、今回の場合であれば上にある情報ほどクエリ配列の種である可能性が高いです。なので、重複する情報を削除をすることで上位〇〇の配列情報を抽出できるはずです。

今回の結果だと、Queray_seqid:9ec3b8948f811b1372768bd06e9f5d96に対して検索結果が最大100配列分記載されていますが、色のついている部分のMugil cephalusについては、100%一致で情報が重複しています。正しく同定された個体に基づく配列なのであれば、どのAccession numberの情報でもMugil cephalusには変わなさそうです。

なのでこの場合、どれか1本の配列情報に代表させてしまってもMugil cephalusが100%一致でHitしたという事実には変わらないでしょう。以降のMismatchが2や3となっている場合や別種についても同様です。

blastres_sep_hd.tsv

以降は、私がbashスクリプトで処理しきれないためpythonベースでデータ処理していきます。

Blastnで出力された結果から上位10の配列情報に整形する

スクリプト&記事作成中…

まとめ

  • Miniconda3を使ってQIIME2をインストールしました
  • SRA toolkitを使ってハイスループットシーケンスデータの取得を実施しました
  • QIIME2を用いてハイスループットシーケンスデータの解析を実施しました
  • BLAST+をインストールしてQIIME2で得られた配列データに対して相同性検索を実施しました

以降、追記予定

  • 相同性検索の結果、検出された配列に対してその由来である可能性のある種の情報 上位10個をテーブルデータに反映させます。

備考

ダウンロードしたfastqファイルのNNNNNNの数の不一致について

fastq.gzファイルのプライマーのNNNNNN数の不一致について、以下はseqfuのviewオプションを使ってMiFish-UnityをDRR321580_R1.fastq.gzに対して検索かけた結果です。

MiFish-Unity-F
MiFish-Unity-R

>>>はMiFish-Unity-F, <<<はMiFish-Unity-Rを表しています。それより左はNNNNNNに該当する部分ですが、2塩基のときもあれば6塩基のときもあるといった感じです。使用するプライマーのNの数を変えることは配列決定の際、同じ配列が多いようなライブラリでも識別しやすくなるので、そういったことをしているのかもしれません。ただ、プライマー除去時に<NNNNNN + プライマー配列長>だとプライマー配列が残ったり識別領域の配列を削除してしまったりするので今回はCutadaptを使ってプライマー配列を除去しました。

また、DADA2のDocumentationにも配列長を指定して削除するのはプライマーが読み取りの開始点にあり、長さが一定である場合であり、最終的にプライマー配列はしっかりと除去しておくようにとのコメントがあります。

How can I remove primers/adapters/etc from my amplicon reads?

If primers are at the start of your reads and are a constant length (a common scenario) you can use the trimLeft = c(FWD_PRIMER_LEN, REV_PRIMER_LEN) argument of dada2’s filtering functions to remove the primers.

For more complex situations (e.g. heterogenity spacers, variable length amplicons such as the ITS region) we recommend the cutadapt tool or the trimmomatic tool. You can see a worked example, including verification of primer removal, in the ITS-specific version of the DADA2 workflow.

Please double-check that your primers have been removed! The ambiguous nucleotides in unremoved primers will interfere with the dada2 pipeline.

https://benjjneb.github.io/dada2/faq.html

ここらへんはデータ処理の配列処理の注意点なので生データ等をしっかり見て処理を進めたいですね。

BLASTn -outfmtのオプション一覧

Ref

免責事項

本記事の情報・内容については可能な限り十分注意は払っていますが、正確性や完全性、妥当性および公平性について保証するものではありません。本記事の利用や閲覧によって生じたいかなる損害について責任を負いません。また、本記事の情報は予告なく変更される場合がありますので、ご理解くださいますようお願い申し上げます。

ABOUT ME
しばた
分析・データ解析なんでもします。環境DNA使って研究したい方、気軽に相談ください。