Rで解析：キメラ配列を検出・除去できる「rchime」パッケージの紹介

公開日 2026-08-24

アンプリコン解析では、PCR増幅の過程で異なる鋳型に由来する断片が途中でつながり、実在しない「キメラ配列」が生成されることがあります。しかし、キメラ配列を見分けるには、参照配列データベースとの照合やde novo法による判定など、煩雑な手順を踏む必要があります。「rchime」パッケージは、アンプリコンシーケンスデータからキメラ配列を検出し除去できるパッケージです。VSEARCHのアルゴリズムを用いたde novo法による検出と、参照配列データベースを用いた照合による検出のどちらにも対応しています。また、検出に使うパラメータを個別に調整することも可能です。

パッケージバージョンは0.1.2。Windows 11 x64 (build 26200)のR version 4.6.1で確認しています。

Rグラフィックスクックブック第2版 ―ggplot2によるグラフ作成のレシピ集 | Winston Chang, 石井弓美子, 河内崇

ggplot2でのグラフ作成をレシピ形式でまとめた一冊。この図をどう描くか、という引き方ができる。

パッケージのインストール

下記コマンドを実行してください。

# パッケージのインストール
install.packages("rchime")

# パッケージの読み込み
library("rchime")

コマンド例

詳細はコメント、パッケージのヘルプを確認してください。

アンプリコン解析では、PCR増幅時に異なる鋳型由来の断片が途中でつながり、実在しない配列（キメラ配列）が生じることがあります。キメラ配列を除かずに解析を進めると、存在しない分類群を検出したように見えたり、多様性を過大に見積もったりする原因になります。キメラ配列を判定する方法には、既知の参照配列データベースと照合するリファレンス法と、データ中で存在量の多い配列を仮の親候補として使うde novo法があります。rchimeコマンドは、referenceオプションを指定しなければde novo法、指定するとリファレンス法で判定します。

以下の例では、十勝地方の3つの畑（ジャガイモ畑、コムギ畑、テンサイ畑）の土壌から得たと仮定した5種類のASV（配列単位）を用意します。asv_chimera1は、asv_potato_soilの前半とasv_wheat_rootの後半をつなぎ合わせて作成したキメラ配列です。

# 十勝地方の3圃場の土壌から得たと仮定した5種類のASVの配列を用意する
# asv_chimera1はasv_potato_soilの前半とasv_wheat_rootの後半をつないだキメラ配列
配列表 <- data.frame(
  sequence_name = c("asv_potato_soil", "asv_wheat_root", "asv_beet_rhizo",
                    "asv_dairy_grass", "asv_chimera1"),
  sequence = c(
    "TGAACCGCAGTCGTAGGAAATACGAGGGGTGTTAGTACTGTTACCGGATCACGTGTTGTGTGTCGGAGGGAGGGCCAACTGCAAGCACGTCTCGATGGGATATGGCCTCCCTTCGATTAAGCGTGAGCCGCTTCAGGACGGCGTACGGCCGGATTTATCCCGCATACAATCTCTGTGTGGCCCGCCGCCTCGGACCCGTCTCGGGGGTTAGTCGAATTTG",
    "GCCAGGAGTTCCATTTCAAATATCCTCTAAATTTACAAGCGTGAGCCGTGGCGCCCTTGTAACTATGGCTGATGTCTTGTCTCGTATGGAGTTAATTCGAAACGACTCCAGACAGGAGACCAACTGATTCCCACCCCAACTGGACCTTCCACCGGCTGCAGTAGATGTTTTATAGCCCGGGCAGCGACTTCTAGTATAGCGTCGCTAATCCCGAAGTTAG",
    "CTCGCGGCACACAGTCTAGAATGACTCCATAAGCATAAAGCAGAATTGCTCCATTTTAGACAATAGGCCGCACGCATATATTTATAGATATTCCCCTCCATAATGTGTCACTAATAAATAACCACACCCAGAAGTTTCGAAGTCTCCACCTACGTAGCTCAGATGAACTCCCGATTGCAGGTTGCTTGTCACTCGATCTTACAAGCCACACCTCACATGT",
    "GGACAGGGCTGGTTCTGTCGACTCCCCGGGATTGGGGGTTCTCAAATTTTGCTGGCCTGGTAAGAACATTGTCACGACAACTACACCACCCGGCGCGGAGTATGATTTGCCCGATCTTCAGACAGTCTGCCCTTCATAATACGGCGGTTTTGTATACCGAAAAGGTTATGTTCAAAATGCTTTTAAGGAACCGGATACTGACCAGTAGACTTCTGTGTTT",
    "TGAACCGCAGTCGTAGGAAATACGAGGGGTGTTAGTACTGTTACCGGATCACGTGTTGTGTGTCGGAGGGAGGGCCAACTGCAAGCACGTCTCGATGGGATATGGCCTCCGACAGGAGACCAACTGATTCCCACCCCAACTGGACCTTCCACCGGCTGCAGTAGATGTTTTATAGCCCGGGCAGCGACTTCTAGTATAGCGTCGCTAATCCCGAAGTTAG"
  ),
  stringsAsFactors = FALSE
)

# 3圃場それぞれでの各ASVの存在量（リード数）を用意する
存在量表 <- data.frame(
  sequence_name = rep(配列表$sequence_name, times = 3),
  sample = rep(c("tokachi_potato", "tokachi_wheat", "tokachi_beet"), each = 5),
  abundance = c(120, 95, 60, 40, 3,
                80, 130, 55, 30, 2,
                70, 60, 140, 35, 4)
)

# 配列と存在量を1つのデータフレームにまとめる
土壌データ <- merge(配列表, 存在量表, by = "sequence_name")

キメラ配列を検出し除去する：rchimeコマンド

データフレームを指定すると、strollurオブジェクトへ変換しキメラ配列を検出します。table_namesの既定値はsequence_name、sequence、abundance、sampleの列名を想定しています。土壌データの列名はこれと同じなため、table_namesは省略できます。

オプション	意味	初期値
data	strollurデータセットオブジェクト、または配列データを含むdata.frame	なし
reference	strollurデータセットオブジェクト、または参照配列データを含むdata.frame	NULL
dereplicate	サンプル単位でキメラを除去するか。FALSEだと、あるサンプルでキメラと判定された配列は全サンプルから削除される。既定のTRUEでは、キメラと判定されたサンプルからのみ削除される（保守的な判定のため推奨）	TRUE
verbose	コンソールへの出力を許可するか	TRUE
remove_chimeras	検出したキメラをデータから除去するか	TRUE
rchime_options	vsearch固有のパラメータを調整するリスト。rchime_optionsコマンドで作成する	NULL
table_names	dataまたはreferenceがdata.frameのとき、配列名・配列・存在量・サンプル名の列名を対応付ける名前付きリスト	list(sequence_name = “sequence_name”, sequence = “sequence”, abundance = “abundance”, sample = “sample”)

# de novo法でサンプルごとにキメラ配列を検出し除去する（既定かつ推奨の設定）
結果 <- rchime(data = 土壌データ, verbose = FALSE)
結果
            starts ends nbases ambigs polymers numns numseqs
Minimum:         1  220    220      0        4     0    1.00
2.5%-tile:       1  220    220      0        4     0   22.88
25%-tile:        1  220    220      0        4     0  228.75
Median:          1  220    220      0        4     0  457.50
75%-tile:        1  220    220      0        5     0  686.25
97.5%-tile:      1  220    220      0        5     0  892.12
Maximum:         1  220    220      0        5     0  915.00
Mean:            1  220    220      0        4     0  457.64

scrap_summary:
  type      trash_code unique total
1 sequence rchime_chimeras      1     9

Number of unique seqs: 4 
Total number of seqs: 915 

Total number of samples: 3 
Total number of custom reports: 1

出力のscrap_summaryから、1種類・延べ9件の配列がキメラとして除去されたことがわかります。除去後の配列ごとの存在量と、キメラ判定の詳細を確認します。

# 除去後の配列ごとの存在量を確認する
strollur::abundance(結果, type = "sequence")
    sequence_name abundance
1  asv_beet_rhizo       255
2 asv_dairy_grass       105
3 asv_potato_soil       270
4  asv_wheat_root       285

# キメラと判定された根拠（親候補と判定結果）を確認する
strollur::report(結果, type = "chimera_report")[
  , c("Query", "ParentA", "ParentB", "Chimeric_Status")]
            Query         ParentA        ParentB Chimeric_Status
1  asv_beet_rhizo               *              *               N
2 asv_potato_soil               *              *               N
3  asv_wheat_root               *              *               N
4 asv_dairy_grass               *              *               N
5    asv_chimera1 asv_potato_soil asv_wheat_root               Y

asv_chimera1は存在量表から消え、chimera_reportではasv_potato_soilとasv_wheat_rootを親としてChimeric_Statusが「Y」と判定されています。作成時につなぎ合わせた2配列が、そのまま親として検出されました。

＜おすすめのRに関する書籍です＞

R言語徹底解説 | Hadley Wickham, 石田基広, 市川太祐, 高柳慎一, 福島真太朗

Advanced Rの邦訳。オブジェクトの仕組みや関数型プログラミングなど、Rの言語仕様を深掘りする。

検出パラメータを設定する：rchime_optionsコマンド

vsearch固有の検出パラメータをまとめて設定するコマンドです。dereplicate以外は既定値のまま使うことが推奨されています。返り値には、既定値から変更した引数だけが反映されます。

オプション	意味	初期値
processors	使用するコア数の整数	parallelly::availableCores()
dereplicate	サンプル単位でキメラを除去するか。FALSEだと、あるサンプルでキメラと判定された配列は全サンプルから削除される。既定のTRUEでは、キメラと判定されたサンプルからのみ削除される	TRUE
abskew	de novo法でのみ使用する最小の存在量比。abskewは「min(親候補1の存在量, 親候補2の存在量) / 検証対象の存在量」で計算される	2
minh	キメラとして報告する最小スコア。値を下げると検出感度が上がるが、誤検出も増える	0.28
mindiv	最小分岐比率。クエリ配列と、親候補のうち最も近いものとの一致率の差	0.8
xn	「no」投票の重み（βパラメータ）	8
dn	「no」投票数に対する疑似カウントの事前値	1.4
maxp	検討する親候補の最大数	3

# 既定値のまま呼び出す
rchime_options()
$processors
system 
28 

$dereplicate
[1] TRUE

# minh（検出感度）だけを既定値から変更して呼び出す
rchime_options(minh = 5)
$processors
system 
28 

$dereplicate
[1] TRUE

$minh
[1] 5

既定のrchime_options()ではprocessorsとdereplicateのみが表示され、abskewなど残りの引数は既定値のままだと出力に含まれません。minhに5を指定すると、その項目だけが追加されます。processorsの値（今回の環境では28）は実行するマシンのコア数に依存するため、環境によって異なります。

minhはキメラと判定するための最小スコアで、値を大きくすると判定が厳しくなります。先ほどの土壌データにminh = 5を指定し、判定がどう変わるか確認します。

# 検出を厳しくして（minh = 5）キメラ配列を検出し除去する
結果2 <- rchime(data = 土壌データ,
              rchime_options = rchime_options(minh = 5),
              verbose = FALSE)

# 除去後の配列ごとの存在量を確認する
strollur::abundance(結果2, type = "sequence")
    sequence_name abundance
1  asv_beet_rhizo       255
2    asv_chimera1         9
3 asv_dairy_grass       105
4 asv_potato_soil       270
5  asv_wheat_root       285

既定では除去されていたasv_chimera1が、判定を厳しくしたのでキメラとして検出されなくなり、存在量表に残っています。

同梱データのファイルパスを取得する：rchime_exampleコマンド

rchimeパッケージには、動作確認用のstrollurオブジェクトや参照配列データがinst/extdataフォルダに同梱されています。rchime_exampleコマンドは、この同梱データのファイルパスを取得します。fileを省略すると、同梱ファイルの一覧が得られます。

オプション	意味	初期値
file	ファイル名。NULLの場合は同梱ファイルの一覧を返す	NULL

# 同梱データフォルダに含まれるファイルの一覧を確認する
list.files(rchime_example())
[1] "miseq_abundance.rds"                  "miseq_abundance_by_sample.rds"       
[3] "miseq_data_frame.rds"                 "miseq_data_frame_by_sample.rds"      
[5] "miseq_data_frame_by_sample_small.rds" "miseq_fasta.rds"                     
[7] "miseq_names_by_sample.rds"            "miseq_sequences_by_sample.rds"       
[9] "silva.gold.rds"                       "single_sample_abundance.rds"         
[11] "strollur_miseq_tiny.rds"              "strollur_multi_sample_small.rds"     
[13] "strollur_multi_sample_tiny.rds"       "strollur_reference.rds"              
[15] "strollur_single_sample.rds"           "strollur_single_sample_tiny.rds"  

# silva.gold.rdsのファイルパスを取得する
参照パス <- rchime_example("silva.gold.rds")

# 取得したパスにファイルが実在するか確認する
file.exists(参照パス)
[1] TRUE

silva.gold.rdsは、次のsilva_goldコマンドが内部で読み込んでいるファイルです。

改訂2版 RユーザのためのRStudio[実践]入門〜tidyverseによるモダンな分析フローの世界 | 松村優哉, 湯谷啓明, 紀ノ定保礼, 前田和寛

スクレイピングでの収集、dplyr/tidyrでの前処理、ggplot2での可視化、R Markdownでのレポートまで、分析フロー全体を通しで学べる。

SILVA参照配列データベースを取得する：silva_goldコマンド

silva_goldコマンドに引数はありません。呼び出すと、SILVAデータベースに由来する5,181件の参照配列を収めたdata.frameを返します。この返り値をrchimeコマンドのreferenceオプションへそのまま渡すと、リファレンス法によるキメラ検出に利用できます。

# SILVA由来の参照配列データベースを取得する
参照配列 <- silva_gold()

# 行数と列数を確認する
dim(参照配列)
[1] 5181    2

# 列名を確認する
names(参照配列)
[1] "sequence_name" "sequence"  

# 配列そのものは全長で長大なため、先頭の配列IDのみ確認する
head(参照配列$sequence_name)
[1] "7000004128189528" "7000004128189537" "7000004128189547" "7000004128189554" "7000004128189557"
[6] "7000004128189575"

Rによる統計データ解析 | 小池祐太, 村田昇, 吉田朋広

確率論の基礎から推定・検定・回帰までを、数式とRのコードの両方で追える大学教科書。理論をきちんと押さえたい人向け。

この記事が誰かの役に立ちますように。