五.种群历史与基因交流
书接上回https://www.jianshu.com/p/0147afc3334c
Ne有效群体数量会受到种群大小波动、雌雄数目不均、对后代贡献不均的影响，从而小于真实群体。因此我们可以从有效群体大小来推断整个种群历史过程：
比如在有效群体降低后迅速恢复说明正在遭遇瓶颈效应（驯化过程）
对于野生物种来说，有效群体的降低往往是由于气候变化
PMSC分析需要极高的测序深度，提供突变速率与世代间隔。SMC++软件被更多应用于群体重测序项目中。
treemix利用群体遗传数据计算协方差来构建最大似然树，然后根据估计值和实际值来判断基因流
01.PSMC分析
首先获取单个样本的consensus序列。在参考基因组中变异位点覆盖过低的和过高的都可以质控掉。过高往往是由于重复序列导致的
bcftools mpileup -C 50 -O u -f Chr01.fa X4.Chr01.bam|\ #变异检测
bcftools call \
-c -Ov - | \#获得一致性的变异序列vcf
vcfutils.pl vcf2fq \
-d 10 -D 100 > sample.fq #质控
## 转换成fasta-like格式
fq2psmcfa -q 20 sample.fq > sample.psmcfa
## psmc分析
psmc -N25 -t15 -r5 \
-p "4+25*2+4+6" \ #适用于人的参数
-o sample.psmc sample.psmcfa
## psmc绘图 世代间隔7.5，突变速率2.9e-8
psmc_plot.pl -g 7.5 \
-u 2.9e-8 \
sample.psmc.plot sample.psmc
## eps图片格式转换成pdf
epstopdf.pl sample.psmc.plot.eps sample.psmc.plot.pdf
出图是一张时间和有效群体大小的折线图。接下来做一个个体的自展分析，最后多样本数据整合绘图。
###### bootstrap分析
## 对长序列进行切分，更长的染色体需要切成更大的窗口
$ splitfa sample.psmcfa 100000 > sample.split.psmcfa
$ mkdir sample.bootstrap
## 生成bootstrap分析脚本
seq 100 | while read aa; do echo " \
psmc -N25 -t15 -r5 -b \
-p \"4+25*2+4+6\" -o sample..bootstrap/round-$aa.psmc \
sample..split.psmcfa " ;done > sample..psmc_bootstrap.sh
sh sample.psmc_bootstrap.sh
## 合并结果
cat sample.psmc sample.bootstrap/round-*.psmc > sample.addboot.psmc
## 重新绘图
psmc_plot.pl -g 7.5 -u 2.9e-8 \
sample.addboot.psmc.plot sample.addboot.psmc
## eps转PDF图片
epstopdf.pl sample.addboot.psmc.plot.eps \
sample.addboot.psmc.plot.pdf
#####多样本整合绘图
$ psmc_plot.pl -R \
-M "X4,X12" \ # 显示样品名称
X4-X12.psmc.plot \ # 输出结果前缀
X4.psmc X12.psmc # 输入多个psmc文件
$ epstopdf.pl X4-X12.psmc.plot.eps X4-X12.psmc.plot.pdf
02.SMC++
SMC++主要用于群体的溯祖分析。使用之前需要一个mask文件来标记一些准确度不高的区域（着丝粒区，端粒区等），当做缺失来处理。SNPable可以用于鉴定基因组中大量重复区域
• 文件准备
1.基因组文件
2.群体 vcf 文件
3.突变速率
4.群体分群信息list
## 基因组上滑动取35bp reads 并切分
splitfa genome.fa 35 > read.fa
split -l 20000000 -d read.fa read.split
## 构建基因组bwa比对 index
bwa index genome.fa
## bwa aln才能获得需要的输出
ls ./read.split* | while read aa
do
echo "
bwa aln -t 4 -R 1000000 -O 3 -E 3 genome.fa $aa > $aa.sai
bwa samse genome.fa $aa.sai $aa |gzip > $aa.sam.gz
" > $aa.bwa.sh
done
## 运行比对脚本
ls read.split*.bwa.sh | while read file
do
sh $aa 1>$aa.log 2>$aa.err
done
## 比对结果处理，转成rawMask文件
gzip -dc read.split*.sam.gz | \ #-参数屏幕输出，从而利用|
gen_raw_mask.pl > rawMask_35.fa
## 设置过滤阈值r=0.5
gen_mask -l 35 -r 0.5 rawMask_35.fa > mask_35_50.fa
## 生成mask的基因组文件
apply_mask_s mask_35_50.fa genome.fa > genome.mask.fa
接下来把mask.fa文件中需要屏蔽的区域转换成bed文件，就可以进行SMC++主程序的分析了
第一步的滑动窗口可以变大一点。以现在的测序长度更长，当屏蔽区域太长时可以适当调大参数。
## 准备样品信息 生成pop:sp1,sp2格式文件
$ cat sample.list | awk '$2=="Msie" { sp = sp $1 ","} \
END{print "Msie:" sp } ' | sed 's/,$//' > Msie.list
$ cat Msie.list
Msie:C86,C87,C88,X3,X4,X5,X6,X8
## 一般选2-10个地理分布广、测序质量好的的，作为一个亚群进行似然分析
$ cat sample.list | awk '$2=="Msie"{print $1}' \ | head -n 4 > Msie_sample.list
$ cat Msie_sample.list
C86
C87
C88
X3
$ mkdir Msie_vcf2smc
# 批量生成vcf转smc脚本
### 分样品分染色体运行，双重循环
$ cat chr.list |while read chr
do
cat Msie_sample.list | while read sample
do
echo "smc++ vcf2smc \
-m genome.mask.bed.gz \
-d $sample $sample\ # 指定样品名称
sample.vcf.gz \ # 输入vcf文件名称
Msie_vcf2smc/$sample.$chr.smc.gz \ # 输出smc文件名称
$chr \ # 染色体名称
`cat Msie.list` # 种群名称及样品名称"
done
done > Msie_vcf2smc.sh
$ sh Msie_vcf2smc.sh
###### 进行种群历史有效群体大小估计
##创建分析结果目录
$ mkdir Msie_analysis
## 进行估计分析，生成结果为Msie_analysis/model.final.json
$ smc++ estimate --cores 10 \ # 线程数
--knots 10 \ # 样条节点数，越小越平滑
--spline cubic \ # 默认是piecewise，分段；cubic 平滑线
-o Msie_analysis \ # 输出结果目录
2.9e-8 \ # 每代突变速率
Msie_vcf2smc/*.smc.gz # 输入smc文件
## 结果绘图
smc++ plot -g 7.5 \ # 世代间隔
--linear \ # 纵坐标不取对数
Msie.plot.pdf \ # 输出图片名称
Msie_analysis/model.final.json # 输入json文件
两个亚群
网友评论