Case Study: Parsing Real-World Text Data — Dr. Patel's FASTA File Processor

The Scenario

Dr. Anika Patel is a molecular biologist at the University of Michigan. Her research focuses on identifying genetic markers associated with antibiotic resistance in bacteria. Every week, she receives sequencing results from the campus genomics core facility — sometimes hundreds of DNA sequences at a time.

The sequences arrive in FASTA format, one of the most widely used text formats in bioinformatics. A FASTA file is plain text with a simple structure:

>sequence_id description text
ATCGATCGATCGATCG
GCTAGCTAGCTAGCTA
>next_sequence_id more description
AAAGGGCCCTTTTAAA

Each sequence starts with a header line beginning with >, followed by one or more lines of nucleotide characters (A, T, C, G). That's it. No XML, no JSON, no binary encoding — just text. And that means the entire processing pipeline is fundamentally string processing.

Dr. Patel's problem: she needs to extract specific information from these files — sequence names, lengths, nucleotide compositions, and GC content — and she needs to do it reliably for files containing anywhere from 10 to 10,000 sequences.

The Manual Approach (Before Python)

Before learning Python, Dr. Patel used a combination of Excel and a web-based sequence-statistics tool to process her data. The workflow looked like this:

1. Open the FASTA file in a text editor
2. Manually count sequences (count the > lines)
3. Copy individual sequences into a web form for analysis
4. Record results in a spreadsheet
5. Repeat for each sequence

For a file with 50 sequences, this took about three hours. For 500 sequences, it was simply not feasible by hand.

Building the Parser: Step by Step

Step 1: Reading the File Structure

The first challenge is understanding what each line means. In FASTA, lines serve different roles:

def classify_lines(fasta_text):
    """Classify each line in FASTA text as header, sequence, or empty."""
    for line_num, line in enumerate(fasta_text.splitlines(), start=1):
        stripped = line.strip()
        if not stripped:
            line_type = "EMPTY"
        elif stripped.startswith(">"):
            line_type = "HEADER"
        else:
            line_type = "SEQUENCE"
        print(f"  Line {line_num:3d}: [{line_type:8s}] {stripped[:50]}")


sample = """>gi|12345|ref|NM_001.2| BRCA1 gene, Homo sapiens
ATCGATCGATCGATCGATCG
GCTAGCTAGCTAGCTAGCTA

>gi|67890|ref|NM_002.1| TP53 tumor suppressor
TTTTAAAACCCCGGGGAAAA
CCCCTTTTGGGGAAAA"""

classify_lines(sample)

Output:
  Line   1: [HEADER  ] >gi|12345|ref|NM_001.2| BRCA1 gene, Homo sapiens
  Line   2: [SEQUENCE] ATCGATCGATCGATCGATCG
  Line   3: [SEQUENCE] GCTAGCTAGCTAGCTAGCTA
  Line   4: [EMPTY   ]
  Line   5: [HEADER  ] >gi|67890|ref|NM_002.1| TP53 tumor suppressor
  Line   6: [SEQUENCE] TTTTAAAACCCCGGGGAAAA
  Line   7: [SEQUENCE] CCCCTTTTGGGGAAAA

The string methods at work here: strip() to handle whitespace, startswith(">") to identify headers, and splitlines() to iterate line by line.

Step 2: Parsing Headers

FASTA headers encode structured information, but the format varies between databases. Older NCBI files, which are still common in lab archives, use pipe-separated fields:

>gi|12345|ref|NM_001.2| BRCA1 gene, Homo sapiens

Reading left to right: the GI number (12345), the source database (ref, for RefSeq), the accession number with its version (NM_001.2), and then a free-text description. (NCBI stopped putting GI numbers in its FASTA headers in 2016; newer files start with just the accession, like >NM_001.2 BRCA1 gene. Real-world parsers have to cope with both. The function below only understands the older pipe format, but it checks how many fields it got, so a newer header won't crash it.) To build and test her parser, Dr. Patel uses a few short records for well-known human genes such as BRCA1 and TP53 before pointing it at her bacterial data. The ID numbers in these examples are shortened placeholders, so don't go looking them up.

Dr. Patel needs to extract the accession number and the description:

def parse_ncbi_header(header_line):
    """Parse an NCBI-format FASTA header into components."""
    # Remove the leading '>'
    header = header_line.lstrip(">").strip()

    # Split on '|' to get fields
    parts = header.split("|")

    result = {
        "raw": header,
        "gi_number": None,
        "database": None,
        "accession": None,
        "description": None,
    }

    # parts is now ['gi', '12345', 'ref', 'NM_001.2', ' BRCA1 gene, Homo sapiens']
    if len(parts) >= 4:
        result["gi_number"] = parts[1].strip()
        result["database"] = parts[2].strip()
        result["accession"] = parts[3].strip()
    if len(parts) >= 5:
        # Everything after the fourth '|' is the description
        result["description"] = parts[4].strip()

    return result


# Test
header = ">gi|12345|ref|NM_001.2| BRCA1 gene, Homo sapiens"
info = parse_ncbi_header(header)
for key, value in info.items():
    print(f"  {key:>15s}: {value}")

Output:
              raw: gi|12345|ref|NM_001.2| BRCA1 gene, Homo sapiens
        gi_number: 12345
         database: ref
        accession: NM_001.2
      description: BRCA1 gene, Homo sapiens

Key string techniques: lstrip(">") to remove the header marker, split("|") to break the pipe-delimited fields into a list, indexing to pick out each field by position, and strip() to clean the stray spaces around them. Notice the len(parts) checks: if a header has fewer fields than expected, the function leaves those entries as None instead of crashing with an IndexError.

Step 3: Computing Sequence Statistics

Now Dr. Patel needs to calculate metrics for each sequence — length, nucleotide counts, and GC content:

def sequence_stats(sequence):
    """Calculate statistics for a DNA sequence string."""
    # Clean the sequence: remove whitespace, convert to uppercase
    clean = sequence.replace("\n", "").replace(" ", "").upper()

    length = len(clean)
    if length == 0:
        return None

    # Count each nucleotide
    counts = {}
    for nucleotide in "ATCG":
        counts[nucleotide] = clean.count(nucleotide)

    # GC content: percentage of G + C
    gc = (counts["G"] + counts["C"]) / length * 100

    # Check for unexpected characters
    known = sum(counts.values())
    unknown = length - known

    return {
        "length": length,
        "counts": counts,
        "gc_content": gc,
        "unknown_chars": unknown,
    }


# Test
seq = "ATCGATCGATCGATCG\nGCTAGCTAGCTAGCTA"
stats = sequence_stats(seq)
print(f"  Length:      {stats['length']}")
print(f"  Nucleotides: {stats['counts']}")
print(f"  GC content:  {stats['gc_content']:.1f}%")
print(f"  Unknown:     {stats['unknown_chars']}")

Output:
  Length:      32
  Nucleotides: {'A': 8, 'T': 8, 'C': 8, 'G': 8}
  GC content:  50.0%
  Unknown:     0

Step 4: The Complete Parser

Putting it all together — a function that reads a FASTA text and returns structured data for every sequence:

def parse_fasta(fasta_text):
    """Parse FASTA-formatted text into a list of sequence records.

    Each record is a dict with 'header', 'sequence', and 'stats' keys.
    """
    records = []
    current_header = None
    current_seq_lines = []

    for line in fasta_text.splitlines():
        stripped = line.strip()
        if not stripped:
            continue  # skip blank lines
        if stripped.startswith(">"):
            # Save the previous record (if any)
            if current_header is not None:
                full_seq = "".join(current_seq_lines)
                records.append({
                    "header": parse_ncbi_header(current_header),
                    "sequence": full_seq,
                    "stats": sequence_stats(full_seq),
                })
            # Start a new record
            current_header = stripped
            current_seq_lines = []
        else:
            current_seq_lines.append(stripped)

    # Don't forget the last record
    if current_header is not None:
        full_seq = "".join(current_seq_lines)
        records.append({
            "header": parse_ncbi_header(current_header),
            "sequence": full_seq,
            "stats": sequence_stats(full_seq),
        })

    return records


# Test with sample data
fasta_data = """>gi|12345|ref|NM_001.2| BRCA1 gene
ATCGATCGATCGATCGATCG
GCTAGCTAGCTAGCTAGCTA

>gi|67890|ref|NM_002.1| TP53 gene
TTTTAAAACCCCGGGG

>gi|11111|ref|NM_003.3| MYC oncogene
GGGGCCCCAAAATTTT
GGGGCCCCAAAATTTT"""

records = parse_fasta(fasta_data)

# Display results in a formatted table
print(f"\n  {'Gene':<20s} {'Accession':<12s} {'Length':>7s} {'GC%':>6s}")
print(f"  {'-'*20} {'-'*12} {'-'*7} {'-'*6}")

for rec in records:
    desc = rec["header"]["description"] or "Unknown"
    acc = rec["header"]["accession"] or "N/A"
    length = rec["stats"]["length"]
    gc = rec["stats"]["gc_content"]
    print(f"  {desc:<20s} {acc:<12s} {length:>7d} {gc:>5.1f}%")

print(f"\n  Total sequences: {len(records)}")

Output:
  Gene                 Accession     Length    GC%
  -------------------- ------------ ------- ------
  BRCA1 gene           NM_001.2          40  50.0%
  TP53 gene            NM_002.1          16  50.0%
  MYC oncogene         NM_003.3          32  50.0%

  Total sequences: 3

String Techniques Used

This case study demonstrates nearly every string technique from Chapter 7:

Discussion Questions

1. Robustness: The parser assumes well-formed FASTA data. What could go wrong with real-world files? Consider: malformed headers, non-standard nucleotide characters (N, R, Y), very long lines, mixed line endings (\n vs \r\n). How would you handle each?

2. Scalability: This parser loads the entire file into memory as a string. For a file with 10 million sequences, this could consume gigabytes of RAM. How might you modify the approach to process one sequence at a time? (Hint: think about reading files line by line in Chapter 10.)

3. Immutability matters: The sequence_stats function calls sequence.replace().replace().upper(). Each call creates a new string. For a 3-billion-character human genome sequence, how many copies would this chain create? Is there a more efficient approach?

4. Abstraction (Ch 6 callback): We built this parser as several small functions, each doing one job. How does this decomposition relate to the function design principles from Chapter 6? Could you test each function independently?

5. Real-world impact: Dr. Patel's manual process took 3 hours for 50 sequences. Estimate how long the Python parser would take for 50 sequences, 500 sequences, and 5,000 sequences. What does this tell you about the value of automation?