前两天帮人下载了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%