Python Fundamentals for Bioinformatics

mkdir -p ~/projects/fasta-stats/data/raw ~/projects/fasta-stats/scripts ~/projects/fasta-stats/outputs
cd ~/projects/fasta-stats
cat > data/raw/example.fasta <<'EOF'
>seq1
ACGTACGTNNACGT
>seq2
acgtacgtacgtacgtacgt
>seq3
AGCTTAGC--NNNN
EOF
Rscript --version

#!/usr/bin/env Rscript
# Read FASTA records, validate A/C/G/T/N after normalization, and write CSV.
# Usage: Rscript compute_fasta_stats.R INPUT.fasta OUTPUT.csv
args <- commandArgs(trailingOnly = TRUE)
if (length(args) != 2L) {
  stop("Usage: Rscript compute_fasta_stats.R INPUT.fasta OUTPUT.csv", call. = FALSE)
}
input_path <- normalizePath(args[[1]], mustWork = TRUE)
output_path <- normalizePath(args[[2]], mustWork = FALSE)
if (identical(input_path, output_path)) {
  stop("Input and output paths must be different", call. = FALSE)
}

fasta_records <- function(path) {
  lines <- readLines(path, warn = FALSE, encoding = "UTF-8")
  records <- list()
  current_id <- NULL
  sequence_parts <- character()

  save_record <- function() {
    if (!is.null(current_id)) {
      records[[length(records) + 1L]] <<- list(
        id = current_id,
        sequence = paste(sequence_parts, collapse = "")
      )
    }
  }

  for (raw_line in lines) {
    line <- trimws(raw_line)
    if (!nzchar(line)) next
    if (startsWith(line, ">")) {
      save_record()
      current_id <- trimws(substring(line, 2L))
      if (!nzchar(current_id)) stop("FASTA header is empty", call. = FALSE)
      sequence_parts <- character()
    } else {
      if (is.null(current_id)) {
        stop("Sequence text found before the first FASTA header", call. = FALSE)
      }
      sequence_parts <- c(sequence_parts, line)
    }
  }
  save_record()
  records
}

normalize_sequence <- function(sequence) {
  toupper(gsub("[[:space:]]+", "", sequence, perl = TRUE))
}
valid_bases <- c("A", "C", "G", "T", "N")
records <- fasta_records(input_path)
rows <- lapply(records, function(record) {
  sequence <- normalize_sequence(record$sequence)
  valid <- nzchar(sequence) && all(strsplit(sequence, "", fixed = TRUE)[[1L]] %in% valid_bases)
  gc_text <- ""
  if (valid) {
    bases <- strsplit(sequence, "", fixed = TRUE)[[1L]]
    gc_text <- sprintf("%.5f", sum(bases %in% c("G", "C")) / length(bases))
  }
  data.frame(
    id = record$id,
    length = nchar(sequence, type = "chars"),
    gc_fraction = gc_text,
    valid = if (valid) "yes" else "no",
    stringsAsFactors = FALSE,
    check.names = FALSE
  )
})
result <- if (length(rows)) {
  do.call(rbind, rows)
} else {
  data.frame(id = character(), length = integer(), gc_fraction = character(),
             valid = character(), stringsAsFactors = FALSE)
}
output_dir <- dirname(output_path)
if (!dir.exists(output_dir) && !dir.create(output_dir, recursive = TRUE)) {
  stop("Could not create output directory: ", output_dir, call. = FALSE)
}
write.csv(result, output_path, row.names = FALSE, na = "", fileEncoding = "UTF-8")

Rscript scripts/compute_fasta_stats.R data/raw/example.fasta outputs/sequence_stats.csv
cat outputs/sequence_stats.csv

What you will build

The lab ends with a CSV row for each FASTA record, containing its identifier, normalized length, GC fraction when the sequence passes the chosen alphabet check, and a validity flag. You will also be able to explain the denominator used for that fraction.

• Use Python variables, strings, lists, dictionaries, loops, conditions, and functions for small bioinformatics tasks.
• Read and validate plain-text biological sequence files (FASTA) safely and compute basic statistics: sequence length and GC content.
• Create a reproducible, modular script with a clear input and output path, and run it inside an isolated virtual environment.
• Inspect results and understand what the computed statistics imply and what they do not.

Prerequisites:

• Basic familiarity with the command line (cd, ls) and an editor (vim, nano, VS Code).
• Python 3.10 or newer installed. This lesson uses explicit commands for Python 3; adjust to your system's python executable if needed.
• Approximately 30–60 minutes of hands-on time.

Why validate the input before summarizing it?

A FASTA file is plain text organized as records: a header line begins with `>`, followed by one or more sequence lines. A small parser lets you inspect the mechanics, but its assumptions must be visible. Here the script accepts only A, C, G, T, and N after uppercasing; it flags other symbols rather than silently removing them. That is a teaching rule, not a universal definition of valid biology.

Lab: turn FASTA records into a checked summary

Work in a project directory so the input, script, and output have distinct paths. The inline FASTA fixture is synthetic teaching data, not a biological result. Commands below target a POSIX shell.

Files you should have before starting (create these or replace with your own):

• data/raw/example.fasta: the synthetic fixture created in step 1, or a FASTA file you have permission to use.
• scripts/compute_fasta_stats.py: the script created in step 3.

1. Set up the project and inspect the fixture

# Make a project folder and create input directories
mkdir -p ~/projects/python-bioinfo/data/raw ~/projects/python-bioinfo/scripts
cd ~/projects/python-bioinfo

# Create a small example FASTA if you do not already have one
cat > data/raw/example.fasta <<'EOF'
>seq1
ACGTACGTNNACGT
>seq2
acgtacgtacgtacgtacgt
>seq3
AGCTTAGC--NNNN
EOF

# Confirm the file exists and show contents
ls -l data/raw/example.fasta
cat data/raw/example.fasta

What this produces: a FASTA file at data/raw/example.fasta containing three sequence records. The fixture has three synthetic records. The third contains `-`, so it should be marked invalid by this script. This makes the validation rule observable without presenting fabricated biological data.

2. Create and activate a virtual environment

# Create an isolated environment (will create ~/projects/python-bioinfo/.venv)
python3 -m venv .venv

# Activate the environment (POSIX shells)
source .venv/bin/activate

# Upgrade pip and install recommended packages (optional)
python -m pip install biopython

# Record exact versions for reproducibility
python -m pip freeze > envs-requirements.txt

What this produces: a virtual environment in .venv, and envs-requirements.txt listing installed package versions. The optional Biopython comparison in step 5 requires an internet connection to install; skip it if you only want to run the standard-library script. To recreate the recorded package environment later, use `python -m pip install -r envs-requirements.txt`.

3. Write a small script with explicit checks

Create scripts/compute_fasta_stats.py. The script uses only Python’s standard library. It reads records one at a time, rejects sequence text before the first header and empty headers, and refuses to overwrite the input file accidentally. Its alphabet rule is explicit; adjust it only when your downstream task supports the added symbols.

#!/usr/bin/env python3
"""
scripts/compute_fasta_stats.py

Read a FASTA file, validate sequences, compute length and GC fraction for each record,
and write a CSV file with summary statistics.

Usage:
  python scripts/compute_fasta_stats.py data/raw/example.fasta outputs/sequence_stats.csv
"""
import sys
import csv
from pathlib import Path

VALID_BASES = set("ACGTN")  # Accept uppercase only after normalization

def normalize(seq):
    "Strip whitespace and convert to uppercase."
    return "".join(seq.split()).upper()

def is_valid_dna(seq):
    "Return True if sequence contains only expected bases A,C,G,T,N."
    return bool(seq) and all(ch in VALID_BASES for ch in seq)

def gc_fraction(seq):
    "Compute GC fraction safely; assumes seq is normalized uppercase."
    if len(seq) == 0:
        return 0.0
    gc = sum(1 for ch in seq if ch in ("G", "C"))
    return gc / len(seq)

def fasta_records(path):
    """Simple record-by-record FASTA parser.
    Yields tuples: (header, sequence_string).
    This parser treats any line starting with '>' as header.
    """
    with open(path, "r", encoding="utf-8") as fh:
        name = None
        parts = []
        for raw in fh:
            line = raw.strip()
            if not line:
                continue
            if line.startswith(">"):
                if name is not None:
                    yield name, "".join(parts)
                name = line[1:].strip()
                if not name:
                    raise ValueError("FASTA header is empty")
                parts = []
            elif name is None:
                raise ValueError("Sequence text found before the first FASTA header")
            else:
                parts.append(line)
        if name is not None:
            yield name, "".join(parts)

def main():
    if len(sys.argv) != 3:
        print("Usage: python scripts/compute_fasta_stats.py INPUT.fasta OUTPUT.csv", file=sys.stderr)
        sys.exit(2)

    in_path = Path(sys.argv[1])
    out_path = Path(sys.argv[2])
    if in_path.resolve() == out_path.resolve():
        print("Input and output paths must be different", file=sys.stderr)
        sys.exit(2)
    out_path.parent.mkdir(parents=True, exist_ok=True)

    with open(out_path, "w", newline="", encoding="utf-8") as out_fh:
        writer = csv.DictWriter(out_fh, fieldnames=["id", "length", "gc_fraction", "valid"])
        writer.writeheader()
        for header, seq in fasta_records(in_path):
            seq_norm = normalize(seq)
            valid = is_valid_dna(seq_norm)
            length = len(seq_norm)
            gc = gc_fraction(seq_norm) if valid else ""
            writer.writerow({
                "id": header,
                "length": length,
                "gc_fraction": f"{gc:.5f}" if gc != "" else "",
                "valid": "yes" if valid else "no"
            })

if __name__ == "__main__":
    main()

What this produces: a script that reads a FASTA file, validates sequences against a simple DNA alphabet (A,C,G,T,N), computes length and GC fraction for each record, and writes data to the output CSV path you provide.

4. Run the script

# Run the script with the example FASTA and write results to outputs/
python scripts/compute_fasta_stats.py data/raw/example.fasta outputs/sequence_stats.csv

# Show the resulting CSV
cat outputs/sequence_stats.csv

Expected file: outputs/sequence_stats.csv with columns: id, length, gc_fraction, valid. The exact values depend on the input. For the fixture, the first record has length 14 and GC fraction 0.42857; the second has length 20 and fraction 0.50000; the third is invalid because it contains `-`, so its GC field is blank. These are arithmetic checks on synthetic strings, not empirical findings. You can load the CSV into pandas if that package is installed. Treat files from untrusted sources cautiously when opening them in spreadsheet software.

5. Optional: parse with Biopython for comparison

If you installed Biopython in the virtual environment, use SeqIO to parse records, then apply the same normalization and alphabet rule as the main script. `SeqIO.parse` yields records from an iterator and supports multiple sequence formats; parsing does not replace validating that symbols and metadata fit your analysis.

# Run an interactive comparison (while env is activated)
python - <<'PY'
from Bio import SeqIO
from pathlib import Path

valid_bases = set("ACGTN")
for rec in SeqIO.parse(str(Path("data/raw/example.fasta")), "fasta"):
    seq = "".join(str(rec.seq).split()).upper()
    valid = bool(seq) and all(base in valid_bases for base in seq)
    gc = sum(1 for base in seq if base in "GC") / len(seq) if valid else ""
    gc_text = f"{gc:.5f}" if valid else ""
    print(rec.id, len(seq), gc_text, "yes" if valid else "no")
PY

What this produces: console output listing record id, normalized length, GC fraction (blank for invalid records), and validity. For the supplied fixture, seq3 is reported invalid with no GC value, matching the standard-library script. Biopython provides a maintained parser; the validity rule remains the workflow's explicit responsibility.

Output inspection and interpretation

When you open outputs/sequence_stats.csv, check the following:

• All expected record IDs are present. Missing IDs may indicate a parsing error or that the input FASTA lines were malformed.
• Length values reflect the number of characters after normalization (spaces and tabs removed). Very short sequences (length 0–3) should be inspected manually because they may be fragments, adapter sequences, or corrupted lines.
• gc_fraction is the proportion of G or C bases (0.0–1.0). This script divides the count of G and C characters by the full normalized sequence length, including `N` in the denominator. Compare values only when the denominator rule and sequence scope match; GC content alone is not a quality verdict.
• The valid column flags records where non-standard characters were found (for this lesson that means characters outside A,C,G,T,N). If valid is "no", the gc_fraction field is left blank and you should inspect the corresponding sequence in the FASTA file.

What a plausible result means and what it does NOT prove:

• If a sequence has gc_fraction = 0.42, this simply means 42% of the characters in the provided string are G or C. It does not prove the sequence is error-free or biologically meaningful; contamination, assembly artifacts, or sequencing biases can change GC content.
• Short sequences or sequences labeled invalid may be genuine biological data (e.g., primers with ambiguous bases) or may indicate file corruption. Manual inspection or comparison to original data sources is required.
• Aggregate summaries (mean GC across records) are descriptive. They are not a substitute for statistical testing if you want to compare groups or draw biological conclusions.

Troubleshooting: symptoms, likely causes, and fixes

Symptom: The script crashes with UnicodeDecodeError when reading the FASTA.

• Cause: The input file uses a non-UTF-8 encoding (e.g., Windows-1252) or contains binary content.
• Fix: Convert the file to UTF-8 or open with the correct encoding. Example: iconv -f WINDOWS-1252 -t UTF-8 input.fasta -o output.fasta. Always inspect binary files (e.g., using file and head commands) before processing.

Symptom: Many sequences are flagged valid=no.

• Cause: The input FASTA contains characters outside the accepted set (e.g., gap character '-', IUPAC codes like R/Y, or quality scores accidentally included).
• Fix: Decide what extended characters you want to accept (IUPAC codes: R, Y, S, W, K, M, B, D, H, V) and extend VALID_BASES accordingly. Alternatively, pre-filter or clean sequences to remove gap characters if appropriate for your analysis.

Symptom: Output CSV is missing records or header only.

• Cause: The file may be empty, contain no header beginning with `>`, or use malformed record structure. The parser now rejects sequence text before the first header, but an empty file still produces a header-only CSV.
• Fix: Confirm the FASTA file contains at least one header line (beginning with '>') and sequence lines. Use wc -l and grep '^>' to count headers: grep -c '^>' data/raw/example.fasta.

Symptom: Incorrect GC values compared to external tools.

• Cause: Differences in normalization: some tools ignore ambiguous bases or gaps when computing fractions, others count them in the denominator.
• Fix: Make your normalization explicit. If you want GC fraction computed as GC / (A+C+G+T) excluding ambiguous bases, implement that rule rather than using the raw length.

Limitations, version and reference-data caveats, and responsible use

Limitations of this lesson and script:

• The simple parser is educational and reads records incrementally, though it holds each record sequence in memory. It does not validate FASTA headers against a controlled schema (for example, embedded tabs or spaces), support FASTQ, or provide indexed access. For production work, use a maintained parser such as Biopython SeqIO and validate records against your workflow’s requirements.
• Validation here is limited to a user-defined alphabet (A, C, G, T, N). Biological sequences sometimes include IUPAC ambiguity codes or gap characters (e.g., '-'). Decide which characters your pipeline should accept and document that choice.
• Floating-point representation of GC fraction uses simple rounding for readability; for downstream statistical analysis you may want more precision or native numeric types (not formatted strings).

Version and reference-data caveats:

• Different versions of Biopython, pandas, or Python may change behavior (for example, encoding defaults or CSV dialect handling). Record exact versions with pip freeze or conda list for reproducibility.
• Reference genome sequences and annotations can change over time. If you are comparing GC content to published values, ensure you are using the same assembly and coordinate convention as the reference.

Responsible use notes:

• Do not assume computed statistics imply biological meaning without appropriate controls and metadata assessment (sample origin, library preparation, sequencing platform).
• Be cautious when sharing sequence-derived files that may include sensitive or personally identifying information. Follow institutional policies and data-use agreements.
• Never run or paste unreviewed code into tools that will give it network access or elevate privileges. The examples here are safe for local file processing and avoid destructive commands (no rm, mv, or system calls that touch other directories).

Next-step checklist

1. Run the script on your own FASTA files and confirm outputs/sequence_stats.csv contains expected rows and values.
2. Decide how to treat ambiguous bases and update VALID_BASES and gc_fraction semantics accordingly.
3. For large inputs, consider a streaming parser such as Biopython SeqIO; this script still holds each individual record sequence in memory.
4. Write small unit tests for gc_fraction and is_valid_dna using the pytest framework to ensure behavior is stable as you refactor.
5. Keep the requirements file with the project and document how to recreate the environment. A freeze records package versions, but does not by itself capture Python, operating-system, or system-library details.

References

• The Python 3 documentation: language reference and library documentation.
• Biopython SeqIO API documentation: parser and record-iterator reference.
• Python venv documentation: environment creation and activation.
• Pandas documentation: for tabular data manipulation, reading/writing CSV/TSV, and summary statistics.
• NCBI FASTA format guidance: nucleotide submission header and sequence conventions.