【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的输出结果
# 刘乐乐 liulele622@163.com
# 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的侧翼序列就可以允许一定的碱基错配，或许能增加你的结果数量。