【R語言代碼速記】利用SSRMMD開發SSR引物及結果整理

海量的高通量測序數據為微衛星（SSR）標記開發提供了便利條件，MISA等一大批SSR引物開發軟件被開發出來，但這些軟件要么不能篩選有多態性的位點，要么僅限於linux平台，要么依賴的其他軟件比較多，總而言之，是對菜鳥新手不友好。

幸運的是，我找到了四川農業大學劉亞西教授等去年發布的SSRMMD軟件，僅需要一個perl腳本（SSRMMD.pl）就可以實現多態性SSR位點的篩選，再需要一個腳本（connectorToPrimer3.pl）就可以實現在primer3里開發引物，而且這兩個腳本還被貼心地編譯成exe文件，無需安裝perl也能使用。

不過，其輸出的位點信息文件和引物信息文件是分離的，不太符合自己的需求，故針對自己的數據分析寫了一段R腳本整理結果。

Gou X, Shi H, Yu S, et al. SSRMMD: A Rapid and Accurate Algorithm for Mining SSR Feature Loci and Candidate Polymorphic SSRs Based on Assembled Sequences. Front Genet. 2020;11:706. Published 2020 Jul 27. doi:10.3389/fgene.2020.00706

 

首先，開發針對兩個個體（AU和EU）開發多態性引物，可用  -h 查看兩個腳本命令的幫助文檔

perl SSRMMD.pl -f1 input/AU.fasta -f2 input/EU.fasta -p 1 -o output -mo (2=7,3=6,4=5)
perl connectorToPrimer3.pl -i output/AU.fasta-and-EU.fasta.compare -s 2 -o co_primers.txt -us 100-250

因為數量不符合我的要求，然后分別開發兩個個體的ssr引物

:: AU
perl SSRMMD.pl -f1 input/AU.fasta -p 0 -o output_AU -mo (3=6, 4=5)
perl connectorToPrimer3.pl -i output_AU/AU.fasta.SSRs -s 1 -o AU_primers.txt -us 100-250
:: EU
perl SSRMMD.pl -f1 input/EU.fasta -p 0 -o output_EU -mo (3=6, 4=5)
perl connectorToPrimer3.pl -i output_EU/EU.fasta.SSRs -s 1 -o EU_primers.txt -us 100-250

把序列信息、位點信息和引物信息匯總到一塊，並提取需要的列；此外，把以上三部分位點合並到一塊寫入文件保存。

# 整理SSRMMD的輸出結果
# 劉樂樂 [e-mail]
# 2021年8月15日

# 載入需要的軟件包
library(tidyverse) #優雅地操縱數據
library(Biostrings)#其readDNAStringSet函數用於讀入fasta文件

# 載入個體AU的fasta序列信息
au_fst <- readDNAStringSet("AU.fasta")
au <- paste(au_fst)
names(au) <- gsub(" ", "_", names(au_fst)) #用_代替空格
# 載入另一個個體EU的fasta序列信息
eu_fst <- readDNAStringSet("EU.fasta")
eu <- paste(eu_fst)
names(eu) <- gsub(" ", "_", names(eu_fst))

# 載入SSRMMD.pl的輸出
au_ssr <- read.delim("AU.fasta.SSRs", sep = "\t", header = TRUE)
eu_ssr <- read.delim("EU.fasta.SSRs", sep = "\t", header = TRUE)
co_ssr <- read.delim("co_SSRs.txt", sep = "\t", header = TRUE)

# 載入connectorToPrimer3.pl的輸出
co_primer <- read.delim("co_primers.txt", sep = "\t", header = TRUE)
au_primer <- read.delim("AU_primers.txt", sep = "\t", header = TRUE)
eu_primer <- read.delim("EU_primers.txt", sep = "\t", header = TRUE)

# 整理多態性SSR位點信息
co_res <-co_primer %>% 
  mutate(number = floor(id),
         name = paste0("CO_", number)) %>%
  left_join(co_ssr, by = "number") %>%
  mutate(motif = fasta1_motif,
         repeat_number = paste(fasta1_repeat_number, fasta2_repeat_number, sep = ","),
         seq = au[fasta1_id]) %>%
  select(name, forward_primer, reverse_primer, motif, repeat_number, produce_size, seq, fasta1_id) 

# 整理個體AU的SSR位點信息，隨機篩選150個位點
au_res <- au_primer %>%
  mutate(number = floor(id),
         name = paste0("AU_", number)) %>%
  slice_sample(n = 150)%>%
  select(-id) %>%
  left_join(au_ssr, by = "number") %>%
  mutate(fasta1_id = id,
         seq = au[id]) %>%
  select(name, forward_primer, reverse_primer, motif, repeat_number, produce_size, seq, fasta1_id,)

# 整理個體EU的SSR位點信息，隨機篩選150個位點
eu_res <- eu_primer %>%
  mutate(number = floor(id),
         name = paste0("EU_", number)) %>%
  slice_sample(n = 150)%>%
  select(-id) %>%
  left_join(eu_ssr, by = "number") %>%
  mutate(fasta1_id = id,
         seq = eu[id]) %>%
  select(name, forward_primer, reverse_primer, motif, repeat_number, produce_size, seq, fasta1_id)

# 匯總到一塊
res <- rbind(co_res, au_res, eu_res)

# 去重
#which(duplicated(res$fasta1_id)==TRUE)
#res_c <- res[-c(16, 80),]
#which(duplicated(res_c$fasta1_id)==TRUE)

write_csv(res_c, "ssr_res.csv")

 

更新2021年8月29日，

感謝軟件作者苟香建勘誤和指導。當我們獲得的位點數量不足時，可以使用參數 -me 進行對SSR側翼的序列保守性要求進行調整。

這里我想給您推薦一個參數（-me），可以用來調整檢查SSR側翼序列保守性的算法。由於您在計算中對-me采用的默認設置（NO），這樣設置后，SSR的側翼序列必須是絕對保守的，不允許任何的鹼基錯配/缺失。您可以修改-me的設置，比如設置為NW，即使用Needlman-Wunsch算法來檢查SSR側翼序列保守性，這樣SSR的側翼序列就可以允許一定的鹼基錯配，或許能增加你的結果數量。