000 Annotation File Processing

In [38]:

# Annotation File Processing
# This notebook processes GENCODE GTF annotation files to extract protein-coding
# transcript information and exports it in a tabular format suitable for downstream analysis.

import polars as pl


# =============================================================================
# Configuration
# =============================================================================
REFERENCE_DIR = "/data/for_paper/data/reference"
GENCODE_VERSION = "v47"

# Input: GENCODE GTF files (downloaded from https://www.gencodegenes.org/)
GTF_FILES = {
    "hg38": f"{REFERENCE_DIR}/gencode.{GENCODE_VERSION}.basic.annotation.gtf.gz",
    "hg19": f"{REFERENCE_DIR}/gencode.{GENCODE_VERSION}lift37.basic.annotation.gtf.gz",
}

# Output: Processed annotation TSV files
OUTPUT_FILES = {
    "hg38": f"{REFERENCE_DIR}/gencode.{GENCODE_VERSION}.basic.annotation.processed.tsv",
    "hg19": f"{REFERENCE_DIR}/gencode.{GENCODE_VERSION}lift37.basic.annotation.processed.tsv",
}

print(f"GENCODE version: {GENCODE_VERSION}")
print(f"Reference directory: {REFERENCE_DIR}")

# Annotation File Processing # This notebook processes GENCODE GTF annotation files to extract protein-coding # transcript information and exports it in a tabular format suitable for downstream analysis. import polars as pl # ============================================================================= # Configuration # ============================================================================= REFERENCE_DIR = "/data/for_paper/data/reference" GENCODE_VERSION = "v47" # Input: GENCODE GTF files (downloaded from https://www.gencodegenes.org/) GTF_FILES = { "hg38": f"{REFERENCE_DIR}/gencode.{GENCODE_VERSION}.basic.annotation.gtf.gz", "hg19": f"{REFERENCE_DIR}/gencode.{GENCODE_VERSION}lift37.basic.annotation.gtf.gz", } # Output: Processed annotation TSV files OUTPUT_FILES = { "hg38": f"{REFERENCE_DIR}/gencode.{GENCODE_VERSION}.basic.annotation.processed.tsv", "hg19": f"{REFERENCE_DIR}/gencode.{GENCODE_VERSION}lift37.basic.annotation.processed.tsv", } print(f"GENCODE version: {GENCODE_VERSION}") print(f"Reference directory: {REFERENCE_DIR}")

GENCODE version: v47
Reference directory: /data/for_paper/data/reference

In [39]:

# Helper functions to validate and adjust CDS (coding sequence) boundaries
# These functions ensure that the exon coordinates properly include the stop codon (3 bp)


def check_start_alignment(row):
    """
    Validates that CDS start aligns with exon start and adjusts for stop codon on minus strand.
    For minus strand genes, the stop codon is at the 3' end (lowest genomic position), so we subtract 3 bp.
    """
    cds_start = row["cds_start"]
    exon_starts = list(map(int, row["exon_starts"].strip(",").split(",")))
    if row["strand"] == "-":
        # Extend first exon by 3 bp to include stop codon (on minus strand)
        exon_starts[0] -= 3
    assert cds_start == exon_starts[0], f"{cds_start} != {exon_starts[0]} {row['transcript_id']}"

    exon_starts = ",".join(map(str, exon_starts)) + ","
    return exon_starts


def check_end_alignment(row):
    """
    Validates that CDS end aligns with exon end and adjusts for stop codon on plus strand.
    For plus strand genes, the stop codon is at the 3' end (highest genomic position), so we add 3 bp.
    """
    cds_end = row["cds_end"]
    exon_ends = list(map(int, row["exon_ends"].strip(",").split(",")))
    if row["strand"] == "+":
        # Extend last exon by 3 bp to include stop codon (on plus strand)
        exon_ends[-1] += 3
    assert cds_end == exon_ends[-1], f"{cds_end} != {exon_ends[-1]} {row['transcript_id']}"

    exon_ends = ",".join(map(str, exon_ends)) + ","
    return exon_ends

# Helper functions to validate and adjust CDS (coding sequence) boundaries # These functions ensure that the exon coordinates properly include the stop codon (3 bp) def check_start_alignment(row): """ Validates that CDS start aligns with exon start and adjusts for stop codon on minus strand. For minus strand genes, the stop codon is at the 3' end (lowest genomic position), so we subtract 3 bp. """ cds_start = row["cds_start"] exon_starts = list(map(int, row["exon_starts"].strip(",").split(","))) if row["strand"] == "-": # Extend first exon by 3 bp to include stop codon (on minus strand) exon_starts[0] -= 3 assert cds_start == exon_starts[0], f"{cds_start} != {exon_starts[0]} {row['transcript_id']}" exon_starts = ",".join(map(str, exon_starts)) + "," return exon_starts def check_end_alignment(row): """ Validates that CDS end aligns with exon end and adjusts for stop codon on plus strand. For plus strand genes, the stop codon is at the 3' end (highest genomic position), so we add 3 bp. """ cds_end = row["cds_end"] exon_ends = list(map(int, row["exon_ends"].strip(",").split(","))) if row["strand"] == "+": # Extend last exon by 3 bp to include stop codon (on plus strand) exon_ends[-1] += 3 assert cds_end == exon_ends[-1], f"{cds_end} != {exon_ends[-1]} {row['transcript_id']}" exon_ends = ",".join(map(str, exon_ends)) + "," return exon_ends

In [40]:

def process_gtf_file(gtf_file, output_file):
    """
    Process a GENCODE GTF file and extract protein-coding transcript annotations.

    This function:
    1. Parses the GTF file and extracts relevant attributes
    2. Filters for protein-coding genes only
    3. Aggregates exon/CDS coordinates per transcript
    4. Adjusts coordinates to include stop codons
    5. Outputs a tab-separated file with transcript annotations

    Args:
        gtf_file: Path to input GENCODE GTF file (can be gzipped)
        output_file: Path for output TSV file
    """
    # Read GTF file (standard 9-column format)
    gtf = pl.read_csv(gtf_file, comment_prefix="#", separator="\t", has_header=False)
    gtf.columns = ["chrom", "source", "feature", "start", "end", "score", "strand", "frame", "attribute"]

    # Parse the attribute column (column 9) to extract key-value pairs
    # GTF attributes are semicolon-separated with format: key "value"
    gtf = gtf.with_columns(
        [
            pl.col("attribute").str.extract('gene_id "(.*?)"', 1).alias("gene_id"),
            pl.col("attribute").str.extract('transcript_id "(.*?)"', 1).alias("transcript_id"),
            pl.col("attribute").str.extract('gene_name "(.*?)"', 1).alias("gene_name"),
            pl.col("attribute").str.extract('gene_type "(.*?)"', 1).alias("gene_type"),
            pl.col("attribute").str.extract('transcript_type "(.*?)"', 1).alias("transcript_type"),
            pl.col("attribute").str.extract("exon_number (.*?);", 1).alias("exon_number"),
        ]
    )

    # Flag canonical transcripts (Ensembl canonical and MANE Select)
    gtf = gtf.with_columns(pl.col("attribute").str.contains("Ensembl_canonical").alias("is_canonical"))
    gtf = gtf.with_columns(pl.col("attribute").str.contains("MANE_Select").alias("is_mane_select"))

    # Filter to protein-coding genes only, exclude gene-level features
    protein_coding_gtf = gtf.filter((pl.col("gene_type") == "protein_coding") & (pl.col("feature") != "gene"))
    # Filter to protein-coding transcripts only
    protein_coding_gtf = protein_coding_gtf.filter(pl.col("transcript_type") == "protein_coding")

    # Convert from 1-based (GTF) to 0-based coordinates (BED-like format)
    protein_coding_gtf = protein_coding_gtf.with_columns(pl.col("start") - 1)

    # Aggregate CDS exon coordinates per transcript
    # Creates comma-separated lists of exon start/end positions (sorted by genomic position)
    exon_starts = (
        protein_coding_gtf.filter(pl.col("feature") == "CDS")
        .group_by("transcript_id")
        .agg(
            (pl.col("start").sort().cast(str).str.join(",") + ",").alias("exon_starts"),
            (pl.col("end").sort().cast(str).str.join(",") + ",").alias("exon_ends"),
            pl.col("exon_number").max().alias("exon_count"),
        )
    )

    # Calculate CDS boundaries with stop codon adjustment
    # GENCODE GTF excludes stop codon from CDS, but we want to include it
    # For + strand: stop codon is after the last CDS position (add 3 to max end)
    # For - strand: stop codon is before the first CDS position (subtract 3 from min start)
    # Note: Using min()/max() instead of first()/last() to avoid dependency on row order
    cds_starts = (
        protein_coding_gtf.filter(pl.col("feature") == "CDS")
        .group_by("transcript_id")
        .agg(
            pl.when(pl.col("strand").first() == "-")
            .then(pl.col("start").min() - 3)  # Include stop codon at 3' end (lowest genomic position)
            .otherwise(pl.col("start").min())  # 5' end, no adjustment needed
            .alias("cds_start"),
            pl.when(pl.col("strand").first() == "-")
            .then(pl.col("end").max())  # 5' end (highest genomic position), no adjustment
            .otherwise(pl.col("end").max() + 3)  # Include stop codon at 3' end
            .alias("cds_end"),
        )
    )

    # Get transcript-level metadata (gene info, coordinates, canonical status)
    tx_starts = (
        protein_coding_gtf.filter(pl.col("feature") == "transcript")
        .group_by("transcript_id")
        .agg(
            pl.col("gene_id").first().alias("gene_id"),
            pl.col("gene_name").first().alias("gene_name"),
            pl.col("chrom").first().alias("chrom"),
            pl.col("strand").first().alias("strand"),
            pl.col("start").min().alias("tx_start"),
            pl.col("end").max().alias("tx_end"),
            pl.col("transcript_type").first().alias("transcript_type"),
            pl.col("is_canonical").first().alias("is_canonical"),
            pl.col("is_mane_select").first().alias("is_mane_select"),
        )
    )

    # Join all transcript information together
    joined_df = tx_starts.join(cds_starts, on="transcript_id", how="inner").join(
        exon_starts, on="transcript_id", how="inner"
    )

    # Validate and adjust exon coordinates to include stop codon
    joined_df = joined_df.with_columns(
        pl.struct(["cds_start", "exon_starts", "strand", "transcript_id"])
        .map_elements(check_start_alignment, return_dtype=pl.Utf8)
        .alias("exon_starts"),
        pl.struct(["cds_end", "exon_ends", "strand", "transcript_id"])
        .map_elements(check_end_alignment, return_dtype=pl.Utf8)
        .alias("exon_ends"),
    )

    # Sort by chromosome and position, then select and rename columns for output
    joined_df = joined_df.sort(["chrom", "tx_start"])
    joined_df = joined_df.select(
        [
            "gene_id",
            "transcript_id",
            "chrom",
            "strand",
            "tx_start",
            "tx_end",
            "cds_start",
            "cds_end",
            "exon_count",
            "exon_starts",
            "exon_ends",
            "gene_name",
            "transcript_type",
            "is_canonical",
            "is_mane_select",
        ]
    ).rename(
        {
            "transcript_id": "name",  # Transcript ID becomes the 'name' field
            "tx_start": "txStart",  # Transcript start position
            "tx_end": "txEnd",  # Transcript end position
            "cds_start": "cdsStart",  # CDS start (including stop codon adjustment)
            "cds_end": "cdsEnd",  # CDS end (including stop codon adjustment)
            "exon_starts": "exonStarts",  # Comma-separated exon start positions
            "exon_ends": "exonEnds",  # Comma-separated exon end positions
        }
    )

    # Write output as tab-separated file
    joined_df.write_csv(output_file, separator="\t")

    return joined_df

def process_gtf_file(gtf_file, output_file): """ Process a GENCODE GTF file and extract protein-coding transcript annotations. This function: 1. Parses the GTF file and extracts relevant attributes 2. Filters for protein-coding genes only 3. Aggregates exon/CDS coordinates per transcript 4. Adjusts coordinates to include stop codons 5. Outputs a tab-separated file with transcript annotations Args: gtf_file: Path to input GENCODE GTF file (can be gzipped) output_file: Path for output TSV file """ # Read GTF file (standard 9-column format) gtf = pl.read_csv(gtf_file, comment_prefix="#", separator="\t", has_header=False) gtf.columns = ["chrom", "source", "feature", "start", "end", "score", "strand", "frame", "attribute"] # Parse the attribute column (column 9) to extract key-value pairs # GTF attributes are semicolon-separated with format: key "value" gtf = gtf.with_columns( [ pl.col("attribute").str.extract('gene_id "(.*?)"', 1).alias("gene_id"), pl.col("attribute").str.extract('transcript_id "(.*?)"', 1).alias("transcript_id"), pl.col("attribute").str.extract('gene_name "(.*?)"', 1).alias("gene_name"), pl.col("attribute").str.extract('gene_type "(.*?)"', 1).alias("gene_type"), pl.col("attribute").str.extract('transcript_type "(.*?)"', 1).alias("transcript_type"), pl.col("attribute").str.extract("exon_number (.*?);", 1).alias("exon_number"), ] ) # Flag canonical transcripts (Ensembl canonical and MANE Select) gtf = gtf.with_columns(pl.col("attribute").str.contains("Ensembl_canonical").alias("is_canonical")) gtf = gtf.with_columns(pl.col("attribute").str.contains("MANE_Select").alias("is_mane_select")) # Filter to protein-coding genes only, exclude gene-level features protein_coding_gtf = gtf.filter((pl.col("gene_type") == "protein_coding") & (pl.col("feature") != "gene")) # Filter to protein-coding transcripts only protein_coding_gtf = protein_coding_gtf.filter(pl.col("transcript_type") == "protein_coding") # Convert from 1-based (GTF) to 0-based coordinates (BED-like format) protein_coding_gtf = protein_coding_gtf.with_columns(pl.col("start") - 1) # Aggregate CDS exon coordinates per transcript # Creates comma-separated lists of exon start/end positions (sorted by genomic position) exon_starts = ( protein_coding_gtf.filter(pl.col("feature") == "CDS") .group_by("transcript_id") .agg( (pl.col("start").sort().cast(str).str.join(",") + ",").alias("exon_starts"), (pl.col("end").sort().cast(str).str.join(",") + ",").alias("exon_ends"), pl.col("exon_number").max().alias("exon_count"), ) ) # Calculate CDS boundaries with stop codon adjustment # GENCODE GTF excludes stop codon from CDS, but we want to include it # For + strand: stop codon is after the last CDS position (add 3 to max end) # For - strand: stop codon is before the first CDS position (subtract 3 from min start) # Note: Using min()/max() instead of first()/last() to avoid dependency on row order cds_starts = ( protein_coding_gtf.filter(pl.col("feature") == "CDS") .group_by("transcript_id") .agg( pl.when(pl.col("strand").first() == "-") .then(pl.col("start").min() - 3) # Include stop codon at 3' end (lowest genomic position) .otherwise(pl.col("start").min()) # 5' end, no adjustment needed .alias("cds_start"), pl.when(pl.col("strand").first() == "-") .then(pl.col("end").max()) # 5' end (highest genomic position), no adjustment .otherwise(pl.col("end").max() + 3) # Include stop codon at 3' end .alias("cds_end"), ) ) # Get transcript-level metadata (gene info, coordinates, canonical status) tx_starts = ( protein_coding_gtf.filter(pl.col("feature") == "transcript") .group_by("transcript_id") .agg( pl.col("gene_id").first().alias("gene_id"), pl.col("gene_name").first().alias("gene_name"), pl.col("chrom").first().alias("chrom"), pl.col("strand").first().alias("strand"), pl.col("start").min().alias("tx_start"), pl.col("end").max().alias("tx_end"), pl.col("transcript_type").first().alias("transcript_type"), pl.col("is_canonical").first().alias("is_canonical"), pl.col("is_mane_select").first().alias("is_mane_select"), ) ) # Join all transcript information together joined_df = tx_starts.join(cds_starts, on="transcript_id", how="inner").join( exon_starts, on="transcript_id", how="inner" ) # Validate and adjust exon coordinates to include stop codon joined_df = joined_df.with_columns( pl.struct(["cds_start", "exon_starts", "strand", "transcript_id"]) .map_elements(check_start_alignment, return_dtype=pl.Utf8) .alias("exon_starts"), pl.struct(["cds_end", "exon_ends", "strand", "transcript_id"]) .map_elements(check_end_alignment, return_dtype=pl.Utf8) .alias("exon_ends"), ) # Sort by chromosome and position, then select and rename columns for output joined_df = joined_df.sort(["chrom", "tx_start"]) joined_df = joined_df.select( [ "gene_id", "transcript_id", "chrom", "strand", "tx_start", "tx_end", "cds_start", "cds_end", "exon_count", "exon_starts", "exon_ends", "gene_name", "transcript_type", "is_canonical", "is_mane_select", ] ).rename( { "transcript_id": "name", # Transcript ID becomes the 'name' field "tx_start": "txStart", # Transcript start position "tx_end": "txEnd", # Transcript end position "cds_start": "cdsStart", # CDS start (including stop codon adjustment) "cds_end": "cdsEnd", # CDS end (including stop codon adjustment) "exon_starts": "exonStarts", # Comma-separated exon start positions "exon_ends": "exonEnds", # Comma-separated exon end positions } ) # Write output as tab-separated file joined_df.write_csv(output_file, separator="\t") return joined_df

In [41]:

# Process GENCODE annotation files for both GRCh38 (hg38) and GRCh37 (hg19) assemblies

for assembly in ["hg38", "hg19"]:
    gtf_file = GTF_FILES[assembly]
    output_file = OUTPUT_FILES[assembly]

    print(f"Processing {assembly}: {gtf_file}")
    _ = process_gtf_file(gtf_file, output_file)
    print(f"Output saved to: {output_file}\n")

# Process GENCODE annotation files for both GRCh38 (hg38) and GRCh37 (hg19) assemblies for assembly in ["hg38", "hg19"]: gtf_file = GTF_FILES[assembly] output_file = OUTPUT_FILES[assembly] print(f"Processing {assembly}: {gtf_file}") _ = process_gtf_file(gtf_file, output_file) print(f"Output saved to: {output_file}\n")

Processing hg38: /data/for_paper/data/reference/gencode.v47.basic.annotation.gtf.gz
Output saved to: /data/for_paper/data/reference/gencode.v47.basic.annotation.processed.tsv

Processing hg19: /data/for_paper/data/reference/gencode.v47lift37.basic.annotation.gtf.gz
Output saved to: /data/for_paper/data/reference/gencode.v47lift37.basic.annotation.processed.tsv

---

Part 1: Quality Control of CDS and Extract CDS Sequence¶

This section extracts CDS sequences from the reference genome and validates them for downstream variant analysis.

In [42]:

# =============================================================================
# Additional Configuration for CDS Extraction
# =============================================================================

# Reference genome (GRCh38/hg38)
REFERENCE_GENOME = f"{REFERENCE_DIR}/hg38/hg38.fa"

# Input: Processed annotation from above
ANNOTATION_FILE = OUTPUT_FILES["hg38"]

# Output: Filtered transcripts with CDS sequences
FILTERED_TRANSCRIPTS_FILE = f"{REFERENCE_DIR}/gencode.{GENCODE_VERSION}.basic.annotation.processed.filtered.tsv"

# Valid chromosomes for analysis
VALID_CHROMS = [f"chr{i}" for i in range(1, 23)] + ["chrX"]

# DNA complement mapping for reverse complement operations
COMPLEMENT = {"A": "T", "T": "A", "G": "C", "C": "G", "N": "N"}

print(f"Reference genome: {REFERENCE_GENOME}")
print(f"Annotation file: {ANNOTATION_FILE}")
print(f"Output file: {FILTERED_TRANSCRIPTS_FILE}")

# ============================================================================= # Additional Configuration for CDS Extraction # ============================================================================= # Reference genome (GRCh38/hg38) REFERENCE_GENOME = f"{REFERENCE_DIR}/hg38/hg38.fa" # Input: Processed annotation from above ANNOTATION_FILE = OUTPUT_FILES["hg38"] # Output: Filtered transcripts with CDS sequences FILTERED_TRANSCRIPTS_FILE = f"{REFERENCE_DIR}/gencode.{GENCODE_VERSION}.basic.annotation.processed.filtered.tsv" # Valid chromosomes for analysis VALID_CHROMS = [f"chr{i}" for i in range(1, 23)] + ["chrX"] # DNA complement mapping for reverse complement operations COMPLEMENT = {"A": "T", "T": "A", "G": "C", "C": "G", "N": "N"} print(f"Reference genome: {REFERENCE_GENOME}") print(f"Annotation file: {ANNOTATION_FILE}") print(f"Output file: {FILTERED_TRANSCRIPTS_FILE}")

Reference genome: /data/for_paper/data/reference/hg38/hg38.fa
Annotation file: /data/for_paper/data/reference/gencode.v47.basic.annotation.processed.tsv
Output file: /data/for_paper/data/reference/gencode.v47.basic.annotation.processed.filtered.tsv

In [43]:

# =============================================================================
# CDS Extraction and Quality Control Functions
# =============================================================================


def extract_cds_sequence(row, fasta):
    """
    Extract the coding sequence (CDS) for a transcript from the reference genome.

    This function:
    1. Iterates through exons and extracts only the CDS-overlapping portions
    2. Concatenates exon sequences in genomic order
    3. Reverse complements for minus strand genes

    Args:
        row: DataFrame row containing transcript annotation (chrom, strand, cdsStart, cdsEnd, exonStarts, exonEnds)
        fasta: Dictionary mapping chromosome names to their sequences

    Returns:
        str: The complete CDS sequence in 5' to 3' orientation (transcript strand)
    """
    chrom = row["chrom"]
    strand = row["strand"]
    cds_start = row["cdsStart"]
    cds_end = row["cdsEnd"]

    # Parse comma-separated exon coordinates from annotation file
    exon_starts = [int(x) for x in row["exonStarts"].rstrip(",").split(",")]
    exon_ends = [int(x) for x in row["exonEnds"].rstrip(",").split(",")]

    # Ensure exon boundaries encompass the full CDS (handles edge cases)
    if exon_starts[0] > cds_start:
        exon_starts[0] = cds_start
    if exon_ends[-1] < cds_end:
        exon_ends[-1] = cds_end

    # Extract CDS sequence by iterating through exons
    cds_sequence = ""

    for start, end in zip(exon_starts, exon_ends):
        # Find overlap between this exon and the CDS region
        overlap_start = max(start, cds_start)
        overlap_end = min(end, cds_end)

        if overlap_start < overlap_end:
            # Extract sequence from this exon segment (0-based coordinates)
            seq = str(fasta[chrom][overlap_start:overlap_end]).upper()
            cds_sequence += seq

    # For minus strand genes, reverse complement to get 5' to 3' orientation
    if strand == "-":
        cds_sequence = "".join(COMPLEMENT[base] for base in cds_sequence[::-1])

    return cds_sequence


def check_cds_quality(sequence):
    """
    Validate CDS sequence quality for downstream variant analysis.

    Quality criteria checked:
    1. Starts with ATG (methionine start codon)
    2. Ends with a stop codon (TAA, TAG, or TGA)
    3. Length is divisible by 3 (complete codons)
    4. No premature stop codons within the coding region

    Args:
        sequence: CDS nucleotide sequence string

    Returns:
        dict: Quality metrics including boolean flags and sequence length
    """
    if not sequence or len(sequence) < 3:
        return {
            "has_start_codon": False,
            "has_stop_codon": False,
            "length_divisible_by_3": False,
            "has_internal_stop_codons": False,
            "length": len(sequence) if sequence else 0,
        }

    # Check for canonical start codon (ATG = Methionine)
    has_start_codon = sequence[:3] == "ATG"

    # Check for stop codon at the end
    has_stop_codon = sequence[-3:] in ["TAA", "TAG", "TGA"]

    # CDS should be in-frame (length divisible by 3)
    length_divisible_by_3 = len(sequence) % 3 == 0

    # Check for internal stop codons (premature termination)
    # These indicate potential annotation errors or pseudogenes
    has_internal_stop_codons = False
    if len(sequence) >= 6:  # Need at least 2 codons to check for internal stops
        # Check all codons except the last one (which should be a stop)
        for i in range(0, len(sequence) - 3, 3):
            codon = sequence[i : i + 3]
            if codon in ["TAA", "TAG", "TGA"]:
                has_internal_stop_codons = True
                break

    return {
        "has_start_codon": has_start_codon,
        "has_stop_codon": has_stop_codon,
        "length_divisible_by_3": length_divisible_by_3,
        "has_internal_stop_codons": has_internal_stop_codons,
        "length": len(sequence),
    }

# ============================================================================= # CDS Extraction and Quality Control Functions # ============================================================================= def extract_cds_sequence(row, fasta): """ Extract the coding sequence (CDS) for a transcript from the reference genome. This function: 1. Iterates through exons and extracts only the CDS-overlapping portions 2. Concatenates exon sequences in genomic order 3. Reverse complements for minus strand genes Args: row: DataFrame row containing transcript annotation (chrom, strand, cdsStart, cdsEnd, exonStarts, exonEnds) fasta: Dictionary mapping chromosome names to their sequences Returns: str: The complete CDS sequence in 5' to 3' orientation (transcript strand) """ chrom = row["chrom"] strand = row["strand"] cds_start = row["cdsStart"] cds_end = row["cdsEnd"] # Parse comma-separated exon coordinates from annotation file exon_starts = [int(x) for x in row["exonStarts"].rstrip(",").split(",")] exon_ends = [int(x) for x in row["exonEnds"].rstrip(",").split(",")] # Ensure exon boundaries encompass the full CDS (handles edge cases) if exon_starts[0] > cds_start: exon_starts[0] = cds_start if exon_ends[-1] < cds_end: exon_ends[-1] = cds_end # Extract CDS sequence by iterating through exons cds_sequence = "" for start, end in zip(exon_starts, exon_ends): # Find overlap between this exon and the CDS region overlap_start = max(start, cds_start) overlap_end = min(end, cds_end) if overlap_start < overlap_end: # Extract sequence from this exon segment (0-based coordinates) seq = str(fasta[chrom][overlap_start:overlap_end]).upper() cds_sequence += seq # For minus strand genes, reverse complement to get 5' to 3' orientation if strand == "-": cds_sequence = "".join(COMPLEMENT[base] for base in cds_sequence[::-1]) return cds_sequence def check_cds_quality(sequence): """ Validate CDS sequence quality for downstream variant analysis. Quality criteria checked: 1. Starts with ATG (methionine start codon) 2. Ends with a stop codon (TAA, TAG, or TGA) 3. Length is divisible by 3 (complete codons) 4. No premature stop codons within the coding region Args: sequence: CDS nucleotide sequence string Returns: dict: Quality metrics including boolean flags and sequence length """ if not sequence or len(sequence) < 3: return { "has_start_codon": False, "has_stop_codon": False, "length_divisible_by_3": False, "has_internal_stop_codons": False, "length": len(sequence) if sequence else 0, } # Check for canonical start codon (ATG = Methionine) has_start_codon = sequence[:3] == "ATG" # Check for stop codon at the end has_stop_codon = sequence[-3:] in ["TAA", "TAG", "TGA"] # CDS should be in-frame (length divisible by 3) length_divisible_by_3 = len(sequence) % 3 == 0 # Check for internal stop codons (premature termination) # These indicate potential annotation errors or pseudogenes has_internal_stop_codons = False if len(sequence) >= 6: # Need at least 2 codons to check for internal stops # Check all codons except the last one (which should be a stop) for i in range(0, len(sequence) - 3, 3): codon = sequence[i : i + 3] if codon in ["TAA", "TAG", "TGA"]: has_internal_stop_codons = True break return { "has_start_codon": has_start_codon, "has_stop_codon": has_stop_codon, "length_divisible_by_3": length_divisible_by_3, "has_internal_stop_codons": has_internal_stop_codons, "length": len(sequence), }

In [44]:

# =============================================================================
# Load Reference Genome (GRCh38/hg38)
# =============================================================================
# Pre-load chromosome sequences into memory for faster access during variant generation.
# Only loading standard chromosomes (1-22, X, Y) - excluding patches and alternate contigs.

import pyfaidx


fasta = {}

with pyfaidx.Fasta(REFERENCE_GENOME) as f:
    for chrom in VALID_CHROMS:
        fasta[chrom] = f[chrom][:].seq
print(f"Loaded {len(fasta)} chromosomes from {REFERENCE_GENOME}")

# ============================================================================= # Load Reference Genome (GRCh38/hg38) # ============================================================================= # Pre-load chromosome sequences into memory for faster access during variant generation. # Only loading standard chromosomes (1-22, X, Y) - excluding patches and alternate contigs. import pyfaidx fasta = {} with pyfaidx.Fasta(REFERENCE_GENOME) as f: for chrom in VALID_CHROMS: fasta[chrom] = f[chrom][:].seq print(f"Loaded {len(fasta)} chromosomes from {REFERENCE_GENOME}")

Loaded 23 chromosomes from /data/for_paper/data/reference/hg38/hg38.fa

In [45]:

# =============================================================================
# Load Annotations, Extract CDS Sequences, and Apply Quality Filters
# =============================================================================

# -----------------------------------------------------------------------------
# Step 1: Load processed GENCODE annotation
# -----------------------------------------------------------------------------
# Input: TSV file from the processing above containing
# transcript coordinates, exon boundaries, and canonical status flags
ann = pl.read_csv(ANNOTATION_FILE, separator="\t")
ann = ann.filter(pl.col("chrom").is_in(VALID_CHROMS))
print(f"Loaded {len(ann):,} transcripts from {ANNOTATION_FILE}")

# -----------------------------------------------------------------------------
# Step 2: Deduplicate transcripts with identical genomic structure
# -----------------------------------------------------------------------------
# Multiple transcript IDs can map to the same CDS coordinates (e.g., RefSeq vs Ensembl)
# Keep MANE Select > Ensembl Canonical when duplicates exist
ann = ann.sort(["is_mane_select", "is_canonical"], descending=True)
ann = ann.unique(subset=["chrom", "strand", "cdsStart", "cdsEnd", "exonStarts", "exonEnds"])
print(f"After deduplicating by genomic structure: {len(ann):,} transcripts")

# -----------------------------------------------------------------------------
# Step 3: Extract CDS sequences from reference genome
# -----------------------------------------------------------------------------
print("Extracting CDS sequences...")
sequences = [extract_cds_sequence(row, fasta) for row in ann.iter_rows(named=True)]
ann = ann.with_columns(pl.Series("cds_sequence", sequences))

# -----------------------------------------------------------------------------
# Step 4: Quality control - validate CDS sequences
# -----------------------------------------------------------------------------
print("Running quality checks...")
quality_checks = [check_cds_quality(row["cds_sequence"]) for row in ann.iter_rows(named=True)]

ann = ann.with_columns(
    [
        pl.Series("has_start_codon", [q["has_start_codon"] for q in quality_checks]),
        pl.Series("has_stop_codon", [q["has_stop_codon"] for q in quality_checks]),
        pl.Series("length_divisible_by_3", [q["length_divisible_by_3"] for q in quality_checks]),
        pl.Series("has_internal_stop_codons", [q["has_internal_stop_codons"] for q in quality_checks]),
        pl.Series("cds_length", [q["length"] for q in quality_checks]),
    ]
)

# Print quality summary
print("\n" + "=" * 60)
print("CDS Quality Summary (before filtering):")
print("=" * 60)
print(f"Total transcripts: {len(ann):,}")
print(f"Has start codon (ATG): {ann['has_start_codon'].sum():,} ({ann['has_start_codon'].mean() * 100:.1f}%)")
print(f"Has stop codon (TAA/TAG/TGA): {ann['has_stop_codon'].sum():,} ({ann['has_stop_codon'].mean() * 100:.1f}%)")
print(
    f"Length divisible by 3: {ann['length_divisible_by_3'].sum():,} ({ann['length_divisible_by_3'].mean() * 100:.1f}%)"
)
print(
    f"Has internal stop codons: {ann['has_internal_stop_codons'].sum():,} ({ann['has_internal_stop_codons'].mean() * 100:.1f}%)"
)
print(
    f"All quality criteria met: {(ann['has_start_codon'] & ann['has_stop_codon'] & ann['length_divisible_by_3'] & ~ann['has_internal_stop_codons']).sum():,}"
)

# -----------------------------------------------------------------------------
# Step 5: Apply quality filters
# -----------------------------------------------------------------------------
# Keep only transcripts that pass all quality checks
ann = ann.filter(
    pl.col("has_start_codon")
    & pl.col("has_stop_codon")
    & pl.col("length_divisible_by_3")
    & ~pl.col("has_internal_stop_codons")
)
print(f"\nAfter quality filtering: {len(ann):,} transcripts")

# -----------------------------------------------------------------------------
# Step 6: Filter to canonical transcripts and deduplicate by CDS sequence
# -----------------------------------------------------------------------------
# Keep only MANE Select or Ensembl Canonical transcripts
# Then deduplicate by CDS sequence (different transcripts can encode identical proteins)
initial_count = len(ann)
ann = ann.filter(pl.col("is_mane_select") | pl.col("is_canonical"))
print(f"After filtering to canonical transcripts: {len(ann):,} transcripts")

ann = ann.sort(["is_mane_select", "is_canonical"], descending=True)
ann = ann.unique(subset=["cds_sequence"], keep="first")
ann = ann.sort(["chrom", "txStart"])

print(f"After canonical filter + CDS deduplication: {len(ann):,} unique transcripts")
print(f"  (Removed {initial_count - len(ann):,} transcripts)")
ann.write_csv(FILTERED_TRANSCRIPTS_FILE, separator="\t")
print(f"Saved {len(ann):,} unique transcripts to {FILTERED_TRANSCRIPTS_FILE}")

# ============================================================================= # Load Annotations, Extract CDS Sequences, and Apply Quality Filters # ============================================================================= # ----------------------------------------------------------------------------- # Step 1: Load processed GENCODE annotation # ----------------------------------------------------------------------------- # Input: TSV file from the processing above containing # transcript coordinates, exon boundaries, and canonical status flags ann = pl.read_csv(ANNOTATION_FILE, separator="\t") ann = ann.filter(pl.col("chrom").is_in(VALID_CHROMS)) print(f"Loaded {len(ann):,} transcripts from {ANNOTATION_FILE}") # ----------------------------------------------------------------------------- # Step 2: Deduplicate transcripts with identical genomic structure # ----------------------------------------------------------------------------- # Multiple transcript IDs can map to the same CDS coordinates (e.g., RefSeq vs Ensembl) # Keep MANE Select > Ensembl Canonical when duplicates exist ann = ann.sort(["is_mane_select", "is_canonical"], descending=True) ann = ann.unique(subset=["chrom", "strand", "cdsStart", "cdsEnd", "exonStarts", "exonEnds"]) print(f"After deduplicating by genomic structure: {len(ann):,} transcripts") # ----------------------------------------------------------------------------- # Step 3: Extract CDS sequences from reference genome # ----------------------------------------------------------------------------- print("Extracting CDS sequences...") sequences = [extract_cds_sequence(row, fasta) for row in ann.iter_rows(named=True)] ann = ann.with_columns(pl.Series("cds_sequence", sequences)) # ----------------------------------------------------------------------------- # Step 4: Quality control - validate CDS sequences # ----------------------------------------------------------------------------- print("Running quality checks...") quality_checks = [check_cds_quality(row["cds_sequence"]) for row in ann.iter_rows(named=True)] ann = ann.with_columns( [ pl.Series("has_start_codon", [q["has_start_codon"] for q in quality_checks]), pl.Series("has_stop_codon", [q["has_stop_codon"] for q in quality_checks]), pl.Series("length_divisible_by_3", [q["length_divisible_by_3"] for q in quality_checks]), pl.Series("has_internal_stop_codons", [q["has_internal_stop_codons"] for q in quality_checks]), pl.Series("cds_length", [q["length"] for q in quality_checks]), ] ) # Print quality summary print("\n" + "=" * 60) print("CDS Quality Summary (before filtering):") print("=" * 60) print(f"Total transcripts: {len(ann):,}") print(f"Has start codon (ATG): {ann['has_start_codon'].sum():,} ({ann['has_start_codon'].mean() * 100:.1f}%)") print(f"Has stop codon (TAA/TAG/TGA): {ann['has_stop_codon'].sum():,} ({ann['has_stop_codon'].mean() * 100:.1f}%)") print( f"Length divisible by 3: {ann['length_divisible_by_3'].sum():,} ({ann['length_divisible_by_3'].mean() * 100:.1f}%)" ) print( f"Has internal stop codons: {ann['has_internal_stop_codons'].sum():,} ({ann['has_internal_stop_codons'].mean() * 100:.1f}%)" ) print( f"All quality criteria met: {(ann['has_start_codon'] & ann['has_stop_codon'] & ann['length_divisible_by_3'] & ~ann['has_internal_stop_codons']).sum():,}" ) # ----------------------------------------------------------------------------- # Step 5: Apply quality filters # ----------------------------------------------------------------------------- # Keep only transcripts that pass all quality checks ann = ann.filter( pl.col("has_start_codon") & pl.col("has_stop_codon") & pl.col("length_divisible_by_3") & ~pl.col("has_internal_stop_codons") ) print(f"\nAfter quality filtering: {len(ann):,} transcripts") # ----------------------------------------------------------------------------- # Step 6: Filter to canonical transcripts and deduplicate by CDS sequence # ----------------------------------------------------------------------------- # Keep only MANE Select or Ensembl Canonical transcripts # Then deduplicate by CDS sequence (different transcripts can encode identical proteins) initial_count = len(ann) ann = ann.filter(pl.col("is_mane_select") | pl.col("is_canonical")) print(f"After filtering to canonical transcripts: {len(ann):,} transcripts") ann = ann.sort(["is_mane_select", "is_canonical"], descending=True) ann = ann.unique(subset=["cds_sequence"], keep="first") ann = ann.sort(["chrom", "txStart"]) print(f"After canonical filter + CDS deduplication: {len(ann):,} unique transcripts") print(f" (Removed {initial_count - len(ann):,} transcripts)") ann.write_csv(FILTERED_TRANSCRIPTS_FILE, separator="\t") print(f"Saved {len(ann):,} unique transcripts to {FILTERED_TRANSCRIPTS_FILE}")

Loaded 64,488 transcripts from /data/for_paper/data/reference/gencode.v47.basic.annotation.processed.tsv
After deduplicating by genomic structure: 51,650 transcripts
Extracting CDS sequences...
Running quality checks...

============================================================
CDS Quality Summary (before filtering):
============================================================
Total transcripts: 51,650
Has start codon (ATG): 51,419 (99.6%)
Has stop codon (TAA/TAG/TGA): 51,341 (99.4%)
Length divisible by 3: 51,573 (99.9%)
Has internal stop codons: 113 (0.2%)
All quality criteria met: 51,061

After quality filtering: 51,061 transcripts
After filtering to canonical transcripts: 19,407 transcripts
After canonical filter + CDS deduplication: 19,310 unique transcripts
  (Removed 31,751 transcripts)
Saved 19,310 unique transcripts to /data/for_paper/data/reference/gencode.v47.basic.annotation.processed.filtered.tsv