a) qiime2で系統推定する場合sklearn/blast/vsearchから選択できる
  b) sklearnの場合は参照データとして分類機が必要。blast/vsearchの場合はfastaと系統データが必要。
  c) qiime2で利用する場合はわざわざqiime2でインポートしてqza形式に変換したものを利用しなければならない。カスタムDBを作る場合fastaと系統データをq2用にフォーマットする必要がある。

1 qiime2 リファレンスデータを入手する

1.1 RDP Classifier Files

# DLしてきたファイル名を変更
mv RefOTUs.fa RDPClassifier_16S_trainsetNo18_RefOTUs.fa
mv Ref_taxonomy.txt RDPClassifier_16S_trainsetNo18_taxonomy.txt

# DLしてきたものをそのままインポートしようとしたらエラー。小文字はだめらしい。
# Invalid character 'g' at position 0 on line 2 (does not match IUPAC characters for this sequence type). Allowed characters are ACGTRYKMSWBDHVN

# 配列部分を大文字に変換する。
cat RDPClassifier_16S_trainsetNo18_RefOTUs.fa \
| awk '/^>/{print}!/^>/{print toupper($0)}' > RDPClassifier_16S_trainsetNo18_RefOTUs.fasta

# qiime2でインポートqza形式にする

## fasta配列
qiime tools import \
--type FeatureData[Sequence] \
--input-path RDPClassifier_16S_trainsetNo18_RefOTUs.fasta \
--output-path RDPClassifier_16S_trainsetNo18_RefOTUs.qza

## 系統データ--source-format は--input-formatに変更
qiime tools import \
--type FeatureData[Taxonomy] \
--input-format HeaderlessTSVTaxonomyFormat \
--input-path RDPClassifier_16S_trainsetNo18_taxonomy.txt \
--output-path RDPClassifier_16S_trainsetNo18_taxonomy.qza

# トレーニング`fit-classifier-naive-bayes`, fit-classifier-sklearn開発目的で公開されている
qiime feature-classifier fit-classifier-naive-bayes \
--i-reference-reads RDPClassifier_16S_trainsetNo18_RefOTUs.qza \
--i-reference-taxonomy RDPClassifier_16S_trainsetNo18_taxonomy.qza \
--o-classifier RDPClassifier_16S_trainsetNo18.classifier.qza

1.2 GreenGene 16S

  • qiime2 Data resourcesから分類機を入手可能
  • fastaと系統データも入手可能
  • GreenGene 13_8 99% OTUs の全長配列でトレーニング済みのbayes分類機 gg-13-8-99-nb-classifier.qza
wget -c ftp://greengenes.microbio.me/greengenes_release/gg_13_5/gg_13_8_otus.tar.gz

1.3 Silva 16S/18S

  • qiime2 Data resourcesから分類機を入手可能
  • fastaと系統データも入手可能
  • Silva (16S/18S rRNA)の全長配列でトレーニング済みのbayes分類機 silva-138-99-nb-classifier.qza
  • Silva (16S/18S rRNA)の全長配列および系統データ(qza形式) https://docs.qiime2.org/2022.2/data-resources/

1.4 カスタム NCBI 16S

# taxonkit 最新版をDL
grep "^>" 220804_16S_refseq.fasta | sed 's/^>//' \
| awk '!/Candidatus/{printf("%s", $1"\t"$2" "$3"\t"); for(i=4;i<NF;++i){printf("%s ", $i)}; print $NF} \
/Candidatus/{printf("%s", $1"\t"$2" "$3" "$4"\t"); for(i=5; i<NF; ++i){printf("%s ", $i)}; print $NF}'

# NCBIから落としてきたfastaのヘッダを編集してtaxonkitに投げてlineageを取得
grep "^>" 220804_16S_refseq.fasta | sed 's/^>//' \
| awk '!/Candidatus/{printf("%s", $1"\t"$2" "$3"\t"); for(i=4;i<NF;++i){printf("%s ", $i)}; print $NF} \
/Candidatus/{printf("%s", $1"\t"$2" "$3" "$4"\t"); for(i=5; i<NF; ++i){printf("%s ", $i)}; print $NF}' \
| taxonkit name2taxid -i 2 | taxonkit lineage -i 4 | taxonkit reformat -i 5 -f "{k};{p};{c};{o};{f};{g};{s}" \
| awk -F"\t" '{x=$2" "$3; print $1"\t"x"\t"$4"\t"$6}' > ncbi16s_tax.txt

# １つの配列に複数のtaxonomic idが振られていたので個別に確認
grep -e "NR_153679.1" -e "NR_132333.1" -e "NR_116359.1" -e "NR_118271.1" -e "NR_115954.1" -e "NR_112559.1" ncbi16s_tax.txt > temp.tsv
cat ncbi16s_tax.txt | awk -F"\t" '$3!~/^(1530043|1602736|415956|170600
0|60447)$/{print}'


# taxidつかなかったもの取り出して編集
cat ncbi16s_tax.txt | awk -F"\t" '$4==""{print}' | sed  -e 's/\[//g' -e 's/\]//g' -e  s/\'//g | awk -F"\t" '{print $2}' | cut -f1,2  -d " " | taxonkit name2taxid | taxonkit lineage -i 2 | taxonkit reformat -i 3 -f "{k};{p};{c};{o};{f};{g};{s}"