生信：种群历史有效群体大小推断 SMC++(二)
smc++ plot 绘图效果
cluster文件，两列第一列样本名称第二列group名称
代码块1
cat ../00.raw/cluster | cut -f 2 | sort -u | while read -r a; do
awk -v var="$a" '$2 == var { sp = sp $1 "," } END { print var ":" sp }' ../00.raw/cluster | sed 's/,$//' | awk -v var="$a" '$0 !~ var"$"' > "${a}.list"
done
目标：生成样品信息文件，格式为
pop:sp1,sp2的形式。结果：对于输入文件
../00.raw/cluster，根据第2列的唯一值进行循环。使用awk命令提取符合条件的样品信息，并使用sed命令去掉末尾的逗号，并使用第2列的值作为文件名保存。
代码块2
cat ../00.raw/cluster | cut -f 2 | sort -u | while read -r a; do
awk -v var="$a" '$2 == var {print $1}' ../00.raw/cluster | head -n 4 > "${a}.sample.list"
done
目标：生成Composite likelihood样品列表文件，一般选取2-10个样品。
结果：对于输入文件
../00.raw/cluster，根据第2列的唯一值进行循环。使用awk命令提取符合条件的样品名，并使用head命令限制数量为4个，并保存到以第2列的值为文件名的文件中。
代码块3
bgzip -c /data/Work/22.selective_sweeps/00.raw/all_snps_imputation.eff.reh.chr.vcf > ../00.raw/all_snps_imputation.eff.reh.chr.vcf.gz
tabix ../00.raw/all_snps_imputation.eff.reh.chr.vcf.gz
目标：将 VCF 文件进行压缩，并生成索引文件。
结果：将输入的 VCF 文件
/data/zhougl/Work/22.selective_sweeps/00.raw/all_snps_imputation.eff.reh.chr.vcf使用bgzip命令进行压缩，并将压缩文件保存为../00.raw/all_snps_imputation.eff.reh.chr.vcf.gz。然后使用tabix命令生成该压缩文件的索引。
代码块4
cat ../00.raw/cluster | cut -f 2 | sort -u | while read -r a; do
mkdir "${a}_vcf2smc"
done
目标：创建文件夹以准备进行 VCF 转 SMC 输入格式的操作。
结果：对于输入文件
../00.raw/cluster，根据第2列的唯一值进行循环。使用mkdir命令创建名为${a}_vcf2smc的文件夹。
代码块5
awk 'BEGIN {for (i=1; i<=12; i++) {printf "%d\t%s%02d\n", i, "chr", i}}' | cut -f 2 > chr.list
目标：生成染色体列表文件。
结果：使用
awk命令生成以chr开头的数字编号的染色体列表，并将结果保存到chr.list文件中。
代码块6
cp tmp.sh G2_vcf2smc.sh
cat chr.list | while read chr; do
cat G2.sample.list | while read sample; do
echo "smc++ vcf2smc -m $mask -d $sample $sample $vcf G2_vcf2smc/$sample.$chr.smc.gz $chr `cat G2.list` " >> G2_vcf2smc.sh
done
done
sbatch G2_vcf2smc.sh
目标：执行
smc++ vcf2smc命令将 VCF 文件转换为 SMC 输入格式。结果：首先复制
tmp.sh文件为G2_vcf2smc.sh。然后，对于每个染色体和样品，使用smc++ vcf2smc命令将 VCF 文件转换为 SMC 输入格式，并保存到对应的输出文件中。最后，使用sbatch命令提交G2_vcf2smc.sh脚本进行执行。
代码块7
mkdir G2_analysis
smc++ estimate --cores 10 --knots 10 --spline cubic -o G2_analysis 6.5e-9 G2_vcf2smc/*.smc.gz
目标：进行种群历史的有效群体大小估计。
结果：创建名为
G2_analysis的文件夹。然后，使用smc++ estimate命令进行种群历史的有效群体大小估计，指定核心数为 10，节点数为 10，样条函数为 cubic，并将结果输出到G2_analysis文件夹中，使用突变速率参数为6.5e-9，输入文件为G2_vcf2smc/*.smc.gz。
代码块8
smc++ plot -g 1 --linear G2_analysis.plot.pdf G2_analysis/model.final.json
目标：绘制种群历史估计结果的图表。
结果：使用
smc++ plot命令，指定世代时间比例为 1，绘制线性图表，并将结果输出到G2_analysis.plot.pdf文件中，输入文件为G2_analysis/model.final.json。
国靓的文案
主要分享：1：基于R语言、python和shell的数据分析，数据可视化及生物信息分析；2：零基础学英语从How are you开始学；3：运动健康。（能量是守恒的，喜欢是互相的，关注我，世界上就多了一个爱你的人！）
国靓的文案
做运动 学英语 练瑜伽 听故事 搞科研
公众号
欢迎留言区or后台提问！
点个小赞鼓励我一下子再走呗！
