教程：https://mp.weixin.qq.com/s/JdpUYgPgnL3EKhVSEQ0_RA

01 安装SMC++

conda create -n smcpp_env
conda activate smcpp_env
conda install -c bioconda smcpp

02 准备 mask 文件标记缺失数据区域

（1）数据准备

# 基因组文件
# vcf文件

（2）mask 文件生成

# 生成50bp模拟reads
splitfa genome.fa 50 > genome_50bp_reads.fa
# 构建bwa索引(之前构建过，跳过)
bwa index genome.fa
# 比对
bwa aln -t 8 -R 1000000 -O 3 -E 3 genome.fa genome_50bp_reads.fa > genome_50bp.aln.sai
# 转SAM格式
bwa samse genome.fa genome_50bp.aln.sai genome_50bp_reads.fa > genome_50bp.sam
# 生成原始mask
gen_raw_mask.pl genome_50bp.sam > genome_raw_mask.txt
# 过滤mask（阈值r=0.5）
gen_mask -l 50 -r 0.5 genome_raw_mask.txt > mask_50.fa
# 生成masked基因组
apply_mask_s mask_50.fa genome.fa > genome.mask.fa
# 提取mask区域为bed（压缩+索引）
perl get_lcase.pl genome.mask.fa > genome.mask.bed
bgzip genome.mask.bed
tabix genome.mask.bed.gz

03 运行 smc++

（1）文件准备

# （1）准备样品信息，生成pop:sp1,sp2格式文件
[root@3eb13a5fb4e6 05SMC++]$head pop.cov.new
AH-1 SAD
AH-2 SAD
BZW-1 XG 
BZW-2 XG
## 生成样品列表文件
cat pop.cov.new| awk '$2=="RB" { sp = sp $1 ","} END{print "RB:" sp } ' | sed 's/,$//' > RB.list

#（2）准备Composite likelihood样品列表，一般选2-10个，选3个
cat pop.cov.new | awk '$2=="RB"{print $1}' | head -n 3 > RB_sample.list

#（3）准备vcf文件(SMC++格式要求)
## 需要建立索引tbi

（2） vcf转smc输入格式

###### vcf转smc输入格式
## 批量生成vcf转smc格式脚本列表
## 需要分样品分染色体运行
cat chr.list |while read chr
do
cat RB_sample.list | while read sample
do
echo "smc++ vcf2smc -m Rmolle_genomic_GCA_025413875.1.mask.bed.gz -d $sample $sample ../re6genome_snp_gap5_filter_DP5_GQ20_missing_SnpCluster_3filter_gap5_DP65_noTE_het.recode.vcf.gz RB_vcf2smc/$sample.$chr.smc.gz $chr `cat RB.list`"
done
done > RB_vcf2smc.sh
# 运行
$ sh RB_vcf2smc.sh

 (3) 进行种群历史有效群体大小估计

###### 进行种群历史有效群体大小估计
## 进行估计分析，生成结果为Msie_analysis/model.final.json
$ smc++ estimate \
--cores 10 \ # 线程数
--knots 10 \ # 样条节点数，越小月平滑
--spline piecewise \ # 默认是piecewise，分段结果；cubic 结果为平滑线
-o RB_analysis \ # 输出结果目录
4e-9 \ # 每代突变速率
RB_vcf2smc/*.smc.gz # 输入smc文件


## 结果绘图
singularity exec ../../software/smcpp_latest.sif smc++ plot \
-g 5 \ # 世代间隔
--linear \ # 纵坐标不取对数
RB.plot.pdf \ # 输出图片名称
RB_analysis/model.final.json # 输入json文件


#### 所有亚群结果一起绘图
smc++ plot -g 5 --linear all.plot.pdf RB_analysis/model.final.json NWG_analysis/model.final.json EFJ_analysis/model.final.json CFJ_analysis/model.final.json CZJ_analysis/model.final.json AJB_analysis/model.final.json DML_analysis/model.final.json NJS_analysis/model.final.json SAD_analysis/model.final.json XG_analysis/model.final.json HN_analysis/model.final.json