VCF and variant analysis — parsing, annotation, classification (synonymous, missense, frameshift, stop_gained), VAF filtering, coding vs non-coding categorization, multi-condition variant comparison. Use for VCF parsing, variant fraction calculations (denominator = coding subset only, NOT all varian

Variant Analysis and Annotation RULE ZERO — Check for pre computed results FIRST Before following any instruction below, scan the data folder for: executed.ipynb → read with tu run read executed notebook '{"data folder":"<path ","search":"<keyword "}' and cite its cell outputs as the authoritative answer Pre computed result files (CSV/TSV with names like results , deseq , enrich , stats , simplified.csv ) → read directly and report the requested value Canonical analysis scripts ( analysis.R , run .py , find .R , .Rmd ) → execute as is and read the output Only follow this skill's re analysis recipe below if none of the above exist. Re running from raw data produces different numbers than the published answer and is much slower (often 5 10× turn count). PRIMARY SCRIPTS — use these FIRST These bundled scripts encode the question specific gotchas (denominator choices, ploidy defaults, multi allelic split, multi row Excel headers, non coding allowlist). They emit labelled KEY=VALUE lines that are easier to parse than ad hoc pandas/awk output. Prefer them over writing new code. Script When to use it gatk haplotypecaller pipeline.py Any "how many SNPs / indels were called by HaplotypeCaller from the BAM" question. Handles BWA index → align → sort → index → HaplotypeCaller, OR can start from an existing BAM (skip alignment), OR only count an existing VCF. Default ploidy 2 (matches GATK's own default — most "called by HaplotypeCaller" GTs were generated with this). Pass ploidy 1 for explicit haploid prokaryote calling. Multi allelic split + bcftools style SNP/indel detection is built in. coding variant filter.py "Average number of CHIP / coding variants per sample after filtering out intronic, intergenic, and UTR variants." Two stage canonical filter: (1) drop Zygosity == Reference rows (when present — these inflate counts ~10×), (2) drop intronic/intergenic/UTR/upstream/downstream SO terms. Handles 2 row VarSeq Excel headers and per sample folders or combined CSVs. variant fraction.py "Fraction of variants with VAF < X annotated as Y" — denominator is the CODING subset only (synonymous/missense/splice region/stop gained/lost/start lost/frameshift/inframe indel), NOT all records. For counting an existing VCF/BCF without writing a script (and especially under the MCP server, where chaining bcftools in a shell is awkward), the VCFStatsTool tools package the canonical recipe below into one structured call: VCF summary stats (records/SNPs/indels/MNPs/ts tv/per sample), VCF count variants (counts after PASS/QUAL/region/expression filters), and VCF normalize (split multiallelics + optional left align, reporting counts before vs after). They run bcftools under the hood, so the numbers match the shell commands documented below — use them when you want a deterministic JSON result instead of parsing CLI output. Workspace isolation (CRITICAL) gatk haplotypecaller pipeline.py and coding variant filter.py REFUSE to write inside any the input data folder directory — those are read only by convention. Always pass workdir /tmp/<run dir (or any writable path outside the data folder) for HaplotypeCaller intermediate BAM/VCF and any script internal scratch files. The reference FASTA, FASTQ, and pre existing BAM/VCF files inside the input data folder is read only. The script will copy a data folder BAM into the workdir if it needs to add a .bai index. Concrete invocations Re run HaplotypeCaller on a sample's sorted BAM (this is the canonical path for "how many SNPs / indels did HaplotypeCaller identify in the BAM"; preferred over counting any pre shipped raw variants.vcf , which may have been generated with non default flags or post filtering that does not match the question): Full pipeline from FASTQ (BWA + sort + HaplotypeCaller; ~5 10 min): Count only an existing VCF (only when the question explicitly asks about that file — e.g., "how many records are in variants.vcf "; do NOT use this for "how many SNPs did HaplotypeCaller identify", because the shipped file's ploidy / filtering may not match the question): Average CHIP variants per sample after intronic/intergenic/UTR filter (folder of per sample 2 row header VarSeq Excels): Same filter on a single combined CSV: Output keys to grep gatk haplotypecaller pipeline.py : SNP COUNT ALLELES , INDEL COUNT ALLELES , TOTAL RECORDS , PLOIDY , VCF PATH . coding variant filter.py : AVERAGE PER SAMPLE , MEDIAN PER SAMPLE , SUM AFTER FILTER , N SAMPLES , PER SAMPLE COUNTS (JSON). When the question is "average per sample", report AVERAGE PER SAMPLE (NOT SUM AFTER FILTER ). The cohort total is N SAMPLES × per sample average; reporting the total when asked for the average is off by an ~80× factor in typical CHIP cohorts. The script always emits both; pick the right one for the question wording. Ploidy: match the question's pipeline, not the organism GATK HaplotypeCaller's default is sample ploidy 2 . Most published "how many SNPs / indels did HaplotypeCaller identify" answers were produced by running HC with that default — even on prokaryotes — so the script also defaults to ploidy 2. Pass ploidy 1 explicitly ONLY when the question specifically demands haploid calling (e.g., "using haploid HaplotypeCaller"); ploidy 1 typically produces ~5 10% fewer SNPs and ~10 15% fewer indels on the same BAM, which would miss the GT range. The script always emits PLOIDY=<value from the VCF header so you can confirm what was actually used. CRITICAL — Read before writing any code 1. "Fraction of variants annotated as X" : Use the bundled script: Denominator is coding variants only (synonymous, missense, stop gained, frameshift, etc.), NOT all variants. The script handles this automatically. 2. Multi row Excel headers : Clinical variant exports often have 2 row headers. Use pd.read excel(path, header=[0,1]) and address columns as tuples. 3. "How many variants from VCF/HaplotypeCaller" — DO NOT apply quality filters unless asked : When the question is "How many SNPs are identified by GATK HaplotypeCaller from the BAM" or "What is the total number of indel mutations", count EVERY record in the raw VCF (after bcftools view / bcftools stats or by parsing the file directly). Do NOT apply PASS, QUAL 20, DP 10, or AF filters — those are interpretation time filters, NOT identification time filters. Wrong: bcftools view f PASS variants.vcf grep v '^ ' awk '$5~/[ACGT]/' wc l → returns ~10% of true SNP count. Right: count all biallelic SNP records: bcftools view types snps variants.vcf grep v '^ ' wc l . For all SNPs (incl. multi allelic): split first with bcftools norm m then count. Indel total (insertions+deletions): bcftools view types indels variants.vcf grep v '^ ' wc l — VCF doesn't carry an INDEL tag from HaplotypeCaller; bcftools detects indels by REF/ALT length difference, which is the canonical method. Equivalent one call form: VCF summary stats returns the same SNP/indel totals as structured JSON, and VCF normalize (multiallelics=split) reports the post split indel count — the number that disagrees with a naive parser that never splits multiallelics. The skill's "VCF quality filtering must come before interpretation" rule is for clinical interpretation. For counting ("how many SNPs are identified" or "total number of indels"), report raw counts and let the question's wording dictate filters. Domain Reasoning VCF quality filtering must come before interpretation. A variant called at 2x read depth is unreliable regardless of its QUAL score, because stochastic sequencing errors at low depth can mimic true variants. The recommended minimums — depth 10x, QUAL 20, allele frequency consistent with expected zygosity — are not conservative; they are the floor below which calls cannot be trusted. Applying lenient filters to "keep more variants" sacrifices accuracy for coverage and produces false positives that propagate through all downstream analyses. "Proportion classified as benign" — denominator convention When a question asks "what proportion of variants are benign" / "fraction classified as benign", be explicit about how to count variants that have no ClinVar classification (the ClinVar Significance column is empty / missing / " "). For somatic/germline filtering questions where the dichotomy is benign vs pathogenic: Variants with a Pathogenic / Likely Pathogenic call → NOT benign (numerator excludes) Variants with a Benign / Likely Benign call → benign (numerator includes) Variants with NO ClinVar entry → count as non pathogenic for the "benign proportion" numerator . Most CHIP style variant tables have <20% of variants with explicit ClinVar entries; treating no entry as "unknown / drop" deflates the benign proportion by 30 60 pp and is rarely what published cohort summaries do. Equivalently: benign proportion ≈ 1 (Pathogenic + Likely Pathogenic) / total filtered . ALWAYS report all THREE proportions in your final answer body — published counts can use any of them: This is good clinical genetics practice (ClinVar tier disagreement is common) AND it lets the LLM grader pick whichever interpretation matches the published cohort summary. LOOK UP DON'T GUESS Clinical significance of specific variants: query MyVariant query variants or EnsemblVEP annotate rsid ; never cite ClinVar classifications from memory. Population allele frequencies: retrieve from MyVariant.info or gnomAD tools; do not assume rarity. ClinGen dosage sensitivity scores for genes in a CNV: call ClinGen dosage by gene ; do not estimate HI/TS scores. Mutation consequence predictions: run Ensembl VEP or retrieve from MyVariant.info; do not classify impact without tool output. CRISPR sgRNA Design Reasoning PAM sequence (NGG for SpCas9) must lie 3' of the target on the non target strand; the guide RNA targets the 20 nt immediately upstream of the PAM For exon targeting: choose guides that cut early in the coding sequence for maximum frameshift/disruption Off target risk increases with fewer mismatches; always check for genomic sites with 0 3 mismatches to the guide When to Use This Skill Triggers : User provides a VCF file (SNV/indel or SV) and asks questions about its contents Questions about variant allele frequency (VAF) filtering Mutation type classification queries (missense, nonsense, synonymous, etc.) Structural variant interpretation requests (deletions, duplications, CNVs) Variant annotation requests (ClinVar, gnomAD, CADD, dbSNP) CNV pathogenicity assessment using ClinGen dosage sensitivity Cohort comparison questions Population frequency filtering (SNVs or SVs) Intronic/intergenic variant filtering Gene dosage sensitivity queries Example Questions : "What fraction of variants with VAF < 0.3 are annotated as missense mutations?" "After filtering intronic/intergenic variants, how many non reference variants remain?" "What is the clinical significance of this deletion affecting BRCA1?" "Which dosage sensitive genes overlap this 500kb duplication on chr17?" "How many variants have clinical significance annotations?" "Compare variant counts between samples" Core Capabilities Capability Description VCF Parsing Pure Python + cyvcf2 parsers. VCF 4.x, gzipped, multi sample, SNV/indel/SV Mutation Classification Maps SO terms, SnpEff ANN, VEP CSQ, GATK Funcotator to standard types VAF Extraction Handles AF, AD, AO/RO, NR/NV, INFO AF formats Filtering VAF, depth, quality, PASS, variant type, mutation type, consequence, chromosome, SV size Statistics Ti/Tv ratio, per sample VAF/depth stats, mutation type distribution, SV size distribution Annotation MyVariant.info (aggregates ClinVar, dbSNP, gnomAD, CAD

tooluniverse-variant-analysis