你还在为beagle填充烦恼吗？是的，老子头大

最新推荐文章于 2024-03-02 19:51:05 发布

金金金鱼

最新推荐文章于 2024-03-02 19:51:05 发布

阅读量653

点赞数 9

文章标签：大数据人工智能数据库

本文链接：https://blog.csdn.net/jinning3825/article/details/136263855

版权

本文讲述了作者如何通过编写脚本，使用贝格尔工具处理VCF文件，逐步完成数据填充，并详细记录了准确性和基因型一致性等指标的过程。

摘要由CSDN通过智能技术生成

为啥过不去心理这坎，哦因为老子不想简单一步步完成。

第一天列计划，第二天开始做，小小beagle能难倒老子吗，一顿操作猛如虎，原来是能的。。。

老子非要写成一个脚本代码，这个脚本也就改了几十遍吧！

#!/bin/bash
#SBATCH -p v6_384
#SBATCH -N 1
#SBATCH -n 60
#SBATCH -c 1

# 输入文件
original_vcf="/public5/home/sch7739/WORKPLACE/beagle/tian/DS_AW_ex_zi.vcf.gz"
reference_vcf="/public5/home/sch7739/WORKPLACE/beagle/tian/AW14_xinpian_tian.vcf.gz"
# 输出文件前缀
output_prefix1="fill_data"
output_prefix3="fill_data1"
output_prefix2="filled_data"

# 初始样本集合
initial_sample_set="initial_sample_set.txt"
bcftools query -l "$original_vcf" > "$initial_sample_set"
sp=$(bcftools query -l sites_gt.vcf.gz)
# 设置参数
step_size=100
max_samples=$(wc -l < "$initial_sample_set")

# 创建结果文件
result_file="accuracy_results.txt"
echo "Iteration Sample_Reduction R_Squared Genotype_Concordance Allele_Concordance" > "$result_file"

# 函数：减少样本集合
reduce_sample_set() {
    local input_set="$1"
    local output_set="$2"
    local reduction="$3"
    local total_samples=$(cat "$input_set" | wc -l)
    local reduced_samples=$((total_samples - reduction))

    # 随机选择一定数量的样本
   shuf "$input_set" | head -n "$reduced_samples" > "$output_set"
 }

# 函数：计算填充准确性
calculate_accuracy() {
    local filled="$1"
    local filled1="$2"
    local filled2="$3"
    #就是填充的准确性，DR2，因为全部的不高就取了排行前1000k的位点的DR
    local r_squared=$(zcat "$filled" | awk '$1!~/#/ {print $8}' | sed 's/[=;]/ /g' | awk '{print $2}'|sort -nr | head -n 1000000 | awk '{ sum += $1 } END { if (NR > 0) print sum / NR }')
    #这个就是全部的
    local r_squared_all=$(zcat "$filled1" | awk '$1!~/#/ {print $8}' | sed 's/[=;]/ /g' | awk '{print $2}'|sort -nr | awk '{ sum += $1 } END { if (NR > 0) print sum / NR }')
    #local r_squared=$(bcftools stats -S /public5/home/sch7739/WORKPLACE/beagle/tian/AW14.txt /public5/home/sch7739/WORKPLACE/beagle/tian/DS_AW14_r.vcf.gz "$filled" | grep "R-squared" | awk '{print $3}')
   #这个是实际的基因型和填充后的基因型的的一致性，这个是每个vcf的每个样本取平均
   local genotype_concordance=$(awk '{sum += $2; count++} END {print sum/count}' "$filled2")
    #local allele_concordance=$(bcftools gtcheck -g - -r /public5/home/sch7739/WORKPLACE/beagle/tian/DS_AW14_r.vcf.gz "$filled" | grep "concordant" | awk '{print $3}')
    echo "$r_squared $r_squared_all $genotype_concordance"
    #echo "$r_squared $genotype_concordance $allele_concordance"
}

# 迭代过程
for ((i = 1; i <= max_samples; i += step_size)); do
    # 减少样本集合
   reduced_sample_set="reduced_sample_set_${i}.txt"
    reduce_sample_set "$initial_sample_set" "$reduced_sample_set" $i

    # 使用 Beagle 或其他填充工具填充数据（替换为实际的填充命令）
    fill_vcf="${output_prefix1}_${i}.vcf.gz"
    fill_vcf1="${output_prefix3}_${i}.vcf.gz"
    filled_vcf="${output_prefix2}_${i}.vcf.gz"
    filled_vcf1="${output_prefix2}_${i}.vcf.gz.vcf.gz"
    filled_vcf2="${output_prefix2}1_${i}.vcf.gz"
    filled_vcf3="${output_prefix2}100k_${i}.vcf.gz"
    filled_vcf4="AW_100k_${i}.vcf.gz"
    #先把减少的个体选出来，$original_vcf一定是自填充后的，下面几步注释掉的就是没自填充走过的湾路，这弯路十八湾 ，像极了隔壁实验室的大黄，都是湾的four
     bcftools view -S "$reduced_sample_set" $original_vcf -Oz -o $fill_vcf
     bcftools index $fill_vcf
    # bcftools annotate --columns +FORMAT/GT $fill_vcf | awk 'BEGIN {OFS="\t"} {for(i=10;i<=NF;i++) {if($i==".") $i="./."} }1' | bcftools view -Oz -o $fill_vcf1
    #vcftools --gzvcf $fill_vcf  --recode --recode-INFO-all --stdout --max-missing 1 | sed 's:/:|:g' | gzip > $fill_vcf1
    #java -Xmx240g -jar ~/software/beagle/bref3.22Jul22.46e.jar $fill_vcf > low_${i}.bref3
   #开始填就行了，上面几步就是俺走过的坑，乖，别管
    java -Xmx240g -jar ~/anaconda3/envs/GWAS/share/beagle-5.4_22Jul22.46e-0/beagle.jar gt="$reference_vcf" ref="$fill_vcf" out="$filled_vcf" nthreads=60
    bcftools index $filled_vcf
   #bcftools view -S /public5/home/sch7739/WORKPLACE/beagle/tian/AW14.txt -Oz -o $filled_vcf2 $filled_vcf1 #这个也是坑别管，浪费俺的时间
   #这个是俺取了DR的最小值，为了取前1000k的那个点，其实我就想说，你看结果也不是不好，填到1000准确性也挺好的，
   min_value=$(zcat $filled_vcf | awk '$1!~/#/ {print $8}' | sed 's/[=;]/ /g' | awk '{print $2}'|sort -nr | head -n 1000000 |awk 'NR == 1 || $1 < min { min = $1 } END { print min }')
   bcftools view -i "INFO/DR2>=${min_value}" $filled_vcf -Oz -o $filled_vcf3
   bcftools index $filled_vcf3
   # $filled_vcf3这个就是俺填充前1000k的点,下面把这1000k的点从真实测序中提出来，到要看看，填充前后到底填了什么玩意
   bcftools isec -n=2 -w1 -Oz -o $filled_vcf4 /public5/home/sch7739/WORKPLACE/beagle/tian/DS_AW14_r.vcf.gz $filled_vcf3
    bcftools index $filled_vcf4
#这个是为了那个一致性真实的和填充后的每个样本做了对比，这个可以拿到每个样本填充的一致性也是百分比。
 for sample in $sp; do
   bcftools gtcheck -g sites_gt.vcf.gz filled_data_1.vcf.gz -s "$sample" |grep "$sample" | awk '{print $2,$4}' | tail -1 | awk '{result = ($2 - $1) / $2; print result}' >> gt_${i}.txt
 done
   # 计算填充准确性指标
    accuracy_metrics=$(calculate_accuracy "$filled_vcf" "$filled_vcf" "gt_${i}.txt") 

    # 记录结果到文件
    echo "$i $i $accuracy_metrics" >> "$result_file"
done

不想写注释。。。。

不写估计我自己都不记得了！我要写！

写好了，注释一写，感觉自己真特么拉，写的都是什么玩意....

嘤嘤嘤嘤