种群历史动态研究常用软件有PSMC、MSMC、SMC++
分析案例
PSMC, MSMS, MSC++
MSMC和SMC++是PSMC代表性的改进方法。这两种方法都是突破了PSMC一次只能分析1个样本的局限，而是可以整合并同时分析多个等位基因序列间的最近共祖时间，从而提高了有效群体大小(Ne)预测的精度并提高了效率。由于是多个个体单倍型之间的比较，因此可以提供更多突变时间更近的变异信息(突变是随着时间累积的，所以分化时间越近的序列之间，变异越稀有)，从而提高了这些方法对较近时间点(1万年以内)的Ne的估算精度。
算法原理示意图
1. MSMC方法特点
（1）运算量非常大，一般一次最多只能分析4个个体(8条单倍型)。如果项目中重测序的样本数非常大，MSMC就无法充分使用全部样本的信息。
（2）用于分析的数据必须是相邻等位基因间相位状态已知的基因型数据(phased genotype)。在重测序结果中，每个位点的基因型都是孤立的。比如图2中的基础基因型数据，我们知道两个位点的基因型分别是A/C和C/T。但实际上，这是个二倍体生物，每个位点的两种等位其实是分别分布在两条同源染色体上。
2. MSMC2安装
https://github.com/stschiff/msmc2/releases/tag/v2.0.0
#或者通过GitHub安装
git clone https://github.com/stschiff/msmc-tools.git
#也可以下载压缩文件包zip [here](https://github.com/stschiff/msmc-tools/archive/master.zip)。
分析中会用到以下包含在MSMC2软件中的程序：
cgCaller.py
generate_multihetsep.py
combinedCrossCoal.py
msmc2
3. 测试文件下载
所有测试的输入文件可以在here下载。
其中cg_data子目录中有六个样本信息。前三个样本来自居住在尼日利亚的西非Yoruba的一个父子三人组。其中，NA19240是后代，NA19238和NA19239是两个父代。 后三个样本构成了来自美国犹他州的父子三人组，具有未指定的欧洲血统。其中，NA12878是后代，NA12891和NA12892是父母。
4. 为每个样本生成VCF文件和mask文件
#!/usr/bin/env bash
MASTERVARDIR=/path/to/sequence_data
OUTDIR=/path/to/output_files
CHR=chr1
for IND in NA19238 NA19239 NA19240 NA12878 NA12891 NA12892; do
MASTERVAR=$(ls $MASTERVARDIR/masterVarBeta-$IND-*.tsv.bz2)
OUT_MASK=$OUTDIR/$IND.$CHR.mask.bed.gz
OUT_VCF=$OUTDIR/$IND.$CHR.vcf.gz
cgCaller.py $CHR $IND $OUT_MASK $MASTERVAR | gzip -c > $OUT_VCF
done
在这里仅对1号染色体（chr1）进行分析，还可以在此脚本中循环1至22号染色体。
完成后每个样本生成一个* .mask.bed和一个* .vcf.gz文件。
MSMC是一种隐马尔可夫模型，它使用杂合位点（1/0基因型）的密度来估计到最近的共同祖先的时间。 对于密度，分母是非杂合位点的数量，因此通常是纯合的参考等位基因。 出于效率考虑，这些都不是该VCF文件的一部分。 Mask文件提供了可以调用基因组中哪些区域的信息，覆盖范围不足或质量太低的区域将被排除。
mask文件如下：
chr1 11093 11101
chr1 11137 11154
chr1 11203 11235
5. 准备MSMC2输入文件
使用一个名为generate_multihetsep.py的脚本，该脚本将VCF和mask文件合并在一起，并且还执行简单的三次定相。 如下命令为单个二倍体样本NA12878生成MSMC输入文件：
#!/usr/bin/env bash
$INDIR=/path/to/VCF/and/mask/files
$OUTDIR=/path/to/output_files
$MAPDIR=/path/to/mappability/mask
generate_multihetsep.py --chr 1 --mask $INDIR/NA12878.mask.bed.gz \
--mask $MAPDIR/hs37d5_chr1.mask.bed $INDIR/NA12878.vcf.gz > $OUTDIR/NA12878.chr1.multihetsep.txt
.multihetsep文件如下所示
1 68306 44 TTTCTCCT,TTTCCTTC
1 68316 10 CCCTTCCT,CCCTCTTC
1 87563 13 CCTTTTTT
这是MSMC的输入文件。前两列表示样本中分离位点的染色体和位置。第四列包含我们输入的四个亲本的8个亲本单倍型中的8个等位基因。当有多个模式被逗号分隔时，这意味着定相信息是含糊的，因此有多个可能的定相。
第三列被叫站点数。在上面的第一行中，第一个隔离位点位于位置68306，但并非直到该位点的所有68306位点都被称为纯合参考，而只有44个。这对于MSMC非常重要，因为否则它将假定有一个从1到68306的巨大的纯合片段。
6. 运行MSMC2
(1) 使用MSMC2来估计仅在四个非洲个体和在仅四个欧洲个体内的一致率。
#!/usr/bin/env bash
INPUTDIR=/path/to/multihetsep/files
OUTDIR=/path/to/output/dir
msmc2 -t 6 -p 1*2+15*1+1*2 -o $OUTDIR/EUR.msmc2 -I 0,1,2,3 $INPUTDIR/EUR_AFR.chr1.multihetsep.txt
msmc2 -t 6 -p 1*2+15*1+1*2 -o $OUTDIR/AFR.msmc2 -I 4,5,6,7 $INPUTDIR/EUR_AFR.chr1.multihetsep.txt
参数解读
-t 6: 使用6个CPU。默认情况下，将使用该计算机上的所有CPU。注意：四个个体的一条染色体上，具有六个可能的配对组合，可使用6个CPU。
-p 1 * 2 + 15 * 1 + 1 * 2: 定义时间段模式。默认情况下，MSMC使用32个时间段，分组为1 * 2 + 25 * 1 + 1 * 2 + 1 * 3，这意味着前2个段已合并（强制两个段的合并率相同），然后是25个细分，每个细分都有自己的汇率，然后分别是2个和3个两组。 MSMC2的运行时间和内存使用量随时间段的数量成二次比例增长。在这里，由于我们仅分析单个染色体，因此应减少片段数以避免过度拟合。因此，设置了18个细分，前后分别设置了两组。分组有助于避免过度拟合，因为它减少了可用参数的数量。
-o: 表示输出前缀。其结尾为.final.txt，.loop.txt和.log。
-I: 表示所分析的单倍型从0开始的索引。我们有8个单倍型，前四个是欧洲血统，后一个是非洲血统。在第一轮中，我们估计了欧洲染色体内的聚结率（索引0、1、2、3），在第二种情况下，我们估计了非洲染色体内的聚结率（索引4、5、6、7）。 msmc2的最后一个参数是multihetsep文件。
(2) 使用MSMC2来估计群体分离历史。
#!/usr/bin/env bash
INPUTDIR=/path/to/multihetsep/files
OUTDIR=/path/to/output/dir
msmc2 -t 4 -I 0,1,5,6 -P 0,0,1,1 -s -p 1*2+15*1+1*2 -o $OUTDIR/AFR_EUR.msmc2 $INPUTDIR/EUR_AFR.chr1.multihetsep.txt
参数解读
-t 4: 4个CPU
-I 0,1,5,6: 在每个子种群的前两个亲本染色体上运行
-P 0,0,1,1: 选择的前两个单倍型分别来自子种群0和1。 如果（将花费更长的时间），则必须说
-I 0,1,2,3,4,5,6,7 -P 0,0,0,0,1,1,1,1: 表示要分析所有八种单倍型（可以省略-I标志，默认情况下它使用所有染色体）。
-s: 跳过相位不明确的位点
R画图
mu <- 1.25e-8
gen <- 30
afrDat<-read.table("~/Data/SpringWorkshop/results/AFR.msmc2.final.txt", header=TRUE)
eurDat<-read.table("~/Data/SpringWorkshop/results/EUR.msmc2.final.txt", header=TRUE)
plot(afrDat$left_time_boundary/mu*gen, (1/afrDat$lambda)/mu, log="x",ylim=c(0,100000),
type="n", xlab="Years ago", ylab="effective population size")
lines(afrDat$left_time_boundary/mu*gen, (1/afrDat$lambda)/mu, type="s", col="red")
lines(eurDat$left_time_boundary/mu*gen, (1/eurDat$lambda)/mu, type="s", col="blue")
legend("topright",legend=c("African", "European"), col=c("red", "blue"), lty=c(1,1))
在20万年前之前，两个祖先的人口都有相似的有效人口规模，此后，欧洲祖先经历了严重的人口瓶颈。 这是相对较低的分辨率，因为我们仅分析一个染色体。 请注意，这里我们使用30年的生成时间和1.25e-8的突变率来缩放时间和比率，这与在MSMC上的初始出版物中使用的值相同。
(3) 使用msmc-tools系统信息库中提供的脚本CombineCrossCoal.py来计算总体合并率
#!/usr/bin/env bash
DIR=/path/to/msmc/results
combineCrossCoal.py $DIR/EUR_AFR.msmc2.final.txt $DIR/EUR.msmc2.final.txt \
$DIR/AFR.msmc2.final.txt > $DIR/EUR_AFR.combined.msmc2.final.txt
R画图
mu <- 1.25e-8
gen <- 30
crossPopDat<-read.table("~/Data/SpringWorkshop/results/EUR_AFR.combined.msmc2.final.txt", header=TRUE)
plot(crossPopDat$left_time_boundary/mu*gen, 2 * crossPopDat$lambda_01 / (crossPopDat$lambda_00 + crossPopDat$lambda_11),
xlim=c(1000,500000),ylim=c(0,1), type="n", xlab="Years ago", ylab="relative cross-coalescence rate")
lines(crossPopDat$left_time_boundary/mu*gen, 2 * crossPopDat$lambda_01 / (crossPopDat$lambda_00 + crossPopDat$lambda_11), type="s")
参考
https://wurmlab.github.io/genomicscourse/2016-SIB/practicals/msmc/msmc-tutorial/guide
网友评论