利用模擬 PCR 評估 16S rRNA 基因特定變異區的長度分布

我在「怎麼處理 16S rRNA 標的基因分析中長度不一的序列？」提到可以利用模擬 PCR 來評估16S rRNA 基因特定變異區的長度分布。本文示範如何使用 SILVA 資料庫的 TestPrime 模擬 PCR 擴增 16S rRNA 基因，還有如何將 V4 區域用 Usearch 擷取出來以便統計長度分布。

1. 模擬擴增

進入 TestPrime 後，先輸入順逆向引子 (primer) 的序列，再選擇要模擬擴增的資料庫及類別，接著點擊 "Run TestPrime" 便會開始擴增。由於目標是了解 16S rRNA 基因 V4 區域的長度分布，所以要選擇 SSU r138 資料庫 (small subset units)。

因為 LTP ("The All-Species Living Tree" Project) 存放的是被鑑定過的細菌序列，所以裡頭的條目不像基於演算法預測分類的資料那麼多。因此，我選擇使用 LTP 以節省時間。

2. 加入購物車

模擬擴增完畢後，點擊 "Inspect Results in Taxonomy Browser" 察看結果，再將想要下載的分類群之基因序列納入購物車。此處我選擇下載真細菌的擴增序列。

3. 設定下載格式

點擊購物車的標示後，右上角會顯示加入購物車的項目數量，點擊 "Download" 即可設定下載格式，設定完畢後點擊 "Start export" 就會開始匯出資料並且轉進下載頁面。

4. 下載資料

等待匯出完成，選擇項目後點擊 "Download file" 即可下載資料。

5. 取出擴增標的區域

從 TestPrime 下載到的是含有 V4 區域的序列全長，需要使用其他軟體把 V4 區域擷取出來。我用的是 Usearch 內建的 search_pcr2。而使用 QIIME2 擷取標的區域的方法可參考 "Can we use extracted reference reads to calculate expected amplicon sizes?"

usearch -search_pcr2 Silva_SSU138_LTP.fasta\
  -fwdprimer GTGCCAGCMGCCGCGGTAA \
  -revprimer GGACTACHVGGGTWTCTAAT \
  -minamp 0 -maxamp 1500 \
  -strand both \
  -fastaout "Silva_SSU138_LTP_V4.fasta"

6. 繪製長度分布

取得 V4 區域的序列後就可以統計長度分布了，我使用的是 R 語言的 seqinr 套件來統計序列長度。

library(seqinr)
library(tidyverse)
f <- read.fasta("Silva_SSU138_LTP_V4.fasta")
data <- as.data.frame(table(getLength(f)))
ggplot(data = data, aes(x = Var1, y = Freq)) +
  scale_y_continuous(breaks = seq(0, 14000, 2000)) +
  geom_col() +
  labs(x = "Length", y = "Counts") +
  theme_classic() +
  theme(legend.position = "none", 
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20),
        axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20)
  ) 

16S rRNA 基因 V4 區域的長度集中在 250-256 bp 之間