2021.01.22
MSMC2流程

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

2021.01.22 MSMC2流程

2021.01.22
MSMC2流程