### Install mapDamage using pip

Source: https://github.com/ginolhac/mapdamage/blob/main/README.md

This command installs the mapDamage package using pip, the Python package installer. Ensure you have Python and pip installed on your system.

```bash
pip install mapdamage
```

--------------------------------

### Install mapDamage and Check R Packages

Source: https://context7.com/ginolhac/mapdamage/llms.txt

Installs the mapDamage tool using pip and verifies the availability of required R packages for plotting and statistical analysis.

```bash
pip install mapdamage

mapDamage --check-R-packages
```

--------------------------------

### DNA Composition Around Reads Output Format (Bash)

Source: https://context7.com/ginolhac/mapdamage/llms.txt

Explains the 'dnacomp.txt' file format, which shows DNA composition around reads, with columns for sample, library, strand, and position relative to the read.

```bash
# dnacomp.txt - DNA composition around reads
# Tab-separated with columns:
# Sample  Library  End  Std  Pos  A  C  G  T  Total
# Negative Pos = upstream of read, Positive Pos = within read
```

--------------------------------

### Calculate and Write Base Composition (Python)

Source: https://context7.com/ginolhac/mapdamage/llms.txt

Calculates genome-wide base composition from a reference FASTA file and writes it to a CSV file. Also shows how to read previously computed base composition.

```python
from mapdamage.composition import write_base_comp, read_base_comp
from pathlib import Path

# Calculate and write base composition from reference FASTA
write_base_comp(
    fasta=Path("reference.fasta"),
    destination=Path("results/dnacomp_genome.csv")
)
# Output CSV format:
# A,C,G,T
# 0.295,0.205,0.205,0.295

# Read previously computed base composition
base_comp = read_base_comp(Path("results/dnacomp_genome.csv"))
# Returns: {'A': '0.295', 'C': '0.205', 'G': '0.205', 'T': '0.295'}
```

--------------------------------

### Fragment Length Distribution Output Format (Bash)

Source: https://context7.com/ginolhac/mapdamage/llms.txt

Details the 'lgdistribution.txt' file format, which summarizes fragment length distributions for paired-end and single-end libraries.

```bash
# lgdistribution.txt - Fragment length distribution
# Tab-separated with columns:
# Sample  Library  Std  Kind  Length  Occurrences
# Kind = pe (paired-end) or se (single-end)
```

--------------------------------

### Bayesian Statistical Estimation of Damage Parameters with mapDamage

Source: https://context7.com/ginolhac/mapdamage/llms.txt

Conducts Bayesian estimation of DNA damage parameters using MCMC methods. This process estimates damage rates, overhang lengths, and generates posterior probability distributions. Custom MCMC parameters and specific damage models (e.g., single-stranded, different overhangs, substitution models) can be applied. It can also be run solely on existing results.

```bash
# Full analysis with Bayesian estimation (default behavior)
mapDamage -i ancient_sample.bam -r reference.fasta

# Custom MCMC parameters for more thorough estimation
mapDamage -i ancient_sample.bam -r reference.fasta \
    --burn 20000 \
    --iter 100000 \
    --rand 50

# Use only 5' end mismatches for damage model
mapDamage -i ancient_sample.bam -r reference.fasta --termini 5p

# Use only 3' end mismatches
mapDamage -i ancient_sample.bam -r reference.fasta --termini 3p

# Single-stranded library protocol (different damage patterns)
mapDamage -i ss_library.bam -r reference.fasta --single-stranded

# Allow different overhang parameters for 5' and 3' ends
mapDamage -i ancient_sample.bam -r reference.fasta --diff-hangs

# Use Jukes-Cantor substitution model instead of HKY85
mapDamage -i ancient_sample.bam -r reference.fasta --jukes-cantor

# Disable Bayesian estimation (only compute damage patterns)
mapDamage -i ancient_sample.bam -r reference.fasta --no-stats

# Run only statistical estimation on existing results folder
mapDamage --stats-only -d results_ancient_sample -r reference.fasta
```

--------------------------------

### Basic Damage Pattern Analysis with mapDamage

Source: https://context7.com/ginolhac/mapdamage/llms.txt

Performs a complete DNA damage analysis, including misincorporation patterns and fragment length distributions. It can optionally include Bayesian estimation and generates PDF plots and tabular output files. Requires a BAM alignment file and a reference genome FASTA file.

```bash
# Basic analysis with BAM file and reference genome
mapDamage -i aligned_reads.bam -r reference.fasta

# Specify output directory and plot title
mapDamage -i ancient_sample.bam -r hg38.fasta -d results_ancient_sample -t "Ancient Sample Analysis"

# Analyze with downsampling to 100,000 reads for faster processing
mapDamage -i large_dataset.bam -r reference.fasta -n 100000 --downsample-seed 42

# Analyze first 100 nucleotides from read termini with 15bp flanking regions
mapDamage -i aligned_reads.bam -r reference.fasta -l 100 -a 15

# Filter bases with quality below Phred 20
mapDamage -i aligned_reads.bam -r reference.fasta -Q 20
```

--------------------------------

### Align Sequence with Quality Filtering (Python)

Source: https://context7.com/ginolhac/mapdamage/llms.txt

Shows how to align a sequence with quality filtering, masking bases below a specified Phred score threshold with 'N'.

```python
from mapdamage.align import align_with_qual

# Align with quality filtering
cigar = [(0, 3), (2, 1), (0, 4)]  # 3M1D4M
seq, qual, refseq = align_with_qual(
    cigarlist=cigar,
    seq="ATCGATCG",
    qual="IIIIIIII",  # ASCII quality scores
    threshold=20,     # Minimum Phred score
    ref="ATCATCG"
)
# Bases below threshold are masked as 'N'
```

--------------------------------

### Extract Coordinates and Reference Sequence Flanking Alignment (Python)

Source: https://context7.com/ginolhac/mapdamage/llms.txt

Demonstrates how to extract the 5' and 3' coordinates of a read and retrieve the reference sequence before and after the alignment using pysam and custom functions.

```python
from mapdamage.coordinates import get_coordinates, get_around
import pysam

# Assuming 'read' is a pysam AlignedSegment object
fivep, threep = get_coordinates(read)

# Example for forward reads: fivep = read.pos, threep = read.aend
# Example for reverse reads: fivep = read.aend, threep = read.pos

# Get reference sequence flanking the alignment
reflengths = {'chr1': 248956422}
pysam_fasta_handle = pysam.FastaFile("reference.fasta")
before, after = get_around(
    coord=(fivep, threep),
    chrom="chr1",
    reflengths=reflengths,
    length=10,
    ref=pysam_fasta_handle
)
# Returns 10bp before and after the alignment
```

--------------------------------

### Reading from Standard Input

Source: https://context7.com/ginolhac/mapdamage/llms.txt

Processes BAM data piped from standard input, enabling integration with other command-line tools like samtools.

```APIDOC
## POST /mapDamage/stdin

### Description
Processes BAM data directly from standard input (stdin), allowing for seamless integration with other bioinformatics tools that can pipe BAM output. Note that the `--rescale` option is not supported with piped input due to its two-pass requirement.

### Method
POST

### Endpoint
/mapDamage/stdin

### Parameters
#### Query Parameters
- **-i** (string) - Required - Specifies input source. Use `-` to indicate standard input.
- **-r** (string) - Required - Reference FASTA file path.
- **-d** (string) - Required - Output directory for results.

### Request Example
```bash
samtools view -b filtered.bam | mapDamage -i - -r reference.fasta -d results_filtered
```

### Response
#### Success Response (200)
- **results_directory** (string) - Path to the directory where analysis results are stored.

#### Response Example
```json
{
  "results_directory": "results_filtered"
}
```
```

--------------------------------

### Damage Probabilities for Rescaling Output Format (Bash)

Source: https://context7.com/ginolhac/mapdamage/llms.txt

Describes the 'Stats_out_MCMC_correct_prob.csv' file format, which contains damage probabilities used for quality score rescaling, including position and C>T/G>A probabilities.

```bash
# Stats_out_MCMC_correct_prob.csv - Damage probabilities for rescaling
# CSV with columns:
# Position,C.T,G.A
# Position ranges from -seq_length to +seq_length
# Values are damage probabilities used for quality rescaling
```

--------------------------------

### Misincorporation Frequency Output Format (Bash)

Source: https://context7.com/ginolhac/mapdamage/llms.txt

Describes the format of the 'misincorporation.txt' file, which details base misincorporation frequencies, including sample, library, strand, position, and mutation types.

```bash
# misincorporation.txt - Base misincorporation frequencies
# Tab-separated with columns:
# Sample  Library  End  Std  Pos  A  C  G  T  Total  G>A  C>T  A>G  T>C  ...
# Pos = position from 5' (5p) or 3' (3p) terminus
# Std = strand (+ or -)
# Mutations listed as Ref>Read (e.g., C>T)
```

--------------------------------

### Parse CIGAR String for Operations (Python)

Source: https://context7.com/ginolhac/mapdamage/llms.txt

Demonstrates how to parse a CIGAR string to extract information about specific operations, such as soft clips, and their positions.

```python
from mapdamage.cigar import parse_cigar

# Parse CIGAR for specific operations
# Operation codes: 0=M, 1=I, 2=D, 3=N, 4=S, 5=H, 6=P, 7==, 8=X
cigar = [(4, 5), (0, 50), (4, 3)]  # 5S50M3S
soft_clips = parse_cigar(cigar, 4) # 4 corresponds to soft clip (S)
# Returns: [[5, 0], [3, 55]] - [(length, position), ...]
```

--------------------------------

### Statistics Classes for DNA Damage Analysis in Python

Source: https://context7.com/ginolhac/mapdamage/llms.txt

This module provides Python classes for collecting and analyzing DNA damage statistics, including misincorporation rates, DNA composition around alignments, and fragment length distributions. These classes require initialized BAMReader and library information. Statistics are updated per read and can be written to files.

```python
from mapdamage.statistics import MisincorporationRates, DNAComposition, FragmentLengths
from mapdamage.reader import BAMReader
from pathlib import Path

# Initialize statistics collectors
reader = BAMReader(filepath=Path("aligned_reads.bam"))
libraries = reader.get_libraries()

# Misincorporation rate tracking (70bp from each terminus)
misincorp = MisincorporationRates(libraries=libraries, length=70)

# DNA composition around alignments
dnacomp = DNAComposition(libraries=libraries, around=10, length=70)

# Fragment length distribution
lgdistrib = FragmentLengths(libraries=libraries)

# Process reads and update statistics
for read in reader:
    library = reader.get_sample_and_library(read)

    # Update fragment length statistics
    lgdistrib.update(read, library)

    # For misincorporation and composition, you also need reference sequence alignment
    # (See main.py for complete alignment workflow)

# Write results to files
misincorp.write(Path("results/misincorporation.txt"))
dnacomp.write(Path("results/dnacomp.txt"))
lgdistrib.write(Path("results/lgdistribution.txt"))

reader.close()
```

--------------------------------

### Align Sequence with Gaps Based on CIGAR String (Python)

Source: https://context7.com/ginolhac/mapdamage/llms.txt

Illustrates how to align a read sequence against a reference sequence, inserting gaps to account for deletions or insertions as defined by a CIGAR string.

```python
from mapdamage.align import align

# Insert gaps according to CIGAR string
seq = "ATCGATCG"  # read sequence
refseq = "ATCATCG"  # reference (has deletion)
cigar = [(0, 3), (2, 1), (0, 4)]  # 3M1D4M
aligned_seq, aligned_ref = align(cigar, seq, refseq)
# aligned_seq: "ATC-ATCG" (gap inserted at deletion)
# aligned_ref: "ATCGATCG" (gap removed)
```

--------------------------------

### Python API - Sequence Utilities

Source: https://context7.com/ginolhac/mapdamage/llms.txt

Provides essential DNA sequence manipulation functions, including reverse complementing, FASTA indexing, and sequence dictionary comparison.

```APIDOC
## Python API - Sequence Utilities

### Description
A collection of utility functions for common DNA sequence manipulations required in ancient DNA analysis, such as calculating reverse complements and comparing sequence dictionaries derived from FASTA and BAM files.

### Functions
#### `revcomp(seq: str) -> str`
Calculates the reverse complement of a DNA sequence.
- `seq` (str): The input DNA sequence.
Returns: The reverse complemented DNA sequence.

#### `read_fasta_index(fai_path: str) -> Dict[str, int]`
Reads a FASTA index file (.fai) and returns a dictionary of reference sequences and their lengths.
- `fai_path` (str): Path to the FASTA index file.
Returns: A dictionary mapping reference names to their lengths.

#### `compare_sequence_dicts(fasta_refs: Dict[str, int], bam_refs: Dict[str, int]) -> bool`
Compares sequence dictionaries from FASTA and BAM files to ensure compatibility.
- `fasta_refs` (Dict[str, int]): Dictionary of reference sequences and lengths from FASTA.
- `bam_refs` (Dict[str, int]): Dictionary of reference sequences and lengths from BAM.
Returns: True if all BAM references are present in the FASTA dictionary with matching lengths, False otherwise. Logs warnings for extra FASTA sequences and errors for missing or mismatched references.

### Request Example
```python
from mapdamage.seq import revcomp, read_fasta_index, compare_sequence_dicts

# Reverse complement a sequence
seq = "ATCGATCG"
rc_seq = revcomp(seq)
# Returns: "CGATCGAT"

# Handles ambiguous bases (IUPAC codes)
ambig_seq = "ATCGMRWSYKVHDB"
rc_ambig = revcomp(ambig_seq)
# Returns: "VDHBMRSWYKATCG"

# Read FASTA index file
fai_dict = read_fasta_index("reference.fasta.fai")
# Returns: {'chr1': 248956422, 'chr2': 242193529, ...}

# Compare BAM and FASTA sequence dictionaries
bam_refs = {'chr1': 248956422, 'chr2': 242193529}
fasta_refs = {'chr1': 248956422, 'chr2': 242193529, 'chrM': 16569}
is_compatible = compare_sequence_dicts(fasta_refs, bam_refs)
# Returns True if all BAM references found in FASTA with matching lengths
# Logs warnings for extra FASTA sequences, errors for missing/mismatched
```

### Response
#### Success Response (200)
- The utility functions return the results of the sequence operations (e.g., reverse complement string, dictionary of references, boolean for compatibility).

#### Response Example
```json
{
  "reverse_complement": "CGATCGAT",
  "fasta_references": {"chr1": 248956422},
  "sequences_compatible": true
}
```
```

--------------------------------

### Sequence Utilities for DNA Manipulation in Python

Source: https://context7.com/ginolhac/mapdamage/llms.txt

Provides essential DNA sequence manipulation functions, including reverse complementing, reading FASTA index files, and comparing sequence dictionaries between BAM and FASTA references. It handles ambiguous bases (IUPAC codes) and checks for compatibility between reference sequences.

```python
from mapdamage.seq import revcomp, read_fasta_index, compare_sequence_dicts

# Reverse complement a sequence
seq = "ATCGATCG"
rc_seq = revcomp(seq)
# Returns: "CGATCGAT"

# Handles ambiguous bases (IUPAC codes)
ambig_seq = "ATCGMRWSYKVHDB"
rc_ambig = revcomp(ambig_seq)
# Returns: "VDHBMRSWYKATCG"

# Read FASTA index file
fai_dict = read_fasta_index("reference.fasta.fai")
# Returns: {'chr1': 248956422, 'chr2': 242193529, ...}

# Compare BAM and FASTA sequence dictionaries
bam_refs = {'chr1': 248956422, 'chr2': 242193529}
fasta_refs = {'chr1': 248956422, 'chr2': 242193529, 'chrM': 16569}
is_compatible = compare_sequence_dicts(fasta_refs, bam_refs)
# Returns True if all BAM references found in FASTA with matching lengths
# Logs warnings for extra FASTA sequences, errors for missing/mismatched
```

--------------------------------

### Process BAM Data from Standard Input with mapDamage

Source: https://context7.com/ginolhac/mapdamage/llms.txt

Allows mapDamage to process BAM data piped from standard input, enabling integration with other command-line tools like samtools. This is useful for performing analyses without writing intermediate BAM files. Note that the --rescale option cannot be used with piped input as it requires two passes.

```bash
# Pipe BAM data from samtools
samtools view -b filtered.bam | mapDamage -i - -r reference.fasta -d results_filtered
```

--------------------------------

### Python API - Statistics Classes

Source: https://context7.com/ginolhac/mapdamage/llms.txt

Classes for collecting and tracking various statistics related to DNA damage, including misincorporation rates, DNA composition, and fragment lengths.

```APIDOC
## Python API - Statistics Classes

### Description
These classes enable the collection and tracking of key statistics during the analysis of ancient DNA sequences. They cover misincorporation rates, DNA base composition around alignments, and fragment length distributions.

### Classes
#### `MisincorporationRates(libraries: List[Tuple[str, str]], length: int = 70)`
Tracks misincorporation rates within a specified window from each terminus of the reads.
- `libraries` (List[Tuple[str, str]]): List of libraries obtained from `BAMReader.get_libraries()`.
- `length` (int, optional): The number of bases from each terminus to consider for misincorporation analysis. Defaults to 70.

#### `DNAComposition(libraries: List[Tuple[str, str]], around: int = 10, length: int = 70)`
Analyzes DNA base composition within a specified window around alignment points.
- `libraries` (List[Tuple[str, str]]): List of libraries.
- `around` (int, optional): Number of bases around the alignment point to consider. Defaults to 10.
- `length` (int, optional): The total length of the window for composition analysis. Defaults to 70.

#### `FragmentLengths(libraries: List[Tuple[str, str]])`
Collects statistics on the distribution of DNA fragment lengths.
- `libraries` (List[Tuple[str, str]]): List of libraries.

### Methods (Common to Statistics Classes)
#### `update(self, read: AlignedSegment, library: Tuple[str, str])`
Updates the statistics collector with data from a single read.
- `read` (AlignedSegment): The aligned read object from `BAMReader`.
- `library` (Tuple[str, str]): The sample and library identifier for the read.

#### `write(self, filepath: Path)`
Writes the collected statistics to a file.
- `filepath` (Path): The path to the output file.

### Request Example
```python
from mapdamage.statistics import MisincorporationRates, DNAComposition, FragmentLengths
from mapdamage.reader import BAMReader
from pathlib import Path

# Initialize statistics collectors
reader = BAMReader(filepath=Path("aligned_reads.bam"))
libraries = reader.get_libraries()

# Misincorporation rate tracking (70bp from each terminus)
misincorp = MisincorporationRates(libraries=libraries, length=70)

# DNA composition around alignments
dnacomp = DNAComposition(libraries=libraries, around=10, length=70)

# Fragment length distribution
lgdistrib = FragmentLengths(libraries=libraries)

# Process reads and update statistics
for read in reader:
    library = reader.get_sample_and_library(read)

    # Update fragment length statistics
    lgdistrib.update(read, library)

    # For misincorporation and composition, you also need reference sequence alignment
    # (See main.py for complete alignment workflow)

# Write results to files
misincorp.write(Path("results/misincorporation.txt"))
dnacomp.write(Path("results/dnacomp.txt"))
lgdistrib.write(Path("results/lgdistribution.txt"))

reader.close()
```

### Response
#### Success Response (200)
- The statistics classes provide methods to update and write collected data. The `write` method generates output files containing the statistics.

#### Response Example
```json
{
  "misincorporation_file": "results/misincorporation.txt",
  "dnacomp_file": "results/dnacomp.txt",
  "fragment_length_file": "results/lgdistribution.txt"
}
```
```

--------------------------------

### Python API - BAMReader Class

Source: https://context7.com/ginolhac/mapdamage/llms.txt

Provides functionality to read and parse BAM files, including filtering of reads and downsampling capabilities.

```APIDOC
## Python API - BAMReader Class

### Description
The `BAMReader` class facilitates the parsing of BAM files. It supports filtering of reads based on standard criteria (unmapped, secondary, QC-failed, duplicates) and offers options for downsampling the dataset and merging libraries.

### Methods
#### `__init__(self, filepath: Path, downsample_to: Union[float, int] = None, downsample_seed: int = None, merge_libraries: bool = False)`
Initializes the `BAMReader`.
- `filepath` (Path): Path to the BAM file.
- `downsample_to` (Union[float, int], optional): Fraction (0.0-1.0) or number of reads to downsample to. Defaults to None.
- `downsample_seed` (int, optional): Seed for random downsampling. Defaults to None.
- `merge_libraries` (bool, optional): If True, merges all libraries into a single analysis. Defaults to False.

#### `get_references(self) -> Dict[str, int]`
Retrieves reference sequences and their lengths from the BAM file header.
Returns: A dictionary mapping reference names to their lengths.

#### `get_libraries(self) -> List[Tuple[str, str]]`
Retrieves library information (sample and library IDs) from the BAM file header.
Returns: A list of tuples, where each tuple contains (SampleName, LibraryID).

#### `get_sample_and_library(self, read: AlignedSegment) -> Tuple[str, str]`
Determines the sample and library for a given read based on BAM tags.
Returns: A tuple containing (SampleName, LibraryID).

#### `__iter__(self)`
Iterates through the filtered reads in the BAM file.

#### `close(self)`
Closes the BAM file reader.

### Request Example
```python
from pathlib import Path
from mapdamage.reader import BAMReader

# Basic BAM reading
reader = BAMReader(filepath=Path("aligned_reads.bam"))

# Get reference sequences and lengths
references = reader.get_references()
# Returns: {'chr1': 248956422, 'chr2': 242193529, ...}

# Get library information
libraries = reader.get_libraries()
# Returns: dict_keys([('SampleName', 'LibraryID'), ...])

# Iterate through filtered reads
for read in reader:
    library = reader.get_sample_and_library(read)
    print(f"Read: {read.query_name}, Library: {library}")

reader.close()

# Downsampling to fraction (10% of reads)
reader_downsampled_fraction = BAMReader(
    filepath=Path("large_dataset.bam"),
    downsample_to=0.1,
    downsample_seed=42
)

# Downsampling to fixed number (100,000 reads)
reader_downsampled_count = BAMReader(
    filepath=Path("large_dataset.bam"),
    downsample_to=100000,
    downsample_seed=42
)

# Merge all libraries into single analysis
reader_merged = BAMReader(
    filepath=Path("multi_library.bam"),
    merge_libraries=True
)
```

### Response
#### Success Response (200)
- The `BAMReader` object provides access to BAM file data and metadata.

#### Response Example
```json
{
  "references": {"chr1": 248956422, "chr2": 242193529},
  "libraries": [["SampleName", "LibraryID"]]
}
```
```

--------------------------------

### BAMReader Class for BAM File Parsing in Python

Source: https://context7.com/ginolhac/mapdamage/llms.txt

The BAMReader class in mapDamage's Python API facilitates parsing BAM files with built-in filtering (unmapped, secondary, QC-failed, duplicates) and downsampling capabilities. It allows retrieval of reference sequences, library information, and iteration over reads. Downsampling can be performed to a fraction or a fixed number of reads, and libraries can be merged.

```python
from pathlib import Path
from mapdamage.reader import BAMReader

# Basic BAM reading
reader = BAMReader(filepath=Path("aligned_reads.bam"))

# Get reference sequences and lengths
references = reader.get_references()
# Returns: {'chr1': 248956422, 'chr2': 242193529, ...}

# Get library information
libraries = reader.get_libraries()
# Returns: dict_keys([('SampleName', 'LibraryID'), ...])

# Iterate through filtered reads (excludes unmapped, secondary, QC-failed, duplicates)
for read in reader:
    library = reader.get_sample_and_library(read)
    print(f"Read: {read.query_name}, Library: {library}")

reader.close()

# Downsampling to fraction (10% of reads)
reader = BAMReader(
    filepath=Path("large_dataset.bam"),
    downsample_to=0.1,
    downsample_seed=42
)

# Downsampling to fixed number (100,000 reads)
reader = BAMReader(
    filepath=Path("large_dataset.bam"),
    downsample_to=100000,
    downsample_seed=42
)

# Merge all libraries into single analysis
reader = BAMReader(
    filepath=Path("multi_library.bam"),
    merge_libraries=True
)
```

--------------------------------

### Generate and Customize Damage Pattern Plots with mapDamage

Source: https://context7.com/ginolhac/mapdamage/llms.txt

Generates or regenerates DNA damage pattern plots from existing analysis results or by processing BAM files. Users can customize plot parameters such as the maximum y-axis value, the number of bases plotted from read termini, the number of reference bases plotted around reads, and the plot title.

```bash
# Regenerate plots from existing results folder
mapDamage --plot-only -d results_ancient_sample
```

```bash
# Customize plot parameters
mapDamage -i ancient_sample.bam -r reference.fasta \
    -y 0.5 \
    -m 30 \
    -b 15 \
    -t "Sample XYZ Damage Profile"
```

--------------------------------

### Python API - Alignment Functions

Source: https://context7.com/ginolhac/mapdamage/llms.txt

Functions for parsing CIGAR strings and extracting alignment coordinates from BAM reads.

```APIDOC
## Python API - Alignment Functions

### Description
This module provides functions for interpreting alignment information from BAM files, specifically focusing on parsing CIGAR strings and determining read coordinates.

### Functions
#### `get_coordinates(read: AlignedSegment) -> Tuple[int, int]`
Calculates the start and end coordinates of a read's alignment on the reference sequence.
- `read` (AlignedSegment): The aligned read object from `BAMReader`.
Returns: A tuple containing the start and end coordinates (0-based).

#### `get_around(read: AlignedSegment, ref_seq: str, window: int) -> str`
Extracts a DNA sequence window around the alignment of a read on the reference.
- `read` (AlignedSegment): The aligned read object.
- `ref_seq` (str): The reference sequence string.
- `window` (int): The size of the window to extract around the alignment.
Returns: The DNA sequence string within the specified window.

#### `align(read: AlignedSegment, ref_seq: str) -> str`
Aligns a read sequence to a reference sequence segment.
- `read` (AlignedSegment): The aligned read object.
- `ref_seq` (str): The reference sequence string.
Returns: The aligned sequence string.

#### `align_with_qual(read: AlignedSegment, ref_seq: str) -> Tuple[str, List[int]]`
Aligns a read sequence to a reference sequence segment and returns the corresponding quality scores.
- `read` (AlignedSegment): The aligned read object.
- `ref_seq` (str): The reference sequence string.
Returns: A tuple containing the aligned sequence string and a list of quality scores.

#### `parse_cigar(cigar: str) -> List[Tuple[int, str]]`
Parses a CIGAR string into a list of operations and their lengths.
- `cigar` (str): The CIGAR string from the BAM alignment.
Returns: A list of tuples, where each tuple is (length, operation_type).

### Request Example
```python
from mapdamage.align import get_coordinates, get_around, align, align_with_qual, parse_cigar
from pysam import AlignmentFile, AlignedSegment

# Assuming 'reader' is an initialized BAMReader object
bamfile = AlignmentFile(reader.filepath, "rb")
read = next(bamfile.fetch())

# Get 5' and 3' coordinates of aligned read
start_coord, end_coord = get_coordinates(read)
print(f"Coordinates: {start_coord}-{end_coord}")

# Get sequence around alignment (requires reference sequence)
# ref_seq = "..."
# window_seq = get_around(read, ref_seq, 50)

# Parse CIGAR string
cigar_ops = parse_cigar(read.cigarstring)
print(f"CIGAR Operations: {cigar_ops}")

bamfile.close()
```

### Response
#### Success Response (200)
- The alignment functions return processed data such as coordinates, sequence segments, quality scores, or parsed CIGAR operations.

#### Response Example
```json
{
  "read_coordinates": [1000, 1500],
  "parsed_cigar": [[50, "M"], [2, "D"]]
}
```
```

--------------------------------

### Rescale BAM Quality Scores (Python)

Source: https://context7.com/ginolhac/mapdamage/llms.txt

Rescales BAM quality scores based on estimated damage probabilities. This function requires a reference FASTA file and an options object containing input/output paths and rescaling parameters.

```python
from mapdamage.rescale import rescale_qual
import pysam
from pathlib import Path
from types import SimpleNamespace

# Prepare options object (normally from argument parsing)
options = SimpleNamespace(
    filename=Path("ancient_sample.bam"),
    folder=Path("results_ancient_sample"),
    rescale_out=Path("results_ancient_sample/ancient_sample.rescaled.bam"),
    rescale_length_5p=12,
    rescale_length_3p=12,
    seq_length=12
)

# Open reference FASTA
ref = pysam.FastaFile("reference.fasta")

# Perform rescaling
return_code = rescale_qual(ref, options)
# Returns 0 on success, 1 on error

ref.close()

# Rescaled reads get an MR (mapDamage Rescaled) tag
# with the expected number of rescaled bases
```

--------------------------------

### Rescale 5' and 3' Ends of BAM Files with mapDamage

Source: https://context7.com/ginolhac/mapdamage/llms.txt

Adjusts quality scores for the first 8 bases from both 5' and 3' ends of BAM files. This operation requires input BAM, a reference FASTA file, and specifies the sequence length for rescaling. The output is a modified BAM file with adjusted quality scores and an MR tag indicating the number of rescaled bases.

```bash
# Rescale only first 8 bases from 5' and 3' ends
mapDamage -i ancient_sample.bam -r reference.fasta --rescale \
    --seq-length 12 \
    --rescale-length-5p 8 \
    --rescale-length-3p 8
```

--------------------------------

### Alignment Functions for Parsing CIGAR Strings in Python

Source: https://context7.com/ginolhac/mapdamage/llms.txt

Offers Python functions for parsing CIGAR strings, extracting alignment coordinates, and retrieving sequence segments around alignments. These functions are crucial for detailed analysis of read alignments within the mapDamage framework.

```python
from mapdamage.align import get_coordinates, get_around, align, align_with_qual, parse_cigar

# Get 5' and 3' coordinates of aligned read

```

--------------------------------

### Quality Score Rescaling with mapDamage

Source: https://context7.com/ginolhac/mapdamage/llms.txt

Rescales base quality scores in BAM files to account for DNA damage, which can improve variant calling accuracy. This can be performed as part of a full analysis or as a standalone step if prior Bayesian estimation has been completed. The output path for the rescaled BAM file can be specified.

```bash
# Full analysis with quality rescaling
mapDamage -i ancient_sample.bam -r reference.fasta --rescale

# Rescale only (requires prior analysis with Bayesian estimation)
mapDamage --rescale-only -i ancient_sample.bam -r reference.fasta -d results_ancient_sample

# Specify output path for rescaled BAM
mapDamage -i ancient_sample.bam -r reference.fasta --rescale \
    --rescale-out rescaled_reads.bam
```

--------------------------------

### Plot Generation

Source: https://context7.com/ginolhac/mapdamage/llms.txt

Generates or regenerates damage pattern plots from existing analysis results. Allows customization of plot appearance and content.

```APIDOC
## POST /mapDamage/plot

### Description
Generates or regenerates DNA damage pattern plots from existing mapDamage analysis results stored in a specified directory. Customization options for plot appearance are available.

### Method
POST

### Endpoint
/mapDamage/plot

### Parameters
#### Query Parameters
- **--plot-only** (boolean) - Required - Flag to only generate plots without re-analyzing data.
- **-d** (string) - Required - Directory containing existing mapDamage analysis results.
- **-i** (string) - Optional - Input BAM file path (required if generating plots from scratch or if results directory is not specified).
- **-r** (string) - Optional - Reference FASTA file path (required if generating plots from scratch).
- **-y, --ymax** (float) - Optional - Maximum y-axis value for misincorporation plots. Defaults to 0.3.
- **-m, --readplot** (integer) - Optional - Number of bases to plot from read termini. Defaults to 25.
- **-b, --refplot** (integer) - Optional - Number of reference bases to plot around reads. Defaults to 10.
- **-t, --title** (string) - Optional - Custom title for the plots.

### Request Example
```bash
# Regenerate plots from existing results folder
mapDamage --plot-only -d results_ancient_sample

# Customize plot parameters
mapDamage -i ancient_sample.bam -r reference.fasta -y 0.5 -m 30 -b 15 -t "Sample XYZ Damage Profile"
```

### Response
#### Success Response (200)
- **plots** (array) - List of generated plot file paths.

#### Response Example
```json
{
  "plots": [
    "results_ancient_sample/misincorporation_plot.png",
    "results_ancient_sample/fragment_length_distribution.png"
  ]
}
```
```

--------------------------------

### Rescale Quality Scores

Source: https://context7.com/ginolhac/mapdamage/llms.txt

Adjusts quality scores of reads based on their position from the 5' and 3' ends. This is useful for correcting biases introduced by DNA damage.

```APIDOC
## POST /mapDamage/rescale

### Description
Rescales quality scores for reads in a BAM file, focusing on the first 8 bases from both the 5' and 3' ends. This process adjusts for potential DNA damage artifacts.

### Method
POST

### Endpoint
/mapDamage/rescale

### Parameters
#### Query Parameters
- **-i** (string) - Required - Input BAM file path.
- **-r** (string) - Required - Reference FASTA file path.
- **--rescale** (boolean) - Optional - Enables rescaling of quality scores.
- **--seq-length** (integer) - Optional - Specifies the sequence length to consider for rescaling. Defaults to 12.
- **--rescale-length-5p** (integer) - Optional - Number of bases from the 5' end to rescale. Defaults to 8.
- **--rescale-length-3p** (integer) - Optional - Number of bases from the 3' end to rescale. Defaults to 8.
- **-d** (string) - Optional - Output directory for results.

### Request Example
```bash
mapDamage -i ancient_sample.bam -r reference.fasta --rescale --seq-length 12 --rescale-length-5p 8 --rescale-length-3p 8
```

### Response
#### Success Response (200)
- **output_files** (array) - List of generated output files, including the rescaled BAM file and a file containing MR tags.

#### Response Example
```json
{
  "output_files": [
    "results_folder/ancient_sample.rescaled.bam",
    "results_folder/ancient_sample.rescaled.bam.mr.tag"
  ]
}
```
```

=== COMPLETE CONTENT === This response contains all available snippets from this library. No additional content exists. Do not make further requests.