Hello, I am trying to compare the quality of the samples within my vcf file. I want to compare the quality of sequencing across two metrics, the different subspecies in my file and the sample type (modern vs. historical). I have been using this script, but when I run it, the output files say that the mean and standard deviation for quality is the same across all the subspecies and in both types of samples which I know is not the case. In the script, I start with a vcf file that contains all the samples, and then filter it based on lists of samples within each subspecies and within each sample type. I am trying to look at the QUAL scores and the QD scores (but if anyone has suggestions for other metrics to look at please let me know).
# Define the directory and file paths
VCF_DIR="/project/sedmands_1143/ekirsch/savannahsparrow/savannahsparrow/41-Passerculus/41-Passerculus/data"
ALL_SAMPLES_VCF="lineremoved_GATK_bgzip_filtered_passerculus.vcf.gz" # Your comprehensive VCF file
SUBSPECIES_DIR="/project/sedmands_1143/ekirsch/savannahsparrow/savannahsparrow/41-Passerculus/41-Passerculus"
TIME_PERIOD_DIR="/project/sedmands_1143/ekirsch/savannahsparrow/savannahsparrow/41-Passerculus/41-Passerculus"
OUTPUT_DIR="/project/sedmands_1143/ekirsch/savannahsparrow/savannahsparrow/41-Passerculus/41-Passerculus/output"
# Subspecies list files
SUBSPECIES=("alaudinus" "brooksi" "beldingi" "rostratus" "nevadensis")
# Load bcftools (modify this based on your environment)
module load bcftools
# Function to filter the VCF file by subspecies and time period
filter_vcf_by_subspecies_and_time() {
local subspecies_list=$1
local time_period_list=$2
local vcf_file=$3
local output_vcf=$4
bcftools view -S $subspecies_list -S $time_period_list $vcf_file -Oz -o $output_vcf
}
# Function to calculate quality stats for QUAL and QD
calculate_quality_stats() {
local vcf_file=$1
local output_file=$2
bcftools query -f '%QUAL\t%QD\n' $vcf_file > $output_file
awk '
BEGIN { sum_qual=0; sumsq_qual=0; sum_qd=0; sumsq_qd=0; count=0 }
{
sum_qual += $1;
sumsq_qual += ($1)^2;
sum_qd += $2;
sumsq_qd += ($2)^2;
count++;
}
END {
if (count > 0) {
mean_qual = sum_qual / count;
variance_qual = (sumsq_qual - (sum_qual^2)/count) / (count - 1);
stddev_qual = sqrt(variance_qual);
mean_qd = sum_qd / count;
variance_qd = (sumsq_qd - (sum_qd^2)/count) / (count - 1);
stddev_qd = sqrt(variance_qd);
print "Mean QUAL:", mean_qual;
print "Standard Deviation QUAL:", stddev_qual;
print "Mean QD:", mean_qd;
print "Standard Deviation QD:", stddev_qd;
} else {
print "No data available.";
}
}
' $output_file
}
# Loop through each subspecies and filter by time period, then calculate quality scores
for subspecies in "${SUBSPECIES[@]}"; do
subspecies_list="${SUBSPECIES_DIR}/${subspecies}_list.txt"
# Process historical samples
historical_output_vcf="${OUTPUT_DIR}/${subspecies}_historical.vcf.gz"
filter_vcf_by_subspecies_and_time $subspecies_list "${TIME_PERIOD_DIR}/historical_samples.txt" "$VCF_DIR/$ALL_SAMPLES_VCF" $historical_output_vcf
calculate_quality_stats $historical_output_vcf "${OUTPUT_DIR}/quality_scores_${subspecies}_historical.txt" > "${OUTPUT_DIR}/stats_${subspecies}_historical_quality.txt"
# Process modern samples
modern_output_vcf="${OUTPUT_DIR}/${subspecies}_modern.vcf.gz"
filter_vcf_by_subspecies_and_time $subspecies_list "${TIME_PERIOD_DIR}/modern_samples.txt" "$VCF_DIR/$ALL_SAMPLES_VCF" $modern_output_vcf
calculate_quality_stats $modern_output_vcf "${OUTPUT_DIR}/quality_scores_${subspecies}_modern.txt" > "${OUTPUT_DIR}/stats_${subspecies}_modern_quality.txt"
done
echo "Quality score calculation completed."
If anyone can find an error in my script or suggest another way to compare quality scores that would be amazing. Thank you!
