本文為CSDN博主「生信菜鳥(niǎo)1號(hào)」的原創(chuàng)文章,原文鏈接:https://blog.csdn.net/weixin_45694863/article/details/126801831
SMC++是一個(gè)用于從全基因組序列數(shù)據(jù)中估計(jì)種群大小歷史的程序总放,其可以進(jìn)行多個(gè)樣本的分析呈宇,v1.15.4 以后smc++不支持conda了,改為docker運(yùn)行局雄。
1.SMC++安裝
SMC++可以作為Docker image直接運(yùn)行甥啄,其安裝方法也及其簡(jiǎn)單,就是運(yùn)行一遍炬搭,當(dāng)docker找不到SMC++image時(shí)自動(dòng)下載蜈漓,需要一定時(shí)間。
#簡(jiǎn)直不要太簡(jiǎn)單
sudo apt-get install -y python3-dev libgmp-dev libmpfr-dev libgsl0-dev #下載依賴庫(kù)宫盔,我忘記用docker安裝需不需要這些庫(kù)了融虽,但是下載了也沒(méi)啥問(wèn)題。
sudo docker run --rm -v $PWD:/mnt terhorst/smcpp:latest #下載
sudo docker images #看一下docker有哪些image
如下如飘言,MSC++ image下載完畢
看一下具體參數(shù)
sudo docker run --rm -v $PWD:/mnt terhorst/smcpp:latest -h
有效種群大小主要用到vcf2msc衣形,eastimate,plot
2.輸入文件
SMC++的輸入文件可以是vcf文件姿鸿,非常棒谆吴。首先需要得到你分析的群體的vcf文件
vcftoosl --vcf xxx.vcf --recode --keep poplation.txt --out smc #population中時(shí)你的群體樣本
bgzip smc.recode.vcf #壓縮
tabix smc.recode.vcf.gz #為vcf構(gòu)建tbi索引
之后就是格式轉(zhuǎn)換,vcf轉(zhuǎn)smc格式苛预,用到了vcf2msc
for i in {1..19}
do
sudo docker run --rm -v $PWD:/mnt terhorst/smcpp:latest vcf2smc ./smc.recode.vcf.gz ./data/chr${i}.smc.gz $i S:S1,S2,S4,S5,S6,S7,S9
done
#輸入vcf文件句狼,輸出smc.gz文件,$i為染色體热某,樣本的排列是 population name:sample1,sample2,sample3
#值得注意的是mask文件腻菇,在MSMC2中需要,在MSC++也可以添加
3.分析
sudo docker run --rm -v $PWD:/mnt terhorst/smcpp:latest estimate --spline cubic --knots 15 --timepoints 1000 1000000 --cores 20 -o ./data/estimate-1/ 2e-8 ./data/*.smc.gz
#–spline 線條類(lèi)型 cubic為平滑曲線吧
#–knots 節(jié)點(diǎn)昔馋,就是繪圖曲線的點(diǎn)數(shù)
#–timepoints 時(shí)間范圍筹吐,單位為多少代,如1000 1000000 為1000代到1000000代
#cores 計(jì)算核數(shù)
#2e-8為突變率
4.繪圖
sudo docker run --rm -v $PWD:/mnt terhorst/smcpp:latest plot ./plot-IN-6e-9.pdf ./data/estimate-1/*.final.json -g 1 --ylim 0 100000000 -c
#.final.json為最終模型
#-c --csv 輸出繪圖文件秘遏,利用該文件在R中繪圖丘薛,具體自己發(fā)揮
#-g 代數(shù),如人為25邦危,鼠為1
#–ylim Y軸范圍
#–xlim X軸范圍
總結(jié):smc++不僅能多個(gè)文件分析洋侨,速度還快舍扰,還支持vcf文件,非常好希坚。
stairway plot 分析
————————————————