Step 6: Evaluation Against the Truth Set - aduckwo4-boop/EPP622-2026-SNPCalling GitHub Wiki

Step 6: Evaluation Against the Truth Set

Commands

Link to the provided evaluate_snps_v4.py and truth_snps.vcf. Find the precision, recall and F1 for gently filtered and stringently filtered SNP sets.

Example of how to run evaluate_snps_v4.py:

python3 evaluate_snps_v4.py \
    --truth truth_snps.vcf --calls <yourvcffile> \
    --all-samples --tag strict --outdir outdir

The raw, gently filtered, and stringently filtered SNP calls were compared against the provided truth SNP dataset. The evaluation calculated true positives (TP), false positives (FP), false negatives (FN), precision, recall, and F1 scores for each sample to assess how filtering affected SNP-calling accuracy. Created the evaluation directory and linked the required files:

mkdir -p 06_evaluate_snps
cd 06_evaluate_snps
ln -s ../check/evaluate_snps_v4.py .
ln -s ../check/truth_snps.vcf .
ln -s ../04_bcftools/all_samples.gentle.vcf.gz .
ln -s ../04_bcftools/all_samples.stringent.vcf.gz .
ln -s ../04_bcftools/all_samples.raw.vcf.gz .



# mkdir -p: creates the evaluation directory if it does not exist.
# cd: changes to the evaluation directory.
# ln -s: creates symbolic links without copying the original files.
# evaluate_snps_v4.py: Python script used to evaluate SNP calls.
# truth_snps.vcf: reference dataset containing known SNPs.
# all_samples.gentle.vcf.gz: gently filtered SNP calls.
# all_samples.stringent.vcf.gz: stringently filtered SNP calls.
# all_samples.raw.vcf.gz: unfiltered SNP calls.
# .: places the symbolic links in the current directory.

Ran the evaluation for all three filtering approaches:

#Gentle filtering
python3 evaluate_snps_v4.py \
  --truth truth_snps.vcf \
  --calls all_samples.gentle.vcf.gz \
  --all-samples \
  --tag gentle \
  --outdir gentle_results

#Stringent Filtering
python3 evaluate_snps_v4.py \
  --truth truth_snps.vcf \
  --calls all_samples.stringent.vcf.gz \
  --all-samples \
  --tag stringent \
  --outdir stringent_results

#Raw SNP Calls
python3 evaluate_snps_v4.py \
  --truth truth_snps.vcf \
  --calls all_samples.raw.vcf.gz \
  --all-samples \
  --tag raw \
  --outdir raw_results


# python3: executes the SNP evaluation Python script.
# --truth: specifies the VCF containing the known true SNPs.
# --calls: specifies the VCF containing the SNP calls being evaluated.
# --all-samples: evaluates each sample in the multisample VCF.
# --tag: assigns a label identifying the filtering approach.
# --outdir: specifies the directory for evaluation results.

Directory structure

From this folder that holds all my analysis, here is the directory structure (paste the output of: tree .):

#output of tree command:
.
├── 01_fastqc
│   ├── fastqc.bash
│   ├── multiqc_data
│   │   ├── fastqc_adapter_content_plot.txt
│   │   ├── fastqc_overrepresented_sequences_plot.txt
│   │   ├── fastqc_per_base_n_content_plot.txt
│   │   ├── fastqc_per_base_sequence_quality_plot.txt
│   │   ├── fastqc_per_sequence_gc_content_plot_Counts.txt
│   │   ├── fastqc_per_sequence_gc_content_plot_Percentages.txt
│   │   ├── fastqc_per_sequence_quality_scores_plot.txt
│   │   ├── fastqc_sequence_counts_plot.txt
│   │   ├── fastqc_sequence_duplication_levels_plot.txt
│   │   ├── fastqc-status-check-heatmap.txt
│   │   ├── fastqc_top_overrepresented_sequences_table.txt
│   │   ├── llms-full.txt
│   │   ├── multiqc_citations.txt
│   │   ├── multiqc_data.json
│   │   ├── multiqc_fastqc.txt
│   │   ├── multiqc_general_stats.txt
│   │   ├── multiqc.log
│   │   ├── multiqc.parquet
│   │   ├── multiqc_software_versions.txt
│   │   └── multiqc_sources.txt
│   ├── multiqc_report.html
│   ├── Sample_A_R1_fastqc.html
│   ├── Sample_A_R1_fastqc.zip
│   ├── Sample_A_R1.fastq.gz -> ../reads/Sample_A_R1.fastq.gz
│   ├── Sample_A_R2_fastqc.html
│   ├── Sample_A_R2_fastqc.zip
│   ├── Sample_A_R2.fastq.gz -> ../reads/Sample_A_R2.fastq.gz
│   ├── Sample_B_R1_fastqc.html
│   ├── Sample_B_R1_fastqc.zip
│   ├── Sample_B_R1.fastq.gz -> ../reads/Sample_B_R1.fastq.gz
│   ├── Sample_B_R2_fastqc.html
│   ├── Sample_B_R2_fastqc.zip
│   ├── Sample_B_R2.fastq.gz -> ../reads/Sample_B_R2.fastq.gz
│   ├── Sample_C_R1_fastqc.html
│   ├── Sample_C_R1_fastqc.zip
│   ├── Sample_C_R1.fastq.gz -> ../reads/Sample_C_R1.fastq.gz
│   ├── Sample_C_R2_fastqc.html
│   ├── Sample_C_R2_fastqc.zip
│   ├── Sample_C_R2.fastq.gz -> ../reads/Sample_C_R2.fastq.gz
│   ├── Sample_D_R1_fastqc.html
│   ├── Sample_D_R1_fastqc.zip
│   ├── Sample_D_R1.fastq.gz -> ../reads/Sample_D_R1.fastq.gz
│   ├── Sample_D_R2_fastqc.html
│   ├── Sample_D_R2_fastqc.zip
│   └── Sample_D_R2.fastq.gz -> ../reads/Sample_D_R2.fastq.gz
├── 02_fastp
│   ├── fastp_trimmed
│   │   ├── Sample_A.fastp.html
│   │   ├── Sample_A.fastp.json
│   │   ├── Sample_A_R1.trimmed.fastq.gz
│   │   ├── Sample_A_R2.trimmed.fastq.gz
│   │   ├── Sample_B.fastp.html
│   │   ├── Sample_B.fastp.json
│   │   ├── Sample_B_R1.trimmed.fastq.gz
│   │   ├── Sample_B_R2.trimmed.fastq.gz
│   │   ├── Sample_C.fastp.html
│   │   ├── Sample_C.fastp.json
│   │   ├── Sample_C_R1.trimmed.fastq.gz
│   │   ├── Sample_C_R2.trimmed.fastq.gz
│   │   ├── Sample_D.fastp.html
│   │   ├── Sample_D.fastp.json
│   │   ├── Sample_D_R1.trimmed.fastq.gz
│   │   └── Sample_D_R2.trimmed.fastq.gz
│   ├── run_fastp.sh
│   ├── Sample_A_R1.fastq.gz -> ../reads/Sample_A_R1.fastq.gz
│   ├── Sample_A_R2.fastq.gz -> ../reads/Sample_A_R2.fastq.gz
│   ├── Sample_B_R1.fastq.gz -> ../reads/Sample_B_R1.fastq.gz
│   ├── Sample_B_R2.fastq.gz -> ../reads/Sample_B_R2.fastq.gz
│   ├── Sample_C_R1.fastq.gz -> ../reads/Sample_C_R1.fastq.gz
│   ├── Sample_C_R2.fastq.gz -> ../reads/Sample_C_R2.fastq.gz
│   ├── Sample_D_R1.fastq.gz -> ../reads/Sample_D_R1.fastq.gz
│   ├── Sample_D_R2.fastq.gz -> ../reads/Sample_D_R2.fastq.gz
│   └── trimmed_fastqc
│       ├── multiqc_trimmed_report_data
│       │   ├── fastqc_adapter_content_plot.txt
│       │   ├── fastqc_per_base_n_content_plot.txt
│       │   ├── fastqc_per_base_sequence_quality_plot.txt
│       │   ├── fastqc_per_sequence_gc_content_plot_Counts.txt
│       │   ├── fastqc_per_sequence_gc_content_plot_Percentages.txt
│       │   ├── fastqc_per_sequence_quality_scores_plot.txt
│       │   ├── fastqc_sequence_counts_plot.txt
│       │   ├── fastqc_sequence_duplication_levels_plot.txt
│       │   ├── fastqc_sequence_length_distribution_plot.txt
│       │   ├── fastqc-status-check-heatmap.txt
│       │   ├── llms-full.txt
│       │   ├── multiqc_citations.txt
│       │   ├── multiqc_data.json
│       │   ├── multiqc_fastqc.txt
│       │   ├── multiqc_general_stats.txt
│       │   ├── multiqc.log
│       │   ├── multiqc.parquet
│       │   ├── multiqc_software_versions.txt
│       │   └── multiqc_sources.txt
│       ├── multiqc_trimmed_report.html
│       ├── Sample_A_R1.trimmed_fastqc.html
│       ├── Sample_A_R1.trimmed_fastqc.zip
│       ├── Sample_A_R2.trimmed_fastqc.html
│       ├── Sample_A_R2.trimmed_fastqc.zip
│       ├── Sample_B_R1.trimmed_fastqc.html
│       ├── Sample_B_R1.trimmed_fastqc.zip
│       ├── Sample_B_R2.trimmed_fastqc.html
│       ├── Sample_B_R2.trimmed_fastqc.zip
│       ├── Sample_C_R1.trimmed_fastqc.html
│       ├── Sample_C_R1.trimmed_fastqc.zip
│       ├── Sample_C_R2.trimmed_fastqc.html
│       ├── Sample_C_R2.trimmed_fastqc.zip
│       ├── Sample_D_R1.trimmed_fastqc.html
│       ├── Sample_D_R1.trimmed_fastqc.zip
│       ├── Sample_D_R2.trimmed_fastqc.html
│       └── Sample_D_R2.trimmed_fastqc.zip
├── 03_bwamem2
│   ├── multiqc_bam_report.html
│   ├── multiqc_data
│   │   ├── llms-full.txt
│   │   ├── multiqc_citations.txt
│   │   ├── multiqc_data.json
│   │   ├── multiqc_general_stats.txt
│   │   ├── multiqc.log
│   │   ├── multiqc.parquet
│   │   ├── multiqc_samtools_coverage.txt
│   │   ├── multiqc_samtools_flagstat.txt
│   │   ├── multiqc_samtools_stats.txt
│   │   ├── multiqc_software_versions.txt
│   │   ├── multiqc_sources.txt
│   │   ├── samtools_alignment_plot.txt
│   │   ├── samtools-coverage_Bases.txt
│   │   ├── samtools-coverage_BQ.txt
│   │   ├── samtools-coverage_Coverage.txt
│   │   ├── samtools-coverage_Mean_depth.txt
│   │   ├── samtools-coverage_MQ.txt
│   │   ├── samtools-coverage_Reads.txt
│   │   ├── samtools-coverage-table.txt
│   │   ├── samtools-flagstat-pct-table.txt
│   │   ├── samtools-flagstat-table.txt
│   │   ├── samtools_insert_size.txt
│   │   └── samtools-stats-dp.txt
│   ├── multiqc_report.html
│   ├── ref -> /lustre/isaac24/proj/UTK0505/aduckwo4/622_test2/ref
│   ├── run_bwa.sh
│   ├── run_sam.sh
│   ├── Sample_A.coverage.txt
│   ├── Sample_A.flagstat.txt
│   ├── Sample_A_R1.trimmed.fastq.gz -> ../02_fastp/fastp_trimmed/Sample_A_R1.trimmed.fastq.gz
│   ├── Sample_A_R2.trimmed.fastq.gz -> ../02_fastp/fastp_trimmed/Sample_A_R2.trimmed.fastq.gz
│   ├── Sample_A.sorted.bam
│   ├── Sample_A.sorted.bam.bai
│   ├── Sample_A.stats.txt
│   ├── Sample_B.coverage.txt
│   ├── Sample_B.flagstat.txt
│   ├── Sample_B_R1.trimmed.fastq.gz -> ../02_fastp/fastp_trimmed/Sample_B_R1.trimmed.fastq.gz
│   ├── Sample_B_R2.trimmed.fastq.gz -> ../02_fastp/fastp_trimmed/Sample_B_R2.trimmed.fastq.gz
│   ├── Sample_B.sorted.bam
│   ├── Sample_B.sorted.bam.bai
│   ├── Sample_B.stats.txt
│   ├── Sample_C.coverage.txt
│   ├── Sample_C.flagstat.txt
│   ├── Sample_C_R1.trimmed.fastq.gz -> ../02_fastp/fastp_trimmed/Sample_C_R1.trimmed.fastq.gz
│   ├── Sample_C_R2.trimmed.fastq.gz -> ../02_fastp/fastp_trimmed/Sample_C_R2.trimmed.fastq.gz
│   ├── Sample_C.sorted.bam
│   ├── Sample_C.sorted.bam.bai
│   ├── Sample_C.stats.txt
│   ├── Sample_D.coverage.txt
│   ├── Sample_D.flagstat.txt
│   ├── Sample_D_R1.trimmed.fastq.gz -> ../02_fastp/fastp_trimmed/Sample_D_R1.trimmed.fastq.gz
│   ├── Sample_D_R2.trimmed.fastq.gz -> ../02_fastp/fastp_trimmed/Sample_D_R2.trimmed.fastq.gz
│   ├── Sample_D.sorted.bam
│   ├── Sample_D.sorted.bam.bai
│   └── Sample_D.stats.txt
├── 04_bcftools
│   ├── all_samples.filtered.stats.txt
│   ├── all_samples.filtered.vcf.gz
│   ├── all_samples.gentle.stats.txt
│   ├── all_samples.gentle.vcf.gz
│   ├── all_samples.raw.fixed.vcf.gz
│   ├── all_samples.raw.stats.txt
│   ├── all_samples.raw.vcf.gz
│   ├── all_samples.stringent.stats.txt
│   ├── all_samples.stringent.vcf.gz
│   ├── bam_list.txt
│   ├── header.txt
│   ├── ref -> ../ref
│   └── run_bcftools.sh
├── 06_evaluate_snps
│   ├── all_samples.gentle.vcf.gz -> ../04_bcftools/all_samples.gentle.vcf.gz
│   ├── all_samples.raw.vcf.gz -> ../04_bcftools/all_samples.raw.vcf.gz
│   ├── all_samples.stringent.vcf.gz -> ../04_bcftools/all_samples.stringent.vcf.gz
│   ├── evaluate_snps_v4.py -> ../check/evaluate_snps_v4.py
│   ├── gentle_results
│   │   ├── Sample_A_gentle_false_negatives.tsv
│   │   ├── Sample_A_gentle_false_positives.tsv
│   │   ├── Sample_A_gentle_genotype_confusion.tsv
│   │   ├── Sample_B_gentle_false_negatives.tsv
│   │   ├── Sample_B_gentle_false_positives.tsv
│   │   ├── Sample_B_gentle_genotype_confusion.tsv
│   │   ├── Sample_C_gentle_false_negatives.tsv
│   │   ├── Sample_C_gentle_false_positives.tsv
│   │   ├── Sample_C_gentle_genotype_confusion.tsv
│   │   ├── Sample_D_gentle_false_negatives.tsv
│   │   ├── Sample_D_gentle_false_positives.tsv
│   │   ├── Sample_D_gentle_genotype_confusion.tsv
│   │   └── summary_all_runs.tsv
│   ├── raw_results
│   │   ├── Sample_A_raw_false_negatives.tsv
│   │   ├── Sample_A_raw_false_positives.tsv
│   │   ├── Sample_A_raw_genotype_confusion.tsv
│   │   ├── Sample_B_raw_false_negatives.tsv
│   │   ├── Sample_B_raw_false_positives.tsv
│   │   ├── Sample_B_raw_genotype_confusion.tsv
│   │   ├── Sample_C_raw_false_negatives.tsv
│   │   ├── Sample_C_raw_false_positives.tsv
│   │   ├── Sample_C_raw_genotype_confusion.tsv
│   │   ├── Sample_D_raw_false_negatives.tsv
│   │   ├── Sample_D_raw_false_positives.tsv
│   │   ├── Sample_D_raw_genotype_confusion.tsv
│   │   └── summary_all_runs.tsv
│   ├── stringent_results
│   │   ├── Sample_A_stringent_false_negatives.tsv
│   │   ├── Sample_A_stringent_false_positives.tsv
│   │   ├── Sample_A_stringent_genotype_confusion.tsv
│   │   ├── Sample_B_stringent_false_negatives.tsv
│   │   ├── Sample_B_stringent_false_positives.tsv
│   │   ├── Sample_B_stringent_genotype_confusion.tsv
│   │   ├── Sample_C_stringent_false_negatives.tsv
│   │   ├── Sample_C_stringent_false_positives.tsv
│   │   ├── Sample_C_stringent_genotype_confusion.tsv
│   │   ├── Sample_D_stringent_false_negatives.tsv
│   │   ├── Sample_D_stringent_false_positives.tsv
│   │   ├── Sample_D_stringent_genotype_confusion.tsv
│   │   └── summary_all_runs.tsv
│   └── truth_snps.vcf -> ../check/truth_snps.vcf
├── check -> /lustre/isaac24/proj/UTK0505/test2/check
├── link_check.sh
├── link_data.sh
├── link_reads.sh
├── reads -> /lustre/isaac24/proj/UTK0505/test2/reads
└── ref -> /lustre/isaac24/proj/UTK0505/test2/ref

18 directories, 221 files

Results

Sample	Filter	TP	FP	FN	Precision	Recall	F1
A	raw	4,875	0	125	1.0000	0.9750	0.9873
A	gentle	4,864	0	136	1.0000	0.9728	0.9862
A	strict	4,772	0	228	1.0000	0.9544	0.9767
B	raw	4,422	23	578	0.9948	0.8844	0.9364
B	gentle	4,419	0	581	1.0000	0.8838	0.9383
B	strict	4,347	0	653	1.0000	0.8694	0.9301
C	raw	4,871	2	129	0.9996	0.9742	0.9867
C	gentle	4,865	0	135	1.0000	0.9730	0.9863
C	strict	4,772	0	228	1.000	0.9544	0.9767
D	raw	4,859	15	141	0.9969	0.9718	0.9842
D	gentle	4,851	1	149	0.9998	0.9702	0.9848
D	strict	4,761	0	239	1.0000	0.9522	0.9755