1. 前言
Pggb输出的VCF文件除了有SV还有SNP,因此要把SV单独提取出来,后续所有的分析都是基于SV的。
2. 软件安装
mamba create -n Truvari
mamba activate Truvari
python3 -m pip install Truvari
3. 运行
3.1 鉴定SV类型和长度
truvari anno svinfo -o SV_Chr_1.vcf Chr_1.vcf
3.2 提取SV
bcftools view -i 'SVTYPE="DEL" || SVTYPE="INS"' SV_Chr_1.vcf -o SV_Chr_1_DEL_INS.vcf
bcftools view -i 'SVTYPE="DEL"' SV_Chr_1.vcf -o SV_Chr_1_DEL.vcf
bcftools view -i 'SVTYPE="INS"' SV_Chr_1.vcf -o SV_Chr_1_INS.vcf
提取后其实可以做很多事,比如使用堆叠柱状图展示每条染色体SV的数量。
3.3 统计SV中插入和缺失的数量
grep -v '^#' SV_Chr_1_DEL.vcf|wc -l
grep -v '^#' SV_Chr_1_INS.vcf|wc -l
3.4 统计插入和缺失SV的长度
bcftools view -i 'abs(SVLEN) < 500' SV_Chr_1_DEL_INS.vcf | grep -v '^#'|wc -l
bcftools view -i 'abs(SVLEN) > 500 && abs(SVLEN) < 1000' SV_Chr_1_DEL_INS.vcf | grep -v '^#'|wc -l
bcftools view -i 'abs(SVLEN) > 1000 && abs(SVLEN) < 1500' SV_Chr_1_DEL_INS.vcf | grep -v '^#'|wc -l
bcftools view -i 'abs(SVLEN) > 1500 && abs(SVLEN) < 2000' SV_Chr_1_DEL_INS.vcf | grep -v '^#'|wc -l
bcftools view -i 'abs(SVLEN) > 2000 && abs(SVLEN) < 2500' SV_Chr_1_DEL_INS.vcf | grep -v '^#'|wc -l
bcftools view -i 'abs(SVLEN) > 2500 && abs(SVLEN) < 3000' SV_Chr_1_DEL_INS.vcf | grep -v '^#'|wc -l
bcftools view -i 'abs(SVLEN) > 3000 && abs(SVLEN) < 3500' SV_Chr_1_DEL_INS.vcf | grep -v '^#'|wc -l
bcftools view -i 'abs(SVLEN) > 3500 && abs(SVLEN) < 4000' SV_Chr_1_DEL_INS.vcf | grep -v '^#'|wc -l
bcftools view -i 'abs(SVLEN) > 4000 && abs(SVLEN) < 4500' SV_Chr_1_DEL_INS.vcf | grep -v '^#'|wc -l
bcftools view -i 'abs(SVLEN) > 4500 && abs(SVLEN) < 5000' SV_Chr_1_DEL_INS.vcf | grep -v '^#'|wc -l
bcftools view -i 'abs(SVLEN) > 5000' SV_Chr_1_DEL_INS.vcf | grep -v '^#'|wc -l
3.5 计算Pi值
Vcftools计算Pi的时候必须是严格的二倍体,也就是有两个等位基因,所以要先把单倍体修改为二倍体,也就是添加一个等位基因。
python3 ~/Script/Py/haploids_to_diploids.py SV_Chr_1_DEL_INS.vcf diploids_SV_Chr_1_DEL_INS.vcf
vcftools --vcf diploids_SV_Chr_1_DEL_INS.vcf --window-pi 1000000 --window-pi-step 200000 --out Test
3.6 计算Fst值
python3 ~/Script/Py/haploids_to_diploids.py SV_Chr_1_DEL_INS.vcf diploids_SV_Chr_1_DEL_INS.vcf
vcftools --vcf diploids_SV_Chr_1_DEL_INS.vcf --window-pi 1000000 --window-pi-step 200000 --out Test
其实还可以计算Fst或者其他衡量种群多样性的指标。但是Fst计算的时候只关注于两个种群之间的比较,但是本研究的研究关注的多个种群,所以这里就不往下做了。