smc++又称为smcpp,
从官方github下载文件
https://github.com/popgenmethods/smcpp
git clone https://github.com/popgenmethods/smcpp.git
注意不要下载右侧发布的版本，而是下载项目，这才是最新版。
下载后进入路径，
make -j6
make install
安装完成后，运行smc++ --help
没有报错，但是运行第一步时，报错找不到libgsl.so.0文件
手动从https://mirror.ibcp.fr/pub/gnu/gsl/gsl-latest.tar.gz下载gsl,
在本地执行编译，
./configure --prefix=/share/home/chai/soft/gsl
make -j6
make install
把下面三行文本添加到~/.bashrc文件里。
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/share/home/chai/soft/gsl/lib
export CFLAGS="-I/share/home/chai/soft/gsl/include"
export LDFLAGS="-L/share/home/chai/soft/gsl/lib"
重新进入终端后，查看/share/home/chai/soft/gsl/lib目录里是否有libgsl.so.0文件
find /share/home/chai/soft/gsl/lib -name "libgsl.so.0"
#如果没有上述文件，则运行下面的命令即可
cd /share/home/chai/soft/gsl/lib
ln -s libgsl.so libgsl.so.0
运行smc++ estimate命令时，出现错误“AttributeError: module 'typing' has no attribute '_ClassVar'”
解决方法是：pip uninstall dataclasses -y运行之后，就恢复正常。
运行smc++包括下面的步骤：
注意：输入的vcf是需要使用bgzip压缩的gz文件
getvcfchr.py 的脚本的内容如下：
#!/usr/bin/env python3
#输入 vcf文件和输出文件名
#输出的是vcf的chr的列表
import sys
import gzip
def get_unique_chromosomes_from_vcf(filename):
chromosomes = set()
opener = gzip.open if filename.endswith('.gz') else open
with opener(filename, 'rt') as file:
for line in file:
if line.startswith('#'):
continue
else:
chromosome = line.split('\t')[0]
chromosomes.add(chromosome)
return sorted(chromosomes)
def write_chromosomes_to_file(chromosomes, output_file):
with open(output_file, "w") as file:
for chromosome in chromosomes:
file.write(f"{chromosome}\n")
if __name__ == "__main__":
if len(sys.argv) != 3:
print("python3 getvcfchr.py inputvcf.gz outputtxt")
else:
vcf_filename = sys.argv[1]
output_filename = sys.argv[2]
unique_chromosomes = get_unique_chromosomes_from_vcf(vcf_filename)
write_chromosomes_to_file(unique_chromosomes, output_filename)
print("Chromosomes written to file:", output_filename)
输入文件
input.vcf (snp位点的vcf文件,vcf文件必须有头部，即包含了染色体名称和长度信息的行)
输入的vcf的文件的头部必须包含类似下面的内容的行，smcpp需要从中读取染色体名称和长度，没有这些内容会报错。
##contig=<ID=Chr01,assembly=ref.fa,length=12346530>Groupinfo.txt （样本的分群信息）
准备工作
inputvcf="youinput.snp.vcf"
#获取chr的列表
python3 getvcfchr.py ${inputvcf} chr.list
#使用bgzip压缩为gz格式
bgzip ${inputvcf} -c >${inputvcf}.gz
#构建索引
tabix ${inputvcf}.gz
开始第一步分析，获取每条染色体的smc.gz文件
这里的chr.list就是一列染色体的编号的文件，就是vcf文件里的染色体号。
此处示例是样本分为2个群，G1是第一个群的名称,后面的sample11，sample12,sample41,sampleM是属于这个群的样本的名称（一定要和vcf文件里的一致）
如果有多个群，也是一样的方法，依次增加即可。注意输出的文件夹G1这种，一定要自己手动新建。如果目录不存在输出的文件夹，smc++不会自动新建。
mkdir G1
mkdir G2
cat chr.list|while read line;
do
smc++ vcf2smc ${inputvcf}.gz ./G1/$line.smc.gz $line G1:sample11,sample12,sample41,sampleM
done
cat chr.list|while read line;
do
smc++ vcf2smc ${inputvcf}.gz ./G2/$line.smc.gz $line G2:sample1,sample2,sample4,sampleN
done
运行后会输出结果分别在G1和G2两个文件夹内，如果你的染色体名称是chr1这种格式，G1内会有chr1.smc.gz,chr2.smc.gz这种文件.G2内也是一样。
第2步 使用smc++分析计算（需要使用较多的cpu数量）
普通版
smc++ estimate -o G1 6.9e-09 G1/*.smc.gz --cores 48
默认情况下只需要
-o 指定输出目录
mu 每一代每一碱基对的突变频率
--cores 指定使用的cpu数量
进阶版
smc++ estimate --spline cubic --knots 15 --timepoints 1000 1000000 --cores 48 -o G1 6.9e-09 G1/*.smc.gz
–spline 线条类型 cubic为平滑曲线,
–knots 节点，就是绘图曲线的点数
–timepoints 时间范围，单位为多少代，如1000 1000000 为1000代到1000000代
第2步运行完成后，会在G1目录输出一个文件model.final.json.
运行的时候可能会报错RuntimeError: erroneous average coalescence time.出现这种情况可能是因为你的group内的遗传距离太近了，修改你的mu值，可能会解决这个问题。
第3步，可视化
普通版
smc++ plot plot.G1.pdf G1/*.final.json
输出文件为plot.G1.pdf
进阶版
smc++ plot plot.pdf *.final.json -g 6 --ylim 0 100000000 -c
-c --csv 输出绘图文件，利用该文件在R中绘图（也许可以画折线图）
-g 代数，如人为25，鼠为1 ，牛为6或5
–ylim Y轴范围
–xlim X轴范围
合并计算亚群的分化时间（一次只能比较2个群之间的信息）
mkdir data
cd data
ls ../G1/*.smc.gz|while read line;do ln -s ${line} G1$(basename ${line}); done
ls ../G2/*.smc.gz|while read line;do ln -s ${line} G2$(basename ${line}); done
cd ../
smc++ split -o split/ G1/model.final.json G2/model.final.json data/*.smc.gz
smc++ plot G1.G2.pdf split/model.final.json
一次绘制3个或多个群的信息
#加入分为3个群，分别单独完成运算，生成model.final.json文件
cp G1/model.final.json G1.model.final.json
cp G2/model.final.json G2.model.final.json
cp G3/model.final.json G3.model.final.json
smc++ plot all.pdf *model.final.json
这样就可以一次绘制多个的结果到一张图上
msmc2的教程
msmc-tools/msmc-tutorial/guide.md ·STSCHIFF/MSMC-tools ·GitHub
可以提供基因组的mask文件
ref_genome="~/Gh_meta/ref.genome.fa"
# https://lh3lh3.users.sourceforge.net/snpable.shtml
~/soft/msmc2/seqbility-20091110/splitfa ${ref_genome} 35 > ref.fq
bwa index -a bwtsw ${ref_genome}
bwa aln -R 1000000 -O 3 -E 3 -t 32 ${ref_genome} ref.fq | bwa samse ${ref_genome} - ref.fq > ref.sam
~/soft/msmc2/seqbility-20091110/gen_raw_mask.pl ref.sam > ref.sam.mask
~/soft/msmc2/seqbility-20091110/gen_mask -l 35 -r 0.5 ref.sam.mask > ref.sam.mask.fa
#makeMappabilityMask.py这个是python2的脚本，原始是msmc2提供的，此处是我修改了支持自定义输入和输出文件
makeMappabilityMask.py ref.sam.mask.fa ref.{}.mask.bed.gz
rm -rf ref.sam
rm -rf ref.fq
rm -rf ref.sam.mask
#输出的mask结果文件是
#ref.{}.mask.bed.gz
以上的分析步骤已经打包成流程runsmcpp,有问题直接提issue。
