Marat's Notes

Validating genomic data files, particularly VCF (Variant Call Format) files and other medical data formats, is crucial for ensuring data quality and reliability in bioinformatics workflows. This guide will walk you through the process of validating these files using SAM Tools and other essential tools.

Prerequisites

Before getting started, ensure you have:

System Requirements:
- Linux/Unix-based system
- Python 3.8+
- SAM Tools installed
- BCF Tools installed
- VCF Tools installed
Required Tools:
- SAM Tools
- BCF Tools
- VCF Tools
- Python packages for bioinformatics

Initial Setup

1. Tool Installation

Install the necessary tools:

# Install SAM Tools
sudo apt-get install samtools

# Install BCF Tools
sudo apt-get install bcftools

# Install VCF Tools
sudo apt-get install vcftools

# Install Python packages
pip install pyvcf pysam pandas numpy

2. Environment Setup

Set up your Python environment:

# Create virtual environment
python -m venv bioinfo-env
source bioinfo-env/bin/activate

# Install required packages
pip install pyvcf pysam pandas numpy

VCF File Validation

1. Basic VCF Validation

Validate a VCF file using BCF Tools:

import subprocess
from pathlib import Path

class VCFValidator:
    def __init__(self):
        self.bcftools = "bcftools"
        self.vcftools = "vcftools"
    
    def validate_vcf(self, vcf_file):
        try:
            # Check if file exists
            if not Path(vcf_file).exists():
                raise FileNotFoundError(f"VCF file not found: {vcf_file}")
            
            # Validate VCF file
            result = subprocess.run(
                [self.bcftools, "view", "-h", vcf_file],
                capture_output=True,
                text=True
            )
            
            if result.returncode != 0:
                raise ValueError(f"Invalid VCF file: {result.stderr}")
            
            return True
        except Exception as e:
            print(f"Error validating VCF file: {e}")
            return False

2. Advanced VCF Validation

Perform comprehensive VCF validation:

def validate_vcf_comprehensive(self, vcf_file):
    try:
        # Check file format
        format_check = subprocess.run(
            [self.bcftools, "view", "-h", vcf_file],
            capture_output=True,
            text=True
        )
        
        # Check for common errors
        error_check = subprocess.run(
            [self.bcftools, "view", "-e", vcf_file],
            capture_output=True,
            text=True
        )
        
        # Validate header
        header_check = subprocess.run(
            [self.bcftools, "view", "-h", vcf_file],
            capture_output=True,
            text=True
        )
        
        return {
            'format_valid': format_check.returncode == 0,
            'has_errors': error_check.returncode != 0,
            'header_valid': header_check.returncode == 0
        }
    except Exception as e:
        print(f"Error in comprehensive validation: {e}")
        return None

BAM/SAM File Validation

1. Basic BAM Validation

Validate BAM files using SAM Tools:

def validate_bam(self, bam_file):
    try:
        # Check if file exists
        if not Path(bam_file).exists():
            raise FileNotFoundError(f"BAM file not found: {bam_file}")
        
        # Validate BAM file
        result = subprocess.run(
            ["samtools", "quickcheck", bam_file],
            capture_output=True,
            text=True
        )
        
        if result.returncode != 0:
            raise ValueError(f"Invalid BAM file: {result.stderr}")
        
        return True
    except Exception as e:
        print(f"Error validating BAM file: {e}")
        return False

2. Advanced BAM Validation

Perform comprehensive BAM validation:

def validate_bam_comprehensive(self, bam_file):
    try:
        # Check file format
        format_check = subprocess.run(
            ["samtools", "view", "-H", bam_file],
            capture_output=True,
            text=True
        )
        
        # Check for common errors
        error_check = subprocess.run(
            ["samtools", "view", "-f", "4", bam_file],
            capture_output=True,
            text=True
        )
        
        # Validate index
        index_check = subprocess.run(
            ["samtools", "index", bam_file],
            capture_output=True,
            text=True
        )
        
        return {
            'format_valid': format_check.returncode == 0,
            'has_errors': error_check.returncode != 0,
            'index_valid': index_check.returncode == 0
        }
    except Exception as e:
        print(f"Error in comprehensive validation: {e}")
        return None

Medical File Validation

1. DICOM Validation

Validate DICOM files:

import pydicom
from pydicom.errors import InvalidDicomError

class MedicalFileValidator:
    def validate_dicom(self, dicom_file):
        try:
            # Check if file exists
            if not Path(dicom_file).exists():
                raise FileNotFoundError(f"DICOM file not found: {dicom_file}")
            
            # Validate DICOM file
            ds = pydicom.dcmread(dicom_file)
            
            # Check required fields
            required_fields = ['PatientName', 'StudyDate', 'Modality']
            missing_fields = [field for field in required_fields if field not in ds]
            
            if missing_fields:
                raise ValueError(f"Missing required fields: {missing_fields}")
            
            return True
        except InvalidDicomError as e:
            print(f"Invalid DICOM file: {e}")
            return False
        except Exception as e:
            print(f"Error validating DICOM file: {e}")
            return False

2. FASTQ Validation

Validate FASTQ files:

def validate_fastq(self, fastq_file):
    try:
        # Check if file exists
        if not Path(fastq_file).exists():
            raise FileNotFoundError(f"FASTQ file not found: {fastq_file}")
        
        # Validate FASTQ file
        with open(fastq_file, 'r') as f:
            line_count = 0
            for line in f:
                line_count += 1
                if line_count % 4 == 1 and not line.startswith('@'):
                    raise ValueError(f"Invalid FASTQ format at line {line_count}")
                elif line_count % 4 == 3 and not line.startswith('+'):
                    raise ValueError(f"Invalid FASTQ format at line {line_count}")
        
        return True
    except Exception as e:
        print(f"Error validating FASTQ file: {e}")
        return False

Quality Control

1. VCF Quality Metrics

Calculate VCF quality metrics:

def calculate_vcf_metrics(self, vcf_file):
    try:
        # Calculate basic statistics
        stats = subprocess.run(
            [self.vcftools, "--vcf", vcf_file, "--freq"],
            capture_output=True,
            text=True
        )
        
        # Calculate quality scores
        quality = subprocess.run(
            [self.bcftools, "query", "-f", "%QUAL\n", vcf_file],
            capture_output=True,
            text=True
        )
        
        return {
            'stats': stats.stdout,
            'quality_scores': quality.stdout
        }
    except Exception as e:
        print(f"Error calculating VCF metrics: {e}")
        return None

2. BAM Quality Metrics

Calculate BAM quality metrics:

def calculate_bam_metrics(self, bam_file):
    try:
        # Calculate basic statistics
        stats = subprocess.run(
            ["samtools", "stats", bam_file],
            capture_output=True,
            text=True
        )
        
        # Calculate coverage
        coverage = subprocess.run(
            ["samtools", "depth", bam_file],
            capture_output=True,
            text=True
        )
        
        return {
            'stats': stats.stdout,
            'coverage': coverage.stdout
        }
    except Exception as e:
        print(f"Error calculating BAM metrics: {e}")
        return None

Error Handling

1. Validation Error Handling

Implement robust error handling:

class ValidationErrorHandler:
    def handle_validation_error(self, error, file_type):
        try:
            error_log = {
                'timestamp': datetime.now().isoformat(),
                'file_type': file_type,
                'error': str(error),
                'severity': 'ERROR'
            }
            
            # Log error
            self.log_error(error_log)
            
            # Notify if critical
            if self.is_critical_error(error):
                self.notify_admin(error_log)
            
            return error_log
        except Exception as e:
            print(f"Error handling validation error: {e}")
            return None

2. Error Logging

Set up comprehensive error logging:

def log_error(self, error_log):
    try:
        # Log to file
        with open('validation_errors.log', 'a') as f:
            f.write(json.dumps(error_log) + '\n')
        
        # Log to database if available
        if hasattr(self, 'db_connection'):
            self.db_connection.insert_error(error_log)
    except Exception as e:
        print(f"Error logging error: {e}")

Best Practices

1. Validation Workflow

Implement a systematic validation workflow:

def validate_workflow(self, input_files):
    try:
        results = {}
        
        for file in input_files:
            file_type = self.detect_file_type(file)
            
            if file_type == 'vcf':
                results[file] = self.validate_vcf_comprehensive(file)
            elif file_type == 'bam':
                results[file] = self.validate_bam_comprehensive(file)
            elif file_type == 'dicom':
                results[file] = self.validate_dicom(file)
            elif file_type == 'fastq':
                results[file] = self.validate_fastq(file)
        
        return results
    except Exception as e:
        print(f"Error in validation workflow: {e}")
        return None

2. Performance Optimization

Optimize validation performance:

def optimize_validation(self, input_files):
    try:
        # Use multiprocessing for parallel validation
        with Pool() as pool:
            results = pool.map(self.validate_file, input_files)
        
        return dict(zip(input_files, results))
    except Exception as e:
        print(f"Error optimizing validation: {e}")
        return None

Conclusion

Proper validation of genomic and medical files is crucial for ensuring data quality and reliability. By following this guide, you can:

Validate VCF files using BCF Tools
Validate BAM/SAM files using SAM Tools
Validate medical files (DICOM, FASTQ)
Implement quality control measures
Handle errors effectively

Remember to:

Regularly validate your data files
Implement comprehensive error handling
Follow best practices for file validation
Optimize validation performance
Keep track of validation results

With proper implementation, you can ensure the quality and reliability of your genomic and medical data files.

VCF and Medical Files Validation: A Comprehensive Guide Using SAM Tools