用ANGSD和Stairway Plot2构建种群历史动态

1、用angsd将bam文件转换为saf格式

比起easySFS.py脚本来说,angsd对于missing data的考虑更加全面,但是处理起来相对而言会更麻烦一点。

1
2
3
4
5
6
7
#用samtools给参考基因组构建索引
samtools faidx ref.fa

angsd -bam bam.bamlist \
-doSaf 1 -out XXX \
-anc ref.fa \ #确保上一步生成的fai索引文件位于参考基因组的同一路径下
-GL 1 -minMapQ 1 -minQ 30

2、用angsd内置程序realSFS将saf转换为sfs

这一步总是显示占用过多内存,尤其是在单个种群的样本量多于10个并且都是测序深度较高的时候,相当麻烦。

1
2
3
4
realSFS \
XXX.saf.idx \
-maxIter 100 -P 2 -fold 1 > XXX_folded.sfs
#这个“-fold 1”的参数是为了生成folded版本的sfs。公认的解释似乎是不清楚祖先态的时候适用folded,而vcf文件被较好polarized的时候适用unfolded版本

3、在sfs基础上生成shell批处理文件

首先需要给每个目标种群编辑一份blueprint后缀的文件,里面标注出这个种群的相关信息。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
#XXX blueprint file
#input setting
popid: XXX # id of the population (no white space)
nseq: XX # number of sequences
L: 1234895725 # total number of observed nucleic sites, including polymorphic and monomorphic
whether_folded: true # whethr the SFS is folded (true or false)
SFS: XXX.XXX XXX.XXX XXX.XXX XXX.XXX XXX.XXX XXX.XXX XXX.XXX XXX.XXX XXX.XXX XXX.XXX XXX.XXX # snp frequency spectrum: number of singleton, number of doubleton, etc. (separated by white space)
#smallest_size_of_SFS_bin_used_for_estimation: 1 # default is 1; to ignore singletons, uncomment this line and change this number to 2
#largest_size_of_SFS_bin_used_for_estimation: 17 # default is nseq/2 for folded SFS
pct_training: 0.67 # percentage of sites for training
nrand: (nseq-2)/4 (nseq-2)/2 (nseq-2)*3/4 nseq-2 # number of random break points for each try (separated by white space)
project_dir: /path/to/XXX/ # project directory
stairway_plot_dir: /path/to/stairway_plot_v2.1.1/stairway_plot_es # directory to the stairway plot files
ninput: 200 # number of input files to be created for each estimation
#random_seed: 6
#output setting
mu: XXXe-X # assumed mutation rate per site per generation
year_per_generation: XXX # assumed generation time (in years)
#plot setting
plot_title: demo_XXX # title of the plot
xrange: 0.1,10000 # Time (1k year) range; format: xmin,xmax; "0,0" for default
yrange: 0,0 # Ne (1k individual) range; format: xmin,xmax; "0,0" for default
xspacing: 2 # X axis spacing
yspacing: 2 # Y axis spacing
fontsize: 12 # Font size

编辑blueprint文件参数的相关说明:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
popid: 这里不能包含空格
nseq: 也就是你这个种群个体数×2
L: 参考基因组总长度(总共有多少个碱基)
whether_folded: 按照上面生成sfs文件的过程填true或者false
SFS: 上一步生成的sfs文件中去除第一个数字和之后所有的0
#smallest_size_of_SFS_bin_used_for_estimation: 1 # default is 1; to ignore singletons, uncomment this line and change this number to 2
#largest_size_of_SFS_bin_used_for_estimation: 17 # default is nseq/2 for folded SFS
pct_training: 似乎默认的就是0.67
nrand: 这里的4个数字分别按照这样计算(除不尽时取整数):(nseq-2)/4, (nseq-2)/2, (nseq-2)*3/4, nseq-2
project_dir: 输出结果文件夹
stairway_plot_dir: stairway_plot_es程序所在的路径
ninput: 每次模拟输入的数据量,一般为200
#random_seed: 6
#output setting
mu: 突变率
year_per_generation: 世代数
#plot setting
plot_title: 输出图片的名字
xrange: x轴(距今时间)的范围
yrange: y轴(Ne大小)的范围,默认为0,0
xspacing: X axis spacing
yspacing: Y axis spacing
fontsize: 字体大小

然后运行stairbuilder程序,在这个blueprint文件的基础上生成后续模拟分析需要的批处理脚本:

1
2
nohup java -cp /path/to/stairway_plot_v2.1.1/stairway_plot_es \
Stairbuilder /path/to/XXX_folded.blueprint 2>&1 &

这一步之后会在当前文件夹里生成一个叫XXX.blueprint.sh的脚本。之后只要nohup执行这个脚本就行了。

4、运行stairwayplot程序

1
nohup sh XXX_folded.blueprint.sh > XXX_folded.log 2>&1 &

生成的.final后缀文件就是最终的结果图。