五.种群历史与基因交流
书接上回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 并切分
opt/seqbility-20091110/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 vcf2smc_2pop
#同上生成Msyl亚群的文件
$ cat Msie_sample.list Msyl_sample.list > pop2_sample.list
# 批量生成vcf转smc脚本
## 分染色体再分样品运行，双重循环
##### pop1 vs pop2 生成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文件名称
vcf2smc_2pop/pop12.$sample.$chr.smc.gz \ # 输出smc文件名称
$chr \ # 染色体名称
`cat Msie.list` `cat Msyl.list`# 种群名称及样品名称"
done
done > pop12.vcf2smc.sh
$ sh pop12.vcf2smc.sh
## pop2 vs pop1 生成smc文件
cat chr.list |while read chr
do
cat pop2_sample.list | while read sample
do
echo "smc++ vcf2smc \
-m genome.mask.bed.gz \
-d $sample $sample \
sample.vcf.gz \
vcf2smc_2pop/pop21.$sample.$chr.smc.gz \
$chr `cat Msyl.list` `cat Msie.list` "
done
done > pop21.vcf2smc.sh
$ sh pop21.vcf2smc.sh
###### 进行种群历史有效群体大小估计
##创建分析结果目录
$ mkdir split
## 亚群分化时间估计，生成结果为.json文件
$ smc++ split --cores 10 \ # 线程数
--knots 10 \ # 样条节点数，越小越平滑
--spline cubic \ # 默认是piecewise，分段；cubic 平滑线
-o split/ \ # 输出结果目录
Msie_analysis/model.final.json \#参考亚群
Msyl_analysis/model.final.json \#分离亚群
2.9e-8 \ # 每代突变速率
vcf2smc_2pop/*.smc.gz # 输入smc文件
## 结果绘图
smc++ plot --csv
-g 7.5 \ # 世代间隔
--linear \ # 纵坐标不取对数
Msie.plot.pdf \ # 输出图片名称
Msie_analysis/model.final.json # 输入json文件
生成脚本列表，然后批量生成脚本。最后运行主程序，绘图。也可以选取不同的样本来构成不同的亚群，用一样的流程来进行联合绘图。在最后绘图生成的csv文件的最后一行可以看到具体的分化时间，以不同的参考亚群和分离亚群会输出不一样的图。
03.treemix基因交流分析
treemix_pop.list中有样品与亚群的对应关系信息。treemix分析会被连锁位点干扰，所以也需要过滤掉LD位点
## 生成cluster文件，用于plink提取亚群freq信息
$ awk '{print $1"\t"$1"\t"$2}' treemix_pop.list > pop.cluster
$ head -n 4 pop.cluster
sAMBO sAMBO WS
S5 S5 WS
sya-1 sya-1 WS
s11-1 s11-1 WS
$ cut -f 3 pop.cluster |sort -u > pop.order
## 过滤LD
$ plink --vcf treemix.vcf \
--indep-pairwise 50 10 0.2 \ # LD过滤阈值
--out tmp.ld \
--allow-extra-chr \
--set-missing-var-ids @:# \
--keep-allele-order
## 生成freq文件
$ plink --vcf treemix.vcf \
--extract tmp.ld.prune.in \ # 保留的位点信息
--freq \ # 进行亚群freq统计
--missing \ # 生成missing统计信息
--within pop.cluster \ # 指定cluster文件
--out input \
--allow-extra-chr \
--set-missing-var-ids @:# \
--keep-allele-order
$ gzip input.frq.strat
## 将freq文件转成treemix输入格式
$ python2.7 /opt/plink2treemix.py \
input.frq.strat.gz \
input_treemix.frq.gz
#基因交流的次数
$ seq 0 1 3 | while read i
do
echo "
## treemix计算，如果有亚群只包含单个样品, 设置-noss
## 如需要计算p值，设置-se, 但会消耗大量运行时间
treemix -global \
-i input_treemix.frq.gz \ # 指定输入文件
-m $i \ # 设置基因流次数
-root WS \ # 定根
-o out.$i > treemix_$i.log" >treemix_$i.sh
done
六.GWAS分析
P值：性状和表型不相关的概率。由于GWAS是基于连锁不平衡的位点，显著的label也可能只是在性状位点附近，并非一定和性状相关。对数量性状来说，要求其分布尽量符合正态分布。如果不符合的话最好进行数据转换，比如z-score等。异常的表型值需要剔除，它可能发生了一些稀有的变异，不利于GWAS出结果
01.利用SNP的GWAS
·文件准备
vcf 文件: snp.vcf
亲缘关系3.Q矩阵（来自structure分析）
性状文件/表型数据：trait.table
表型值检验
在GWAS和BLUP（最佳线性无偏估计）之前可以对数量性状进行正态性检验，以判断是否出现了异常值。剔除异常值之后可以使用PhenotypeR包来分析多年多点观测的表型平均值、校正值和遗传力。
种群基因型推断
$ mkdir TMP
$ java -Xmx16g -Djava.io.tmpdir=./TMP -jar /opt/beagle.jar\
gt=../data/all.vcf out=all.impute\
impute=true \##是否做基因型推断，F则只输出分型结果
window=10 nthreads=4
beagle本身是进行基因型填充和分型的软件。现在的测序深度足够高，基因型填充反而变得没那么重要了。
GWAS分析
emmax需要tped格式的输入数据，我们使用plink来进行格式转换。
$ plink --vcf ../data/Real_snp.vcf
--maf 0.05 #maf的过滤阈值不能高于0.01
--geno 0.1 #过滤掉缺失率大于0.1的
--recode 12 transpose #需要tped格式的输出
--output-missing-genotype 0
--out snp_filter --set-missing-var-ids @:#
--allow-extra-chr --keep-allele-order
相应的表型数据也需要经过排序，使之和tped文件里面的样品顺序一样，tfam文件中记载了tped的样品顺序。
$ emmax-kin-intel64 -v -d 10 -o ./pop.kinship snp_filter
用emmax里的程序来计算亲缘关系矩阵。输入的前缀名称为snp_filter
$ emmax-intel64 -v -d 10 \
-t snp_filter -p trait.sort.txt -k pop.kinship \
-o emmax.out 1> emmax.log 2>emmax.err
输入亲缘关系矩阵、总snp和表型。输出结果在emmax.ps文件中，主要关注第一列的snp名称和最后一列的p-value。
$ less -S emmax.out.ps|sort -k 4g|less
大多数位点的p值是不显著的。曼哈顿图需要四个输入：snp名称、p-value、snp所在位点和染色体名称。
$ awk '{print $1"\t"$2"\t"$3"\t"$4"\t"}' snp_filter.tped|\
paste - emmax.out.ps|\
awk 'BEGIN{print "SNP\tSHR\tBP\tPvalue"}{if($2 == $5)print {$2"\t"1"\t"$4"\t"$NF} }' >emmax.out.ps.manht_input
从tped文件中取出前四列来获得染色体名称和snp名称、位点。把它和emmax.ps文件按列粘贴起来后，验证$2和$5中的snp名称是否一致后输出$NF列
BEGIN模式添加表头，最后用CMplot包绘制曼哈顿图。这张图如果是矢量图的话，图片的pdf会过大，可以同时输出png用来查看。emmax分析的时候会自己验证亲缘关系和群体结构矩阵，也可以通过-c参数来手动加进去一个Q矩阵
QQplot图也可以在R中输出，来反映理论P-value和实际情况来验证位点。
除了emmax外，tassel软件也可以进行GWAS分析。不过还是emmax更加合适些，省心。
