针对全基因组数据设计SSR引物

应用软件：samtools（获得染色体全长信息）、bwa（建立索引）、GATK（变异查找与基因分型，需要安装jdk8）、seqkit（fasta/q处理）、Misa（SSR位点搜索）、bedtools（提取序列bed）、primer3（设计引物）、e-PCR（筛选）
软件安装用的conda或者wget命令，以前的文章里面有。

根据基因组大小考虑是否拆分染色体
bedtools2.29在分析大于512M的染色体时，会报错。当染色体过大时，需要拆分染色体。
获取每个染色体行数：

cat Final.re.fa |grep -n ">Chr*" >PLIB.txt             #Final.re.fa是全基因组数据名字

拆分基因组成染色体：
①未拆分染色体：假如得到的第一二行是你1号染色体

sed -n '1,2p' Final.re.fa >Chr1.fa

②拆分了染色体：假如得到的第一二行是你1号染色体的第一部分。三四行是1号染色体的第二部分

sed -n '1,2p' Final.re.fa >Chr1_part1.fa
sed -n '3,4p' Final.re.fa >Chr1_part2.fa

以此类推

建立索引（samtools）
①先cat一个Chr1.bed的空白文件：

cat >Chr1.bed

②统计染色体全长信息：

samtools faidx Chr1.fa

把其中的全长序列信息放入Chr1.bed的空白文件，Tab键分隔。如：

注意，输入.bed文件时，全长要减1，因为是从0计数。

分析SSR（原文链接：//www.greatytc.com/p/d9a9b0f3f09f）
①目标染色体（区段）SSR分析：

perl misa.pl Chr1.fa

②提取SSR上下游各150bp的序列：

cat Chr1.fa.misa |awk 'NR>1 {print $1"\t"$6-150"\t"$7+150}' >Chr1_ssr.bed

③使用bedtools工具提取重复序列：

bedtools getfasta -fi Chr1.fa -bed chr1_ssr.bed -fo Chr1_ssr.fa

④再进行一次SSR查找

perl misa.pl Chr1_ssr.fa

⑤修改p3_in.pl文件，primer3 2.4.0条件下：
print OUT "PRIMER_SEQUENCE_ID= $id"."_$ ssr_nr\nSEQUENCE=$seq\n";

改为

print OUT "SEQUENCE_ID= $id"."_$ ssr_nr\nSEQUENCE_TEMPLATE=$seq\n";
⑥运行misa

perl p3_in.pl Chr1_ssr.fa.misa

⑦primer3进行引物设计：

primer3_core --default_version=1 --output=Chr1_ssr.fa.p3out Chr1_ssr.fa.p3in

对结果进行处理：

perl p3_out.pl Chr1_ssr.fa.p3out Chr1_ssr.fa.misa

⑧e-pcr处理数据：

cat Chr1_ssr.fa.results |awk 'NR>1 {print $1"\t"$8"\t"$11"\t"$1}' >Chr1_ssr.txt

e-PCR输出文件（时间会很久）：

e-PCR  Chr1_ssr.txt D=100-500 /mnt/hgfs/重测序数据/Final.re.fa N=2 G=2 T=3 > Chr1_ssr_ePCR.txt

⑨R studio：
将输出的文件放入本地电脑的文件夹下，我放的H盘
运行以下命令：

setwd("H:\\")          #移动到H盘
data0<-read.table("H:/Chr1_ssr_ePCR.txt",sep="\t",header=F,na.strings="?")        #打开文件
data<-as.matrix(data0)
dup<-unique(data[duplicated(data[,2]),2])
num<-c()
for (k in 1:length(dup)) {                         #开始循环
  n<-which(data[,2]==dup[k])
  num<-c(num,n)
}
res0<-data[-num,]
write.csv(res0,"chr1_ssr_ePCR.csv",row.names=F)          输出结果

合并文件
R studio跑出来的文件就是SSR引物所在位置，需要和原本输出的.result文件合并，就可获得序列
在Excel内用VLOOKUP命令
=VLOOKUP(目标ssr，总ssr+引物，2，0)
然后补充好表格，该染色体ssr数据获取结束。
网页比对数据，检测特异性
http://202.194.139.32/PrimerServer/

2020.8 by ZXG

最后编辑于：2020.08.20 15:16:11