bio-vcf-manipulation

What this skill does

## Version Compatibility

Reference examples tested with: GATK 4.5+, bcftools 1.19+

Before using code patterns, verify installed versions match. If versions differ:
- Python: `pip show <package>` then `help(module.function)` to check signatures
- CLI: `<tool> --version` then `<tool> --help` to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed
package and adapt the example to match the actual API rather than retrying.

# VCF Manipulation

Merge, concat, sort, and compare VCF files using bcftools.

## Operations Overview

| Operation | Command | Use Case |
|-----------|---------|----------|
| Merge | `bcftools merge` | Combine samples from multiple VCFs |
| Concat | `bcftools concat` | Combine regions from multiple VCFs |
| Sort | `bcftools sort` | Sort unsorted VCF |
| Intersect | `bcftools isec` | Compare/intersect call sets |
| Subset | `bcftools view` | Extract samples or regions |

## Merge vs Concat Decision

| Operation | Input structure | Result | Requires index |
|-----------|----------------|--------|----------------|
| `bcftools merge` | Different **samples** at (potentially) same positions | Multi-sample VCF | Yes |
| `bcftools concat` | Same **samples** from different **regions** | Combined-region VCF | No (unless `-a`) |

Common mistake: using concat when merge is needed (or vice versa). `concat --allow-overlaps` (`-a`) can handle some overlap cases but does not merge genotypes across samples -- it only resolves duplicate records from the same sample.

## bcftools merge

**Goal:** Combine VCF files from different samples into a single multi-sample VCF.

**Approach:** Use bcftools merge to join files with different sample columns at shared genomic positions.

**"Merge my per-sample VCFs into one file"** -> Combine variant records from multiple samples into a single multi-sample VCF.

Combine multiple VCF files with **different samples** at the same positions.

### Basic Merge

```bash
bcftools merge sample1.vcf.gz sample2.vcf.gz -Oz -o merged.vcf.gz
```

### Merge Multiple Files

```bash
bcftools merge *.vcf.gz -Oz -o all_samples.vcf.gz
```

### Merge from File List

```bash
# files.txt: one VCF path per line
bcftools merge -l files.txt -Oz -o merged.vcf.gz
```

### Handle Missing Genotypes

```bash
# Output missing genotypes as ./. (default)
bcftools merge sample1.vcf.gz sample2.vcf.gz -Oz -o merged.vcf.gz

# Output missing as reference (0/0)
bcftools merge --missing-to-ref sample1.vcf.gz sample2.vcf.gz -Oz -o merged.vcf.gz
```

### Force Sample Names

When sample names conflict:

```bash
bcftools merge --force-samples sample1.vcf.gz sample2.vcf.gz -Oz -o merged.vcf.gz
```

### Multiallelic Handling During Merge

When merging VCFs, if the same site has different ALT alleles across files, bcftools merge creates a multiallelic record by default. This is usually correct behavior. Use `--merge none` to keep records separate, or `--merge snps`/`--merge indels` to control merging granularity.

```bash
# Default: merge all variant types at same position
bcftools merge sample1.vcf.gz sample2.vcf.gz -Oz -o merged.vcf.gz

# Keep SNPs and indels at same position as separate records
bcftools merge --merge snps sample1.vcf.gz sample2.vcf.gz -Oz -o merged.vcf.gz

# No merging -- every record stays separate
bcftools merge --merge none sample1.vcf.gz sample2.vcf.gz -Oz -o merged.vcf.gz
```

### Merge Specific Regions

```bash
bcftools merge -r chr1:1000000-2000000 sample1.vcf.gz sample2.vcf.gz -Oz -o merged.vcf.gz
```

## bcftools concat

**Goal:** Concatenate VCF files that cover different genomic regions for the same samples.

**Approach:** Use bcftools concat to join region-split files (e.g., per-chromosome VCFs) in order.

Combine VCF files with **same samples** from different regions.

### Concatenate Chromosomes

```bash
bcftools concat chr1.vcf.gz chr2.vcf.gz chr3.vcf.gz -Oz -o genome.vcf.gz
```

### Concatenate All Chromosomes

```bash
bcftools concat chr*.vcf.gz -Oz -o genome.vcf.gz
```

### From File List

```bash
# files.txt: one VCF path per line (in order)
bcftools concat -f files.txt -Oz -o concatenated.vcf.gz
```

### Allow Overlapping Regions

```bash
bcftools concat -a chr1_part1.vcf.gz chr1_part2.vcf.gz -Oz -o chr1.vcf.gz
```

### Remove Duplicates

```bash
bcftools concat -a -d all file1.vcf.gz file2.vcf.gz -Oz -o merged.vcf.gz
```

Options for `-d`:
- `snps` - Remove duplicate SNPs
- `indels` - Remove duplicate indels
- `both` - Remove duplicate SNPs and indels
- `all` - Remove all duplicates
- `exact` - Remove exact duplicates only

## bcftools sort

**Goal:** Sort a VCF file by chromosome and position.

**Approach:** Use bcftools sort with optional temp directory and memory limits for large files.

Sort VCF by chromosome and position.

### Basic Sort

```bash
bcftools sort input.vcf -Oz -o sorted.vcf.gz
```

### With Temporary Directory

For large files:

```bash
bcftools sort -T /tmp input.vcf.gz -Oz -o sorted.vcf.gz
```

### Memory Limit

```bash
bcftools sort -m 4G input.vcf.gz -Oz -o sorted.vcf.gz
```

## bcftools isec

**Goal:** Identify shared and private variants between two or more VCF files.

**Approach:** Use bcftools isec to partition variants into private-to-each-file and shared subsets.

**"Find variants called by both GATK and bcftools"** -> Intersect two call sets to identify concordant and discordant variants.

Intersect and compare VCF files.

**Normalize before comparing:** Before intersecting VCF files from different callers, ALWAYS normalize both files first. Different callers may represent the same variant differently (e.g., different indel left-alignment, multiallelic vs biallelic representation). Without normalization, bcftools isec will report false discordance. See variant-normalization for details.

```bash
bcftools norm -m-any -f ref.fa gatk.vcf.gz -Oz -o gatk_norm.vcf.gz
bcftools norm -m-any -f ref.fa bcftools_calls.vcf.gz -Oz -o bcftools_norm.vcf.gz
bcftools isec -p comparison gatk_norm.vcf.gz bcftools_norm.vcf.gz
```

### Find Shared Variants

```bash
bcftools isec -p output_dir sample1.vcf.gz sample2.vcf.gz
```

Creates:
- `0000.vcf` - Private to sample1
- `0001.vcf` - Private to sample2
- `0002.vcf` - Shared (sample1 records)
- `0003.vcf` - Shared (sample2 records)

### Output Compressed

```bash
bcftools isec -p output_dir -Oz sample1.vcf.gz sample2.vcf.gz
```

### Intersection Only

```bash
bcftools isec -p output_dir -n=2 sample1.vcf.gz sample2.vcf.gz
# Only outputs variants present in exactly 2 files
```

### Comparison Options

| Flag | Description |
|------|-------------|
| `-n=2` | Present in exactly 2 files |
| `-n+2` | Present in 2 or more files |
| `-n-2` | Present in fewer than 2 files |
| `-n~11` | Boolean: file1 AND file2 |
| `-n~10` | Boolean: file1 AND NOT file2 |

### Two-File Intersection

The `-w` flag selects which file's records appear in the output. `-w1` outputs records from the first file, `-w2` from the second. This matters when the two files carry different INFO/FORMAT annotations -- the output inherits whichever file's records are selected.

```bash
# Variants in both files, output records from file 1 (with file 1 annotations)
bcftools isec -n=2 -w1 sample1.vcf.gz sample2.vcf.gz -Oz -o shared.vcf.gz

# Same intersection, but output records from file 2
bcftools isec -n=2 -w2 sample1.vcf.gz sample2.vcf.gz -Oz -o shared_file2_annot.vcf.gz

# Variants only in sample1
bcftools isec -n~10 -w1 sample1.vcf.gz sample2.vcf.gz -Oz -o only_sample1.vcf.gz
```

### Complement Mode

```bash
# Variants in file1 not in file2
bcftools isec -C sample1.vcf.gz sample2.vcf.gz -Oz -o unique.vcf.gz
```

## Subsetting VCF Files

**Goal:** Extract a subset of samples or regions from a multi-sample VCF.

**Approach:** Use bcftools view with -s (samples) or -r/-R (regions) flags to create targeted subsets.

### Extract Samples

```bash
bcftools view -s sample1,sample2 input.vcf.gz -Oz -o subset.vcf.gz
```

### Exclude Samples

```bash
bcftools view -s