Reading and Writing Data

Core Principle: Data ingestion in pandas relies on pd.read_csv() parameterized with explicit sep, dtype, comment, and na_values arguments for delimited text, alongside pd.read_parquet() for high-throughput typed binary storage.

read_csv, properly

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

• A chromosome column containing 1–22, 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:

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

Exercise

Write vcf_summary(chrom, min_qual), which reads the file variants.vcf from the working directory. Stdin gives two lines: a chromosome name and a minimum quality. Return a single string.

The file is a VCF, so it has no usable header row: every line before the data starts with #, and the file's own missing-value marker is .. Read it with pandas.read_csv, supplying the separator, comment="#", the eight column names CHROM POS ID REF ALT QUAL FILTER INFO, and the file's missing marker so that an absent ID arrives as a missing value rather than the string ".". Keep the records on the requested chromosome whose QUAL is greater than or equal to min_qual, and return <chrom> n=<count> mean_qual=<mean to 3 dp> named=<count of records with an ID> pass=<count of records whose FILTER is PASS>. If no record survives: which includes a chromosome that does not appear in the file at all: return <chrom> n=0 and nothing more.