repopack-output.txt

================================================================
Repopack Output File
================================================================

This file was generated by Repopack on: 2024-10-23T14:48:41.265Z

Purpose:
--------
This file contains a packed representation of the entire repository's contents.
It is designed to be easily consumable by AI systems for analysis, code review,
or other automated processes.

File Format:
------------
The content is organized as follows:
1. This header section
2. Repository structure
3. Multiple file entries, each consisting of:
  a. A separator line (================)
  b. The file path (File: path/to/file)
  c. Another separator line
  d. The full contents of the file
  e. A blank line

Usage Guidelines:
-----------------
1. This file should be treated as read-only. Any changes should be made to the
  original repository files, not this packed version.
2. When processing this file, use the separators and "File:" markers to
  distinguish between different files in the repository.
3. Be aware that this file may contain sensitive information. Handle it with
  the same level of security as you would the original repository.

Notes:
------
- Some files may have been excluded based on .gitignore rules and Repopack's
  configuration.
- Binary files are not included in this packed representation. Please refer to
  the Repository Structure section for a complete list of file paths, including
  binary files.



For more information about Repopack, visit: https://github.com/yamadashy/repopack

================================================================
Repository Structure
================================================================
.gitignore
.repopackignore
centromeres.json
Dockerfile
loh.py
README.md
run_analysis.sh

================================================================
Repository Files
================================================================

================
File: .gitignore
================
*.cns
*.csv
*.pileup

================
File: .repopackignore
================
*.txt
*.cns
*.csv
*.pileup
*.bed

================
File: centromeres.json
================
{
    "chr1": {"centromere": 123250000, "length": 249250621},
    "chr2": {"centromere": 93300000, "length": 243199373},
    "chr3": {"centromere": 91000000, "length": 198022430},
    "chr4": {"centromere": 50400000, "length": 191154276},
    "chr5": {"centromere": 48400000, "length": 180915260},
    "chr6": {"centromere": 61000000, "length": 171115067},
    "chr7": {"centromere": 59900000, "length": 159138663},
    "chr8": {"centromere": 45600000, "length": 146364022},
    "chr9": {"centromere": 49000000, "length": 141213431},
    "chr10": {"centromere": 40200000, "length": 135534747},
    "chr11": {"centromere": 53700000, "length": 135006516},
    "chr12": {"centromere": 35800000, "length": 133851895},
    "chr13": {"centromere": 17900000, "length": 115169878},
    "chr14": {"centromere": 17600000, "length": 107349540},
    "chr15": {"centromere": 19000000, "length": 102531392},
    "chr16": {"centromere": 36600000, "length": 90354753},
    "chr17": {"centromere": 24000000, "length": 81195210},
    "chr18": {"centromere": 17200000, "length": 78077248},
    "chr19": {"centromere": 26500000, "length": 59128983},
    "chr20": {"centromere": 27500000, "length": 63025520},
    "chr21": {"centromere": 13200000, "length": 48129895},
    "chr22": {"centromere": 14700000, "length": 51304566},
    "chrX": {"centromere": 60600000, "length": 155270560}
  }

================
File: Dockerfile
================
# Dockerfile
FROM debian:12-slim

# Prevent interactive prompts during installation
ENV DEBIAN_FRONTEND=noninteractive

# Install system dependencies
RUN apt-get update && apt-get install -y \
    python3 \
    python3-pip \
    python3-venv \
    samtools \
    openjdk-17-jre-headless \
    curl \
    --no-install-recommends \
    && rm -rf /var/lib/apt/lists/*

# Set up Python virtual environment
ENV VIRTUAL_ENV=/opt/venv
RUN python3 -m venv $VIRTUAL_ENV
ENV PATH="$VIRTUAL_ENV/bin:$PATH"

# Install Python dependencies in virtual environment
COPY requirements.txt /app/requirements.txt
RUN pip3 install --no-cache-dir -r /app/requirements.txt

# Create necessary directories
RUN mkdir -p /app/ref /app/samples /app/results /app/beds

# Copy files into container
COPY maf30_snps.txt /app/
COPY loh.py /app/
COPY centromeres.json /app/
COPY run_analysis.sh /app/
COPY beds/R210.bed /app/beds/

# Download VarScan from GitHub
RUN curl -L https://github.com/dkoboldt/varscan/releases/download/v2.4.6/VarScan.v2.4.6.jar -o /app/VarScan.v2.4.6.jar

# Set working directory
WORKDIR /app

# Make the run script executable
RUN chmod +x run_analysis.sh

# Default command
ENTRYPOINT ["./run_analysis.sh"]

================
File: loh.py
================
import sys
import json
import csv
from dataclasses import dataclass
from typing import List, Tuple, Dict, Optional, Set
from pathlib import Path
import logging
from enum import Enum

# Configure logging
logging.basicConfig(
    level=logging.INFO,
    format='%(asctime)s - %(levelname)s - %(message)s'
)

class Sex(Enum):
    MALE = 'Male'
    FEMALE = 'Female'
    UNKNOWN = 'Unknown'

@dataclass
class GenomicRegion:
    chromosome: str
    start: int
    end: int
    homozygous_count: int
    total_count: int
    
    @property
    def size(self) -> int:
        return self.end - self.start + 1
    
    @property
    def confidence(self) -> float:
        return (self.homozygous_count / self.total_count * 100) if self.total_count > 0 else 0

class LOHAnalyzer:
    def __init__(
        self,
        min_streak: int = 5,
        loh_threshold: float = 35.0,
        min_region_size: int = 1_000_000,
        max_gap: int = 2,
        sex_determination_threshold: float = 0.2
    ):
        self.min_streak = min_streak
        self.loh_threshold = loh_threshold
        self.min_region_size = min_region_size
        self.max_gap = max_gap
        self.sex_determination_threshold = sex_determination_threshold
        
    def _process_position(
        self, 
        chrom: str, 
        pos: int, 
        var_freq: float,
        current_region: Optional[GenomicRegion],
        gap_count: int,
        regions: List[GenomicRegion]
    ) -> Tuple[Optional[GenomicRegion], int]:
        """Process a single genomic position and update LOH regions."""
        is_homozygous = self._is_homozygous(var_freq)
        
        # Handle chromosome change
        if current_region and current_region.chromosome != chrom:
            if current_region:
                regions.append(current_region)
            current_region = None
            gap_count = 0
        
        # Process homozygous position
        if is_homozygous:
            if current_region is None:
                current_region = GenomicRegion(chrom, pos, pos, 1, 1)
            else:
                current_region = GenomicRegion(
                    chrom,
                    current_region.start,
                    pos,
                    current_region.homozygous_count + 1,
                    current_region.total_count + 1
                )
            gap_count = 0
        
        # Process heterozygous position
        else:
            if current_region is not None:
                gap_count += 1
                current_region = GenomicRegion(
                    chrom,
                    current_region.start,
                    pos,
                    current_region.homozygous_count,
                    current_region.total_count + 1
                )
                
                # Close region if gap is too large
                if gap_count > self.max_gap:
                    regions.append(current_region)
                    current_region = None
                    gap_count = 0
        
        return current_region, gap_count

    @staticmethod
    def load_centromeres(json_path: Path) -> Dict[str, int]:
        """Load centromere positions from JSON file."""
        try:
            with open(json_path) as f:
                data = json.load(f)
            return {k: v['centromere'] for k, v in data.items()}
        except Exception as e:
            logging.error(f"Failed to load centromeres file: {e}")
            raise

    @staticmethod
    def load_bed_regions(bed_path: Path) -> Dict[str, List[Tuple[int, int, str]]]:
        """Load gene regions from BED file."""
        regions = {}
        try:
            with open(bed_path) as f:
                for line in f:
                    chrom, start, end, gene = line.strip().split('\t')
                    chrom = f"chr{chrom}" if not chrom.startswith('chr') else chrom
                    regions.setdefault(chrom, []).append((int(start), int(end), gene))
            return regions
        except Exception as e:
            logging.error(f"Failed to load BED file: {e}")
            raise

    def _is_homozygous(self, var_freq: float) -> bool:
        """Determine if a variant frequency indicates homozygosity."""
        return var_freq <= self.loh_threshold or var_freq >= (100 - self.loh_threshold)

    def analyze_file(self, file_path: Path, centromeres: Dict[str, int]) -> Tuple[List[GenomicRegion], Sex]:
        """Analyze a CNS file for LOH regions and determine sample sex."""
        regions: List[GenomicRegion] = []
        current_region: Optional[GenomicRegion] = None
        gap_count = 0
        chrX_stats = {'het': 0, 'total': 0}
        
        try:
            with open(file_path) as f:
                next(f)  # Skip header
                for line in f:
                    chrom, pos, *_, var_freq = line.strip().split('\t')[:7]
                    pos = int(pos)
                    var_freq = float(var_freq.strip('%'))
                    
                    # Track chrX heterozygosity
                    if chrom == 'chrX':
                        chrX_stats['total'] += 1
                        if not self._is_homozygous(var_freq):
                            chrX_stats['het'] += 1
                    
                    # Process LOH regions
                    current_region, gap_count = self._process_position(
                        chrom, pos, var_freq, current_region, gap_count, regions
                    )
            
            # Add the last region if it exists
            if current_region:
                regions.append(current_region)
            
            # Determine sex
            sex = self._determine_sex(chrX_stats['het'], chrX_stats['total'])
            
            return self._filter_regions(regions, centromeres), sex
            
        except Exception as e:
            logging.error(f"Error processing file {file_path}: {e}")
            raise

    def _determine_sex(self, het_count: int, total_count: int) -> Sex:
        """Determine sample sex based on X chromosome heterozygosity."""
        if total_count == 0:
            return Sex.UNKNOWN
        ratio = het_count / total_count
        return Sex.FEMALE if ratio > self.sex_determination_threshold else Sex.MALE

    def _filter_regions(self, regions: List[GenomicRegion], centromeres: Dict[str, int]) -> List[GenomicRegion]:
        """Filter and split regions based on size and centromere positions."""
        filtered_regions = []
        for region in regions:
            if (region.homozygous_count >= self.min_streak and 
                region.size >= self.min_region_size):
                centromere_pos = centromeres.get(region.chromosome)
                if centromere_pos and region.start < centromere_pos < region.end:
                    # Split region at centromere
                    filtered_regions.extend([
                        GenomicRegion(
                            region.chromosome, region.start, centromere_pos - 1,
                            region.homozygous_count // 2, region.total_count // 2
                        ),
                        GenomicRegion(
                            region.chromosome, centromere_pos, region.end,
                            region.homozygous_count // 2, region.total_count // 2
                        )
                    ])
                else:
                    filtered_regions.append(region)
        return filtered_regions

    def find_affected_genes(self, region: GenomicRegion, bed_regions: Dict[str, List[Tuple[int, int, str]]]) -> Set[str]:
        """Find unique genes that overlap with a given genomic region."""
        affected_genes = set()  # Changed from list to set
        for start, end, gene in bed_regions.get(region.chromosome, []):
            if start <= region.end and end >= region.start:
                affected_genes.add(gene)  # Using add instead of append
        return affected_genes

class ResultsWriter:
    @staticmethod
    def save_to_csv(output_path: Path, results: List[Tuple[str, int, int, str]]):
        """Save analysis results to CSV file."""
        try:
            with open(output_path, 'w', newline='') as csvfile:
                writer = csv.DictWriter(csvfile, 
                    fieldnames=['Chromosome', 'Start', 'End', 'Affected_Genes'])
                writer.writeheader()
                for chrom, start, end, genes in results:
                    writer.writerow({
                        'Chromosome': chrom,
                        'Start': start,
                        'End': end,
                        'Affected_Genes': genes
                    })
        except Exception as e:
            logging.error(f"Failed to save results to {output_path}: {e}")
            raise

def main():
    if len(sys.argv) < 2 or len(sys.argv) > 3:
        logging.error("Usage: python floh.py <file_path> [<bed_file>]")
        sys.exit(1)

    try:
        analyzer = LOHAnalyzer()
        file_path = Path(sys.argv[1])
        bed_path = Path(sys.argv[2]) if len(sys.argv) == 3 else None
        centromeres_path = Path('centromeres.json')

        # Load required files
        centromeres = analyzer.load_centromeres(centromeres_path)
        bed_regions = analyzer.load_bed_regions(bed_path) if bed_path else None

        # Analyze file
        loh_regions, sex = analyzer.analyze_file(file_path, centromeres)

        # Process results
        results = []
        for region in loh_regions:
            if 'X' not in region.chromosome and region.homozygous_count > 40 and region.confidence > 90:
                if bed_regions:
                    affected_genes = analyzer.find_affected_genes(region, bed_regions)
                    if affected_genes:  # Only include regions with affected genes
                        results.append((
                            region.chromosome,
                            region.start,
                            region.end,
                            ','.join(sorted(affected_genes))  # Sort genes for consistent output
                        ))
                else:
                    results.append((
                        region.chromosome,
                        region.start,
                        region.end,
                        ''  # No genes information when BED file is not provided
                    ))

        # Save results
        output_path = file_path.with_suffix('.loh.csv')
        ResultsWriter.save_to_csv(output_path, results)

        # Log results
        logging.info(f"Analysis complete. Sex: {sex.value}")
        logging.info(f"Results saved to {output_path}")

    except Exception as e:
        logging.error(f"Analysis failed: {e}")
        sys.exit(1)

if __name__ == "__main__":
    main()

================
File: README.md
================
# siLOH - 🛢️🔎 SNP Inferred Loss of Heterozygosity Analysis Pipeline

siLOH is a Docker-based pipeline for analysing Loss of Heterozygosity (LOH) from NGS data. The pipeline integrates samtools, VarScan2, and a custom Python analysis script to identify potential regions of Loss of Heterozygosity (LOH).

The siLOH pipeline uses a targeted approach to identify candidate LOH regions by analysing the allele frequencies of predefined common SNPs (curated SNPs from dbSNP with Minor Allele Frequency (MAF) ≥ 30%) across the capture regions contained within an input BAM file.

## Pipeline Steps

1. **Pileup Generation**: Uses samtools mpileup to generate pileup data at specified positions
2. **Variant Calling**: VarScan2 analyzes the pileup to identify variants
3. **LOH Analysis**: Custom Python script processes variants to identify LOH regions
4. **Output Generation**: Results are saved in CSV format

## Prerequisites

- Docker installed and running
- Input BAM files
- Reference genome (ucsc_hg19.fa)
- BED file containing regions of interest
- Sufficient disk space for analysis

## Installation

1. Clone the repository:
```bash
git clone https://github.com/g-pyxl/siLOH.git
cd siLOH
```

2. Build the Docker image:
```bash
docker build -t siloh .
```

## Required Files

The pipeline expects the following files:

- **Reference Genome**: `ref/ucsc_hg19.fa` (GRCh37 example)
- **Sample BAM Files**: Contig naming should follow "chr1" format
- **BED File**: `beds/R210.bed` - The pipeline comes provided with an R210.bed designed for Lynch Syndrome calling. Tool can be ran with or without BED file calling.
- **maf30.txt**: Over 1.8m SNPs with high MAF, curated from dbSNP
- **centromeres.json**: Centromere positions relative to GRCh37 (provided in repo)

## Directory Structure

```
siLOH/
├── Dockerfile
├── requirements.txt
├── run_analysis.sh
├── loh.py
├── maf30.txt
├── centromeres.json
├── ref/
│   └── ucsc_hg19.fa
├── samples/
│   └── your_sample.bam
├── results/
└── beds/
    └── R210.bed
```

## Usage

1. Run the analysis:
```bash
docker run -v /path/to/ref:/app/ref \
           -v /path/to/samples:/app/samples \
           -v /path/to/results:/app/results \
           -v /path/to/beds:/app/beds \
           siloh your_sample_name
```

Replace `/path/to/` with your actual paths and `your_sample_name` with your BAM filename (without the .bam extension).

## Output

The pipeline generates the following files in the results directory:
- `{sample}.pileup`: Raw pileup data
- `{sample}.cns`: VarScan2 consensus output
- `{sample}.loh.csv`: Final LOH analysis results

The LOH CSV file contains:
- Chromosome
- Start position
- End position
- Affected genes (if BED file provided)

## Debugging

To enter the container without running the analysis:
```bash
docker run -it --entrypoint=/bin/bash siloh
```

## Resource Requirements

- Memory: Depends on input BAM file size
- Disk Space: ~3x the size of input BAM file
- CPU: Single-threaded processing

## Known Limitations

- Currently supports hg19 reference genome
- Single-sample processing
- Requires sorted, indexed BAM files - .bai should be within same directory as BAM

## Citation



## License



## Contact



## Contributing

================
File: run_analysis.sh
================
#!/bin/bash
# run_analysis.sh

# Exit on any error
set -e

# Check if input sample name was provided
if [ -z "$1" ]; then
    echo "Error: Please provide input sample name"
    echo "Usage: $0 <sample_name>"
    exit 1
fi

SAMPLE=$1
echo "Processing sample: $SAMPLE"

# Check if required files exist
if [ ! -f "ref/ucsc_hg19.fa" ]; then
    echo "Error: Reference file not found at ref/ucsc_hg19.fa"
    exit 1
fi

if [ ! -f "samples/${SAMPLE}.bam" ]; then
    echo "Error: BAM file not found at samples/${SAMPLE}.bam"
    exit 1
fi

if [ ! -f "beds/R210.bed" ]; then
    echo "Error: BED file not found at beds/R210.bed"
    exit 1
fi

echo "Running samtools mpileup..."
samtools mpileup -l maf30_snps.txt -f ref/ucsc_hg19.fa "samples/${SAMPLE}.bam" > "results/${SAMPLE}.pileup"

echo "Running VarScan..."
java -jar VarScan.v2.4.6.jar pileup2cns "results/${SAMPLE}.pileup" \
    --min-coverage 20 \
    --min-reads2 0 \
    --min-var-freq 0 > "results/${SAMPLE}.cns"

echo "Running LOH analysis..."
python3 loh.py "results/${SAMPLE}.cns" "beds/R210.bed"

echo "Analysis complete. Results can be found in results/${SAMPLE}.loh.csv"