Github: https://github.com/stschiff/msmc2
introduction (Just for reference): https://link.springer.com/protocol/10.1007/978-1-0716-0199-0_7
call SNPs
...
generate mask file
bcftools genomecov --ibam my.bam -bg|awk '{$4=NULL;print $0}'>genomecov.cal
or use vcf2bed implemented in bedops to generate mask file
https://github.com/bedops/bedops/releases/tag/v2.4.39
separate chromosomes
nohup sh -c 'for i in $(<chr.list);do { bcftools view -r ${i} ../vcf_bcftools.SNPs.1.vcf.gz|gzip >vcf_bcftools.SNPs.${i}.vcf.gz; }& done' &
mask
nohup sh -c 'for i in $(<chr.list);do { generate_multihetsep.py --mask ../genomecov.cal --chr ${i} vcf_bcftools.SNPs.${i}.vcf.gz>vcf_bcftools.SNPs.${i}.multihetsep.txt; }& done' &
run
msmc2 -t 27 -o all.msmc vcf_bcftools.SNPs.*.multihetsep.txt
plot
data<-read.table("all.msmc.final.txt",header=T)
mu<-1e-8
gen<-10
pdf("msmc.plot.pdf")
plot(data$left_time_boundary/mu*gen, (1/data$lambda)/(2*mu),log="x",ylim=c(0,100000),
type="n", xlab="Years ago", ylab="effective population size")
lines(data$left_time_boundary/mu*gen, (1/data$lambda)/(2*mu), type="s", col="red")
legend("topright",legend=c("NLE"), col=c("red"), lty=c(1,1))
dev.off()
A good tutorial: https://wurmlab.github.io/genomicscourse/2016-SIB/practicals/msmc/msmc-tutorial/guide
