pythonforbio.
[WASM idle]
pandas11/21

Reading and Writing Data

Starting sandbox…
Beginnerlesson

Reading and Writing Data

In one line: read_csv is the workhorse and it has ~50 parameters — the five that matter are sep, dtype, comment, na_values and usecols.

read_csv, properly

df = pd.read_csv(
    "variants.tsv",
    sep="\t",
    comment="#",                       # skip VCF-style header lines
    dtype={"chrom": "category", "pos": "int32", "sample_id": "str"},
    na_values=[".", "NA", "-", ""],    # what counts as missing in THIS file
    usecols=["chrom", "pos", "ref", "alt", "qual"],
    nrows=1000,                        # peek at a huge file first
    low_memory=False,
)

Always pass dtype= for genomic files. Type inference gets these wrong constantly:

  • A chromosome column containing 122, X, Y, MT becomes str only because of the letters. Read a chr1-only subset and it silently becomes int64, and your merges stop matching.
  • Sample IDs like 0012 lose their leading zeros if inferred as integers.
  • A column with one "." in it becomes object instead of numeric.

comment="#" matters because VCF, GFF and SAM all use # header lines.

Other formats

pd.read_csv(...)                      # also handles .gz, .bz2, .zip, .xz automatically
pd.read_parquet("data.parquet")       # fast, typed, compressed — the best intermediate
pd.read_excel("samples.xlsx", sheet_name="metadata")
pd.read_json("results.json", lines=True)
pd.read_sql("SELECT * FROM variants", conn)
pd.read_hdf("store.h5", key="expr")
pd.read_fwf("fixed_width.txt", widths=[10, 8, 6])
pd.read_clipboard()                   # genuinely useful for one-off pastes

read_csv detects compression from the extension, so pd.read_csv("counts.tsv.gz", sep="\t") just works — no gzip.open needed.

Parquet is the right intermediate format

df.to_parquet("intermediate.parquet")       # preserves dtypes exactly
df = pd.read_parquet("intermediate.parquet")
CSV/TSV Parquet
Dtypes preserved no — re-inferred every read yes
Size baseline ~5–10× smaller
Read speed baseline ~10× faster
Column subset read must read all reads only what you ask
Human readable yes no
Universal tool support yes good and improving

Use TSV at the boundaries (input from collaborators, output for the supplement) and Parquet for everything in between. Re-parsing a 2 GB TSV on every pipeline run is a self-inflicted wound.

Large files

# chunked processing
results = []
for chunk in pd.read_csv("huge.tsv", sep="\t", chunksize=100_000):
    results.append(chunk.query("qual > 30").groupby("chrom").size())
total = pd.concat(results).groupby(level=0).sum()

# only the columns you need
pd.read_csv("huge.tsv", sep="\t", usecols=["chrom", "pos", "af"])

# check the shape before committing
pd.read_csv("huge.tsv", sep="\t", nrows=5)

If chunking is not enough, the answer is polars or DuckDB, both of which query files larger than RAM directly:

import duckdb
duckdb.sql("SELECT chrom, count(*) FROM 'variants.tsv' GROUP BY chrom").df()

Writing

df.to_csv("out.tsv", sep="\t", index=False)     # index=False almost always
df.to_csv("out.tsv.gz", sep="\t", index=False)  # compression from extension
df.to_parquet("out.parquet")
df.to_excel("out.xlsx", sheet_name="results", index=False)
df.to_markdown()                                 # for a README or a notebook
df.to_latex()                                    # for a manuscript table

index=False unless the index is real data. Otherwise you get a spurious unnamed column that reappears as Unnamed: 0 on the next read, and it accumulates every round-trip.

Round-tripping properly:

df.to_csv("out.tsv", sep="\t")                  # writes the index
df = pd.read_csv("out.tsv", sep="\t", index_col=0)   # reads it back

Bioinformatics parsing patterns

VCF — skip the meta lines, find the header:

def read_vcf(path):
    import gzip
    opener = gzip.open if str(path).endswith(".gz") else open
    with opener(path, "rt") as fh:
        for line in fh:
            if line.startswith("#CHROM"):
                cols = line.lstrip("#").rstrip().split("\t")
                break
    return pd.read_csv(path, sep="\t", comment="#", names=cols,
                       dtype={"CHROM": "category", "POS": "int32"})

See Applied - Reading VCF with pandas for the INFO-field work.

GFF/GTF — nine fixed columns, no header:

GFF_COLS = ["seqid", "source", "type", "start", "end",
            "score", "strand", "phase", "attributes"]
gff = pd.read_csv(path, sep="\t", comment="#", names=GFF_COLS,
                  dtype={"seqid": "category", "type": "category",
                         "start": "int32", "end": "int32"})

BED — 0-based, variable column count:

BED_COLS = ["chrom", "start", "end", "name", "score", "strand"]
bed = pd.read_csv(path, sep="\t", header=None,
                  names=BED_COLS[:n_cols], dtype={"chrom": "category"})

FASTA/FASTQ — not tabular. Use SeqIO.

Validation on load

Worth building into any pipeline:

def load_variants(path):
    df = pd.read_csv(path, sep="\t", dtype=SCHEMA, na_values=["."])

    assert not df.empty, f"{path} is empty"
    assert df["pos"].gt(0).all(), "non-positive positions"
    assert df["chrom"].isin(VALID_CHROMS).all(), \
        f"unexpected: {set(df['chrom']) - VALID_CHROMS}"
    if df.duplicated(["chrom", "pos", "ref", "alt"]).any():
        warnings.warn("duplicate variants present")
    return df

Five lines that catch the wrong file, a coordinate-system mix-up, and a chromosome-naming mismatch (chr1 vs 1) before they reach your results. That naming mismatch in particular — Ensembl uses 1, UCSC uses chr1 — silently produces empty joins.

Common mistakes

  • No dtype= on genomic files.
  • No comment="#" on VCF/GFF, so header lines become data rows.
  • Wrong sep. The default is ,; genomics is \t. Symptom: one giant column.
  • index=False forgotten, giving Unnamed: 0 columns that multiply.
  • Not handling the file's own missing-value convention. . in VCF, NA in R output, - in some annotation files.
  • Re-parsing large TSVs instead of caching to Parquet.
  • chr1 vs 1 across files → empty merges with no error.
  • Reading a whole 10 GB file to compute one summary. Chunk, or use DuckDB.

See also

DataFrame · pandas dtypes · Missing Data in pandas · Applied - Reading VCF with pandas · Applied - Genomic File Formats · NumPy File IO

lesson example

No output yet — run the code to populate this drawer.