28.2 Gb基因组SSR序列知多少：Misa+Primer3流程

前两天帮人下载了28.2Gb的蝾螈基因组GCA_002915635.3 ，这么大的基因组，还是第一次近距离接触。

由于之前我优化了下Misa+Primer3 流程脚本，使之可以耗费较小的服务器资源用于分析核心基因组SSR并设计引物。目前为止，我接过的付费分析中还没有超过3Gb的基因组，因此，萌生了对蝾螈基因组分析SSR并设计引物的想法。

就现在网上公开的脚本和修改方法，并不能直接完美将Misa+Primer3 流程应用于核心基因组的分析。如果你有需要分析的，可以直接联系我做付费分析。

流程bash脚本

• 脚本名称run_ssr.sh
• 输入文件为fasta序列，非压缩文件
• 先对fasta序列进行备份儿，加前缀bak_
• fasta序列ID部分空格后的内容需要先删掉。

$ cat run_ssr.sh
#！/usr/bin/env bash

usage() {
cat << EOF
Function: Misa + Primer3
Usage:  bash `basename $0` <fasta file> 
e.g. bash `basename $0` test.fasta
EOF
    exit 1
}

[[ $# == 1 ]] || usage

FASTA=$1

[[ -f bak_${FASTA} ]] || cp ${FASTA} bak_${FASTA} 

sed -i 's/ .*//g' $FASTA

echo "----misa.pl `date`----"
perl misa.pl ${FASTA} 
echo "----p3_in.pl `date`----"
perl p3_in.pl ${FASTA}.misa 
echo "----primer3 `date`----"
primer3_core -default_version=1 -p3_settings_file=my_default_settings.txt -output=${FASTA}.p3out  ${FASTA}.p3in
echo "----p3_out.pl `date`----"
perl p3_out.pl  ${FASTA}.p3out  ${FASTA}.misa
echo "----done `date`----"

运行

• 下载蝾螈基因组GCA_002915635.3_AmbMex60DD_genomic.fna.gz 并解压

nohup bash run_ssr.sh GCA_002915635.3_AmbMex60DD_genomic.fna &

• 结果

时间跨度5月29日00:35至6月1日21:14，将近4天的时间。要是未优化的脚本，如果忽略PB级的存储需求，分析时间估计得以月计算。

$ ll *  |cut -d ' ' -f 5-
 359 5月  27 00:49 run_download.sh
 12K 5月  27 00:52 misa.pl
6.6K 5月  27 00:52 my_bak_default_settings.txt
3.7K 5月  27 00:52 my_default_settings.txt
2.0K 5月  27 00:52 p3_in.pl
4.8K 5月  27 00:52 p3_out.pl
 159 5月  27 00:52 misa.ini
   0 5月  27 00:58 chr.txt
 620 5月  29 00:35 run_ssr.sh
 27G 5月  29 00:38 GCA_002915635.3_AmbMex60DD_genomic.fna
177M 5月  29 05:47 GCA_002915635.3_AmbMex60DD_genomic.fna.misa
971K 5月  29 05:47 GCA_002915635.3_AmbMex60DD_genomic.fna.statistics
2.4G 5月  29 06:45 GCA_002915635.3_AmbMex60DD_genomic.fna.p3in
 14G 5月  29 21:45 GCA_002915635.3_AmbMex60DD_genomic.fna.p3out
1.6G 6月   1 21:14 GCA_002915635.3_AmbMex60DD_genomic.fna.results
246M 6月   1 21:14 nohup.out

结果的简单统计

SSR位点总数

$ cat GCA_002915635.3_AmbMex60DD_genomic.fna.misa |sed '1d' |wc -l
3493294

SSR各基序类型数量

$ cat GCA_002915635.3_AmbMex60DD_genomic.fna.misa |sed '1d' |cut -f 3 |sort |uniq -c
 314475 c
  13481 c*
1970674 p1
 766334 p2
 248903 p3
 144760 p4
  31958 p5
   2709 p6

各染色体上SSR基序数量

• 这里只查看前几条染色体

$ cat GCA_002915635.3_AmbMex60DD_genomic.fna.misa |sed '1d' |cut -f 1 |sort |uniq -c |head
 174815 CM010927.2
 179389 CM010928.2
 168717 CM010929.2
 173822 CM010930.2
 108918 CM010931.2
 189743 CM010932.2
 143619 CM010933.2
 146335 CM010934.2
 156108 CM010935.2
 156510 CM010936.2

统计各染色体的各基序类型的SSR数量

• 这里只展示2条染色体的各基序类型

$ cat GCA_002915635.3_AmbMex60DD_genomic.fna.misa |sed '1d' |cut -f 1,3|sort |uniq -c |sort -k2,2V -k3,3 |head -n16
  15211 CM010927.2	c
    704 CM010927.2	c*
 100919 CM010927.2	p1
  37598 CM010927.2	p2
  12015 CM010927.2	p3
   6777 CM010927.2	p4
   1487 CM010927.2	p5
    104 CM010927.2	p6
  15688 CM010928.2	c
    714 CM010928.2	c*
 100350 CM010928.2	p1
  40087 CM010928.2	p2
  13213 CM010928.2	p3
   7467 CM010928.2	p4
   1734 CM010928.2	p5
    136 CM010928.2	p6

SSR大小分布

$ cat GCA_002915635.3_AmbMex60DD_genomic.fna.misa |sed '1d' |cut -f 5 |sort |uniq -c |sort -k1,1n  |tail -n 10
  64707 24
  94449 18
 103391 16
 106531 20
 137457 13
 204082 15
 221915 14
 442271 11
 593660 12
 911889 10

统计成功设计引物的SSR位点

• 引物列去掉非空行即为成功设计引物的SSR位点

$ cat  GCA_002915635.3_AmbMex60DD_genomic.fna.results |sed '1d' |cut -f 8 |grep -v  '^$'  |wc -l
3420769

引物设计成功率97.92%