0 背景：

除了前面笔记所记录的SNP过滤条件外，在张仁纲的毕业论文里发现还有多核苷酸多态性、覆盖深度、TE区域和杂合度这条件需要进行过滤。

前面的过滤笔记链家：

SNP过滤：（1）indel邻近区域的SNP和10bp范围内的SNP cluster

SNP过滤：（2）基因型的质量GQ和缺失率missing

SNP过滤：（3）次等位基因频率MAF

1 过滤变异类型为多核苷酸多态性 (multiple nucleotide polymorphism, MNP) 的位点

# 变异类型为多核苷酸多态性就是REF列和ALT列长度大于1的地方；
# 核对文件里是否有这种类型的行
zcat 3genome_snp_filtered_deleted.vcf.gz | awk 'BEGIN{OFS="\t"} /^##/ {print; next} /^#/ {print; next} length($4)>1 || length($5)>1 {print}' > output_indels.vcf
# 查看结果文件，发现并没有这样的行，因此跳过该步骤

2 过滤覆盖深度高于65的位点

# 注意：过滤这部分的时候，VCF文件必须没有经过LD过滤，否则会丢失DP这种信息
## 统计深度分布
vcftools --vcf 3genome_snp_gap5_filter_DP5_GQ20_missing_maf_SnpCluster.gap5.vcf --site-mean-depth --out my_depth
## 从深度统计文件中，提取大于65的位点个数（一般最大值为测序深度的两倍（严格过滤），但是参考张仁刚毕业论文15x的数据4<DP<34；因此我选择5<DP<65）
awk '$3 > 65' my_depth.ldepth.mean > DP65.txt
## 使用 VCFtools 从原始文件中彻底删除这些位点
vcftools --vcf 3genome_snp_gap5_filter_DP5_GQ20_missing_maf_SnpCluster.gap5.vcf --exclude-positions  DP65.txt --recode --recode-INFO-all --out 4genome_snp_gap5_filter_DP5_GQ20_missing_maf_SnpCluster.gap5.DP65.vcf

3 排除了 TE 区域内的位点

# 1. 从 GFF 注释中提取所有 TE 的区域坐标，存为 BED 文件
## 下载TE注释的gff文件
##核对gff文件 (通常 GFF 第3列的特征名会包含 transposable_element, TE, 或具体家族如 Gypsy, Copia)
cut -f 3 Rmolle_repeat.gff3 | sort -u # 结果发现存在TE_Unknown
## 提取第 1, 4, 5 列
awk '{print $1"\t"$4"\t"$5}' Rmolle_repeat.gff3 > Rmolle_TE_regions.bed
	
# 2. 使用 vcftools 精准“排雷”，剔除落在 TE 区域内的所有 SNP
## 发现bed文件的染色体id 和 vcf文件的id不一致，导致数据过滤失败，因此修改gff文件里面的id，并重新生成bed文件
vcftools --vcf 4genome_snp_gap5_filter_DP5_GQ20_missing_maf_SnpCluster.gap5.DP65.vcf.recode.vcf --exclude-bed Rmolle_TE_regions.bed --recode --recode-INFO-all --out 5genome_snp_gap5_filter_DP5_GQ20_missing_maf_SnpCluster.gap5.DP65.noTE

4 过滤杂合度大于0.7的位点

# 1.计算所有位点的杂合度统计：
vcftools --vcf your_data.vcf --hardy --out hwe_info
# 生成的 hwe_info.hwe 文件中，有一列名为 Obs-Het（观察杂合度），格式通常为 10/30（代表 30 个样本中有 10 个杂合）。

# 2. 利用 awk 筛选出杂合度 > 0.7 的位点位置：
# 假设 Obs-Het 在第 6 列，利用 / 分隔计算比例
awk 'NR>1 { split($3, a, "/"); sum = a[1] + a[2] + a[3];  if (sum > 0 && (a[2]/sum) > 0.7) print $1, $2 }' hwe_info.hwe > high_het_sites.txt  #共有42117行
	 
# 3. 剔除这些位点：
vcftools --vcf your_data.vcf --exclude-positions high_het_sites.txt --recode --out final_filtered

5 过滤最小等位基因计数 (minor allele count, MAC) 大于 2 

# vcftools 参数：--mac 3
vcftools --vcf 6genome_snp_gap5_filter_DP5_GQ20_missing_maf_SnpCluster.gap5.DP65.noTE.het.recode.vcf --mac 3 --recode --recode-INFO-all --out 7genome_snp_gap5_filter_DP5_GQ20_missing_maf_SnpCluster.gap5.DP65.noTE.het.mac3

总结：

经过上述所有过滤步骤后，可以质量较差、测序导致偏移、不可靠等位点删除，保证后续分析的可靠性。