msmc2 | 遗传群体史推测

本教程假设相关的软件已经安装好，参考：
https://github.com/stschiff/msmc2
https://github.com/stschiff/msmc-tools/blob/master/msmc-tutorial/guide.md
https://lh3lh3.users.sourceforge.net/snpable.shtml

准备mask文件

编码区mask

为了避免编码区CDS对群体推测的影响，我们一般要准备一个mask文件，指明哪些基因组区域是CDS区。

输入文件：参考基因组的GTF或者GFF文件

GTF="~/idata/input_dat/ref_genome/dmelanogaster/dmel-all-r6.42.gtf"
CHRS=('2L' '2R' '3L' '3R')

awk '{if($3=="CDS"){print $1,$4,$5}}' ${GTF} | gzip > cds.bed.gz
for chr in ${CHRS[@]};do
    echo ${chr}
    zcat cds.bed.gz | awk -v cc=${chr} '$1==cc' | gzip > cds_${chr}.bed.gz
done

注意，为了方便后续操作，我们是把各个染色体独立分开存储的，后面的有些程序不支持含有多个染色体的文件。

可比对区mask

这也主要是针对参考基因组的，比如，某些序列可能是重复序列，那么可比对性就比较差一下。我们此处需要生成一个bed文件，说明参考基因组那些区域是比较好的。

参考：https://lh3lh3.users.sourceforge.net/snpable.shtml，但是他这个教材写的很不完整。大概就是把参考基因组分成小段的reads，然后把这些reads比对回参考基因组，看看哪些区域能够比较特异性的比对。

makeMappabilityMask.py主要是将mask的基因组序列，转换成bed格式的文件，使用时注意修改一下脚本中的输入和输出文件路径，该脚本源代码把路径写死到代码中了！我修改了该脚本，使用sys.argv把输入文件改为命令行参数。

最后，mask文件以bed格式分染色体保存。

export REF="~/idata/input_dat/ref_genome/dmelanogaster/dmel-all-chromosome-r6.42.fasta"

# https://lh3lh3.users.sourceforge.net/snpable.shtml
~/seqbility-20091110/splitfa ${REF} 35 > ref.fq  
bwa aln -R 1000000 -O 3 -E 3 -t 32 ${REF} ref.fq | \
    bwa samse ${REF} - ref.fq > ref.sam

~/seqbility-20091110/gen_raw_mask.pl ref.sam > ref.sam.mask

~/seqbility-20091110/gen_mask -l 35 -r 0.5 ref.sam.mask > ref.sam.mask.fa

# edited input and output path in source file - makeMappabilityMask.py
# regions with '3' in bed file: unique mapped sites
~/msmc2/msmc-tools/makeMappabilityMask.py ref.sam.mask.fa

ls | grep '^ref.sam.mask.fa.*.mask.bed.gz' | xargs cat > ref.mask.bed.gz
ls | grep '^ref.sam.mask.fa.*.mask.bed.gz' | xargs rm 
for chr in ${CHRS[@]};do
    echo ${chr}
    zcat ref.mask.bed.gz | awk -v cc=${chr} '$1==cc' | gzip > ref_${chr}.mask.bed.gz
done

测序质量mask

上面两个mask主要是针对参考基因组的。此处的mask是针对测序样本的，包括测序深度，比对质量评分。这儿测序深度，我们使用样本的平均深度，0.5 - 2倍平均深度的区域是我们要保留的区域；比对质量评分使用默认参数，MapQ=20.

样本可以是单个unphased的二倍体；也可以是trio，通过子代对亲代进行Phase。在推测群体变化的时候，这个情况可能使用到了不同的算法。具体的不同，本人尚未深入研究。

此处，我们使用trio测序数据。为了方便计算，此处将过程写成了bash函数，进行了parallel，所以代码看起来复杂一下，具体可以简化。

输入文件：各个样本的bam文件 输出文件：mask.bed.gz的文件，vcf.gz文件

export SPECIES='dmeN'
export FAM='15'

export FID1_Oi="sid_${SPECIES}_${FAM}_M1" # 子代ID
export FID1_M0="sid_${SPECIES}_${FAM}_M0" # 父亲ID
export FID1_F0="sid_${SPECIES}_${FAM}_F0" # 母亲ID

BAMS=("~/idata/call_variants/${SPECIES}/fam${FAM}/${FID1_Oi}.sort.bam" \
    "~/idata/call_variants/${SPECIES}/fam${FAM}/${FID1_M0}.sort.bam" \
    "~/idata/call_variants/${SPECIES}/fam${FAM}/${FID1_F0}.sort.bam" )

# bamCaller.py does not support multiple chromsomes
# 创建参数文件，方便parallel: 样本 + 染色体
rm bam_chr.txt

for bam in ${BAMS[@]};do
    for chr in ${CHRS[@]};do
        echo ${bam} ${chr} >> bam_chr.txt
    done
done


Step1 ()
{
    SID=$(basename ${1} | cut -d '.' -f 1)
    CHR=${2}
    mean_depth=$(samtools depth -r ${CHR} ${1} | awk '{sum += $3} END {print sum / NR}')
    bcftools mpileup -q 20 -Q 20 -C 50 -r ${CHR} -f ${REF} ${1} | \
        bcftools call -c -V indels | \
        ~/msmc2/msmc-tools/bamCaller.py ${mean_depth} ${SID}_${CHR}_mask.bed.gz | \
        gzip -c > ${SID}_${CHR}.vcf.gz
}

export -f Step1 

cat bam_chr.txt | parallel --colsep ' ' -j 8 Step1

输入文件

上面的过程中，准备好了三个mask文件：编码区，可比对区，测序质量

此处，我们通过generate_multihetsep.py脚本来创建msmc2的输入文件。

--mask：声明我们的mask文件，包括各个样本的测序质量的mask；以及参考基因组可比对区的mask

--negative_mask：这个放编码区的mask. 主要是因为bed文件中声明的编码区是我们想剔除掉的；而测序质量的mask和参考基因组可比对区mask是我们想保留的。

--trio：声明样本编号：<child_index>,<father_index>,<mother_index>

下面的脚本中，我们同样对染色体进行了parallel.

# <child_index>,<father_index>,<mother_index>
Step2 ()
{
    chr=${1}
    ~/msmc2/msmc-tools/generate_multihetsep.py \
        --mask ${FID1_Oi}_${chr}_mask.bed.gz \
        --mask ${FID1_M0}_${chr}_mask.bed.gz \
        --mask ${FID1_F0}_${chr}_mask.bed.gz \
        --mask ../ref_${chr}.mask.bed.gz \
        --negative_mask ../cds_${chr}.bed.gz \
        --trio 0,1,2 \
        ${FID1_Oi}_${chr}.vcf.gz ${FID1_M0}_${chr}.vcf.gz ${FID1_F0}_${chr}.vcf.gz \
        > ${chr}.multihetsep.txt
}

export -f Step2

echo ${CHRS[@]} | tr ' ' '\n' | parallel -j 4 Step2

运行

msmc2居然是用D语言编写的？？我下载了一个编译完成的版本，msmc2_Linux。下载下来，修改文件为可执行文件，即可直接用。

上述分染色体运行的结果是每条染色体都有一个${chr}.multihetsep.txt文件。注意！此处不要把这些染色体合并到一个文件中！！

而应该以不同的独立文件在命令行中，传递给脚本。

-t表示使用CPU核心数；-p表示时间节段，此处用了 2+25+2=29个时间节段，前两个和最后两个时间节段合并估计。-I参数是我们在*.multihetsep.txt文件中要使用的单倍体的数量。本例中只有两个亲本的4条单倍体（0-3）,也可以不用。还有其他对于unphase数据的用法，具体参考软件使用教程。

~/msmc2/msmc2_Linux \
    -t 6 \
    -p 1*2+25*1+1*2 \
    -I 0,1,2,3 \
    -o ${SPECIES}_${FAM} \
    *.multihetsep.txt

bootstrap

为了检验群体史推测的稳健性，可以在上面生成的multihetsep.txt中做bootstrap，通过这些随机抽样，检验模型结果是否稳定。

multihetsep_bootstrap.py的具体参数不再赘述，本例使用如下。

## bootstrap
mkdir bootstrap && cd bootstrap 
~/msmc2/msmc-tools/multihetsep_bootstrap.py \
    --nr_bootstraps 20 \
    --chunk_size 1000000 \
    --chunks_per_chromosome 20 \
    --nr_chromosomes 4 \
    test \
    ../*.multihetsep.txt

# 收集bootstrap的结果文件，然后运行msmc2
for ll in $(find . -regex ".*/bootstrap_multihetsep.chr.*.txt$");do
    ss=$(echo ${ll/.\//} | tr '/' '_')
    ~/msmc2/msmc2_Linux \
        -t 6 \
        -p 1*2+25*1+1*2 \
        -I 0,1,2,3 \
        -o ${ss} \
        ${ll}
done

可视化

这儿需要两个参数，mu, 物种的自发突变率；gen，世代时间。具体数值应该查阅相关文献。对于不同物种的突变率情况，可以参考 https://www.biorxiv.org/content/10.1101/2023.01.24.525323v3：

Wang, Y., & Obbard, D. J. (2023). Experimental estimates of germline mutation rate in eukaryotes: a phylogenetic meta-analysis. bioRxiv, 2023-01.

library(data.table)
mu <- 4.8e-9
gen <- 0.1
# 收集boostrap运行的结果文件
ffs <- list.files(pattern='test.*final.txt', recursive=TRUE)
dmel <- fread("dmel_30.final.txt", header=TRUE) # 修改成你的结果文件
plot(dmel$left_time_boundary/mu*gen, (1/dmel$lambda)/(2*mu), 
  log="x",ylim=c(0,2000000),
     type="n", xlab="Years ago", 
     ylab="effective population size")

for(ff in ffs){
    aa <- fread(ff, header=TRUE)
    lines(aa$left_time_boundary/mu*gen, (1/aa$lambda)/(2*mu), type="s", col="grey")
}

lines(dmel$left_time_boundary/mu*gen, (1/dmel$lambda)/(2*mu), type="s", col="red")

于是我们有了如下结果：

转自：https://mp.weixin.qq.com/s/01WFwCkwz_FbeYuWSAre1Q