Repository: tskit-dev/msprime (language: python)

Issue:
Text file format for pedigree data input.
As discussed in #1913, the standard text files like plink's FAM don't contain enough information and we need some custom file format to specify our pedigrees.

@LukeAndersonTrocme can you give some examples and maybe sketch out what you'd like to have as the text file format input?


Discussion / hints:
I may have to follow up on the first part once I do a little bit of poking around to see what I can find regarding examples of standard fam formats with extra columns attached. The genealogies that I have worked with or seen in other papers often resemble fam formats with the first three columns being individual and parental IDs, some have sex and some space-time metadata. 

As far as what I'd like to see in a text format, I would imagine that something containing the bare minimum information needed for the simulations would be a good starting place. (on that note, please let me know if I'm missing anything here). Here's what I have in mind:

```
IID, MID, FID, TIME, IS_SAMPLE
1,     3,   4,  0,  1
2,     3,   4,  0,  1
3,    NA,  NA,  1, -1
4,    NA,  NA,  1, -1
```

As mentioned in #1913, a **population** column could be useful for simulating beyond the pedigree, but I think this would be best implemented as an _advanced_ feature that can only be done using the API. This has the added benefit of avoiding the need to verify the validity of the population IDs. (unless there are good reasons to include population IDs in this basic format?)

One thought I have is whether we embrace the three column ped format for this vanilla simulator and offer some pre-processing for these users. I suspect that some users may not know how to get the time column. I have a short R function that recursively climbs a tree to find the deepest node depth for each individual and takes as input a three column ped file and a list of probands. (Runs on 1.2 million probands of the French-Canadian genealogy in 15 seconds on my laptop, happy to paste the code over if it's useful). 

Just one note about the list of probands, the only reason why I specify probands in my case is because we have "dead-ends" in the French-Canadian genealogy and wanted a minimum date to be considered a _present day_ proband, otherwise we could just automatically consider probands to be all individuals who are not fathers or mothers. 

The real advantage with this is that we would start from the basic three column ped file that many people would be used to, and we would add `time` and `is_sample` by running this pre-processing step. Of course, users may also specify probands/samples themselves. 

edit: Am I right in assuming that all we need apart from the three ped file columns are `time` and `is_sample` ?
Just following up on the last comments from https://github.com/tskit-dev/msprime/pull/1913#issuecomment-968347694

Completely agree with the issues about lacking column names, and not being able to use a 3-column FAM file. However, as I outlined above, it seems feasible to include a little FAM file parser to include the additional columns needed (and likely useful to include column names to be clear about what's what). This would be useful to those who just want to run their FAM files as is.

One thing that I'm less clear on is the usefulness of adding population IDs. Is this just useful for the founders so we know which population to assign them beyond the known-pedigree? In this case, we may be able to include another column `is_founder` for individuals without known-parents. This flag could then be used to help users identify which individuals they may want to add population labels to (if they wanted to). Although I still think that population labels for fancy out-of-pedigree demographics would make sense to reserve for the PedigreeBuilder API. 
> One thing that I'm less clear on is the usefulness of adding population IDs.

People may want to do this sometimes, but it's not more useful than other sorts of metadata (birth location; name; sptial coordinates; age at death...). Of course, it's easy to just stick `0` in for the population column, as for phenotype in the FAM files.

And - sorry, I've not paid enough attention - but how exactly is the `time` field used? [The docs](https://github.com/tskit-dev/msprime/blob/main/msprime/pedigrees.py#L105) say that time is in units of "generations in the past", but in practice this is ill-defined, as different paths up to the same ancestor can easily have different numbers of hops. How do you deal with this, @LukeAndersonTrocme (or do you have to?)? I'm guessing that what we'd really like for `time` is the individual's birth time, in units of (mean generation time) ago; and in practice what the times mean is that the number of mutations separating offspring from parents is proportional to the time difference?
> And - sorry, I've not paid enough attention - but how exactly is the time field used?

This is the time of the two nodes associated with the individual @petrelharp. It's used to (a) sort the individuals and (b) decide when stuff gets added in to the simulation. (We originally used topogical sort, but then had to sort the edges by time later anyway to make a valid ts, so it's simpler to sort by time up front.) It's far simpler for the msprime simulation code if time is a fixed thing, so for now we're making it a required input.
Thanks for the input @LukeAndersonTrocme. I want to decouple this from time imputation because (a) there will be cases where we  *do* know the times, and we want to support this; and (b) it's a non-trivial problem, which we don't need to solve for the first version. (Of course it's useful having an implementation that works for your application, but we have to think hard about what would be the right thing to do more generally, and have to think through all the corner cases where it might fail, write tests for these, etc etc etc.) 

I've made the issue #1916 for time imputation which hopefully we can do as a follow up once we have the basics working.
Would you mind digging out a few examples of the pedigrees that you've seen @LukeAndersonTrocme? I think that'll help us a lot trying to make a file format that's familiar (hah!) for people, and allows them to create input data with the minimal amount of file hacking. 

(Don't post anything confidential, obvs!)
>  other sorts of metadata (birth location; name; sptial coordinates; age at death...)

Is there really a need to handle metadata? I have lots of metadata for the pedigree I'm currently working on, but I don't usually keep track of it while running pedigree based statistics. In these cases I tend to work with just the three columns (IID, MID, FID). It's easy to link to the metadata later with some information being specific to an individual (marriage date, marriage location) and other information specific to locations and regions.

Suggestion: let the user define the set of probands at generation 0. 
This information can be used to initialize the tree climbing algorithms that define generation time. This is also useful for working with a subset of individuals (e.g. the ascending genealogy form all probands in town_A). 

> different paths up to the same ancestor can easily have different numbers of hops. How do you deal with this?

Given a set of probands, we use the maximum genealogical depth for each historical individual. For our purposes we wanted to ensure that we traversed all children before reaching their parents. (code here: https://github.com/tskit-dev/msprime/issues/1916#issuecomment-968924261 )

> (a) there will be cases where we do know the times, and we want to support this

For example, the French-Canadian genealogy has two times for each individual: marriage date and parental marriage date

> (b) it's a non-trivial problem

Indeed, and the definition we're working with is just to help us sort the tree (backwards in time). Would this definition be useful for the msprime simulation code? Or would this need to be a forwards in time definition?
> Would you mind digging out a few examples of the pedigrees that you've seen @LukeAndersonTrocme? I think that'll help us a lot trying to make a file format that's familiar (hah!) for people, and allows them to create input data with the minimal amount of file hacking.
> 
> (Don't post anything confidential, obvs!)

Starting with the [BALSAC](https://ehps-net.eu/article/overview-balsac-population-database-past-developments-current-state-and-future-prospects/2020_vezinabournival_onlinefirst.pdf) genealogy because this is the one I've been working with the most. It comes with a few extra columns of metadata. 

Here's a toy example (I've manually changed these values so its not confidential)
```
ind,father,mother,sex,datemp,saisonmp,lieump,geo_mp,datem,saisonm,lieum,geo_m
2999,1368,1450,2,1871,2,1631,POINT(-73.3522 46.0952),1796,2,1269,POINT(-73.5606 45.5426)
3081,8780,8772,1,1866,1,1269,POINT(-73.5606 45.5426),1798,1,1269,POINT(-73.5606 45.5426)
3073,87877,7869,2,1882,1,1269,POINT(-73.5606 45.5426),1798,1,1269,POINT(-73.5606 45.5426)
3065,49212,49204,1,1846,2,1025,POINT(-72.43 46.3338),1783,1,1269,POINT(-73.5606 45.5426)
3057,6846,6838,2,1839,3,1269,POINT(-73.5606 45.5426),1783,1,1269,POINT(-73.5606 45.5426)
```
`ind,father,mother,sex` are fairly standard columns, where each individual `ind` has a `mother` and `father` who are also in the `ind` column.   

`datemp,saisonmp,lieump,geo_mp` these are the metadata columns for the parents (`mp` refers to parental marriage) where date, saison, lieu and geo refer to the marriage date, marriage season, marriage location_ID and marriage location lat/lon respectively.  The same set of columns are then repeated for the individuals marriage (`datem,saisonm,lieum,geo_m`)

Over the course of my project I've added more metadata to this, mostly pertaining to the geographical information (subregion, region, super-region) but as I mentioned in my previous [comment](https://github.com/tskit-dev/msprime/issues/1915#issuecomment-969040535), I don't need all the metadata all the time. Especially if I have a key to link back to. 

(I'll post a separate comment about genealogies I've seen in other papers)
The first example that came to mind was [this paper] that uses a very basic definition. (https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1008638#pcbi.1008638.s005)
![image](https://user-images.githubusercontent.com/12281604/141816503-27f23ba7-bc4d-426e-ab77-ac8a10bb369c.png)


Then I took a look at this [book about biobanks](https://link.springer.com/book/10.1007/978-1-59745-423-0) that has a section about genealogies (Ch. 10-11). I'll note that I think this information is from the perspective of biobanking (i.e information collecting) rather than from the perspective of research. I would expect that some of this information would not be privy to research scientists. 

Here's the contents of the Swedish Multi-generation Register:
```
National registration number of index person 
Sex of index person
Country of birth of index person 
Registration number of biological father 
Date of birth of biological father
Country of birth of biological father 
Registration number of biological mother 
Date of birth of biological mother 
Country of birth of biological mother 
Child’s position in the mother’s family 
Number of children – mother
Number of children – father 
Registration number of adoptive father 
Date of birth of adoptive father 
Country of birth of adoptive father 
Registration number of adoptive mother 
Date of birth of adoptive mother 
Country of birth of adoptive mother 
Date of immigration
Date of adoption
```
And here's the Icelandic Genealogy Database :
```
Date of birth (yymmdd) 
Place of birth
Twinning
Legitimacy
Date of death
Identification number of father 
Date of birth of father (yymmdd) 
Name of father
Identification number of mother 
Date of birth of mother (yymmdd) 
Name of mother
Note
Other identification
Unique personal identification number 
Name
Gender
```


### In Conclusion

I still feel like handling metadata is not necessary and may be confusing for the average user. My opinion would be to restrict the simulation to only accept (a) a three column pedigree and (b) an input set of probands. Users may want to use metadata before the simulations (to identify a specific set of probands as input) or afterwards (to compare the simulated genomes of sets of individuals). But AFAIK, both these cases don't require the need to handle metadata inside the msprime simulations.  
Barebones pedigree is often stored as

IID FID MID

sometimes FAMILY, SEX, BIRTH_DATE, BIRTH_LOCATION, POPULATION, etc. are provided in quite different formats (most don't vary much), so I would not worry about handling them - it's too big of a job and you only want to provide a way to pass a pedigree to msprime, so people should be able to cast their pedigree into a minimal format for msprime.
OK, how about
```
# IID FID MID TIME POPULATION METADATA
```
- We require the header
- Whitespace separated
- TIME/POPULATION/METADATA are optional (but TIME is required for now)
- METADATA must be the last column (can contain whitespace)
- ID values are arbitrary strings (not including whitespace, I guess)
LGTM!
OK, I'll code this up when I get a chance.

Fix patch:
diff --git a/docs/_toc.yml b/docs/_toc.yml
index 5b6278b99..18982cefd 100644
--- a/docs/_toc.yml
+++ b/docs/_toc.yml
@@ -18,6 +18,7 @@ parts:
 - caption: Utilities
   chapters:
   - file: rate_maps
+  - file: pedigrees
   - file: likelihoods
   - file: logging
 - caption: Miscellaneous
diff --git a/docs/ancestry.md b/docs/ancestry.md
index 041010f84..9d9079170 100644
--- a/docs/ancestry.md
+++ b/docs/ancestry.md
@@ -19,8 +19,28 @@ import msprime
 from IPython.display import SVG, set_matplotlib_formats
 from matplotlib import pyplot as plt
 import numpy as np
+import networkx as nx
 
 set_matplotlib_formats("svg")
+
+
+# Note: duplicating this from the version used in the pedigrees.md file.
+def draw_pedigree(ped_ts):
+
+    G = nx.DiGraph()
+    for ind in ped_ts.individuals():
+        time = ped_ts.node(ind.nodes[0]).time
+        pop = ped_ts.node(ind.nodes[0]).population
+        G.add_node(ind.id, time=time, population=pop)
+        for p in ind.parents:
+            if p != tskit.NULL:
+                G.add_edge(ind.id, p)
+    pos = nx.multipartite_layout(G, subset_key="time", align="horizontal")
+    colours = plt.rcParams['axes.prop_cycle'].by_key()['color']
+    node_colours = [colours[node_attr["population"]] for node_attr in G.nodes.values()]
+    nx.draw_networkx(G, pos, with_labels=True, node_color=node_colours)
+    plt.show()
+
 ```
 
 (sec_ancestry)=
@@ -81,6 +101,9 @@ this.
 {class}`.DiscreteTimeWrightFisher`
 : Generation-by-generation Wright-Fisher
 
+{class}`.FixedPedigree`
+: Simulated conditioned on a given pedigree
+
 {class}`.BetaCoalescent`
 : Beta coalescent multiple-merger
 
@@ -1918,6 +1941,7 @@ See the {ref}`sec_ancestry_multiple_chromosomes` section for information
 on how to simulate multiple chromosomes using the DTWF.
 :::
 
+
 (sec_ancestry_models_multiple_mergers)=
 
 ### Multiple merger coalescents
@@ -2306,6 +2330,279 @@ multiple models. See the {ref}`sec_logging` section for more
 details.
 :::
 
+(sec_ancestry_models_fixed_pedigree)=
+
+### Fixed pedigree
+
+A pedigree defines how *individuals* relate to each other, and these individual
+relationships constrain how genomes are inherited. We can simulate ancestry
+conditioned on a fixed pedigree using the {class}`.FixedPedigree` ancestry
+model. In this model the transfer of genetic material from parent to
+child is determined by the pedigree. Recombination
+occurs randomly such that each of a diploid individual's two genomes
+is a mixture of one of it
+parents two genomes, but everything else is determined by the pedigree.
+
+:::{warning}
+It is important to note that there is no concept of a population genetic
+model "outside" of the pedigree in this model, and that there
+are significant caveats that one needs to be aware of in terms
+of dealing with incomplete pedigree information and the treatment
+of founders.
+Please see the {ref}`sec_ancestry_models_fixed_pedigree_completing`
+section for more information.
+:::
+
+#### Simple example
+
+First we define a small pedigree, which follows a sample of three
+individuals through four generations:
+
+```{seealso}
+See the {ref}`sec_pedigrees` section for details on this pedigree
+text format and information on how pedigrees are encoded in msprime
+via the individual and node tables.
+```
+
+```{code-cell}
+import io
+ped_txt = """\
+# id parent0 parent1 time is_sample
+0   5   4   0.0 1
+1   3   3   0.0 1
+2   3   4   0.0 1
+3   7   7   1.0 0
+4   8   8   1.0 0
+5   8   6   1.0 0
+6   9   10  2.0 0
+7   10  10  2.0 0
+8   9   9   2.0 0
+9   11  12  3.0 0
+10  11  12  3.0 0
+11  NA  NA  4.0 0
+12  NA  NA  4.0 0
+"""
+pedigree = msprime.parse_pedigree(io.StringIO(ped_txt), sequence_length=100)
+
+draw_pedigree(pedigree.tree_sequence())
+```
+:::{note} See the {ref}`sec_pedigrees_visualisation` section for a
+definition of the ``draw_pedigree`` function used here.
+:::
+
+We then use this pedigree information as the
+{ref}`initial state<sec_ancestry_initial_state>` for the simulation
+(note that we must set the sequence length on the pedigree tables
+before calling ``sim_ancestry``):
+
+```{code-cell}
+ped_ts = msprime.sim_ancestry(
+    initial_state=pedigree, model="fixed_pedigree", random_seed=41,
+    recombination_rate=0.001)
+node_labels = {node.id: f"{node.individual}({node.id})" for node in ped_ts.nodes()}
+SVG(ped_ts.draw_svg(y_axis=True,  node_labels=node_labels, size=(600,200)))
+```
+
+The output tree sequence contains the simulated genetic ancestry conditioned
+on this input pedigree, as shown in the trees above. Each node is labelled
+with it the individual ID it's associated with, and its ID.
+The simulation stops at time 4, and we are left with multiple
+roots, which correspond to the founder individuals. If you only
+need simulations of the genetic ancestry through the pedigree,
+then these trees can be used perfectly well in downstream
+applications (such as {ref}`sec_mutations`).
+See the {ref}`sec_ancestry_models_fixed_pedigree_completing`
+section for information (and caveats about) completing these simulations.
+
+:::{note}
+Because the input pedigree fully describes the simulation many features such as
+{ref}`demography<sec_demography>` are not available when we use the {class}`.FixedPedigree`
+model. The {class}`.FixedPedigree` model also cannot be combined with other
+{ref}`ancestry models<sec_ancestry_models>`.
+:::
+
+:::{warning}
+The output tree sequence also contains the full pedigree information
+provided as input. This may be important if you are performing
+simulations on pedigrees that are not freely available.
+:::
+
+It is straightforward to remove the pedigree information from an
+output tree sequence, if required:
+
+```{code-cell}
+tables = ped_ts.dump_tables()
+tables.individuals.clear()
+tables.nodes.individual = np.full_like(tables.nodes.individual, tskit.NULL)
+censored_ts = tables.tree_sequence()
+censored_ts
+```
+
+(sec_ancestry_models_fixed_pedigree_completing)=
+
+#### Completing fixed pedigree simulations
+
+If you wish to simulate ancestry older than what is defined in the
+pedigree, you can "complete" (or "recapitate") the within-pedigree
+simulation using the standard methods described in the
+{ref}`sec_ancestry_initial_state` section.
+
+Continuing the example in the previous section, we provide the simulated
+ancestry in ``ped_ts`` as the initial state for a {class}`.DiscreteTimeWrightFisher`
+simulation with a population size of 10:
+
+```{code-cell}
+final_ts = msprime.sim_ancestry(
+    initial_state=ped_ts, model="dtwf", population_size=10, random_seed=42)
+SVG(final_ts.draw_svg(y_axis=True,  size=(600,200)))
+```
+
+We can see that the nodes 22 and 23 representing the founder genomes
+from the input pedigree are used as the starting point for the DTWF simulation,
+which then finishes at generation 6 when the MRCA (node 26) is found.
+
+(sec_ancestry_models_fixed_pedigree_missing_data)=
+
+#### Missing data
+
+The previous example is very straightforward because the pedigree is
+complete and we can trace all individuals back to the same
+generation in the pedigree. The pedigrees we have in practice
+are usually not like this, in that there is a great deal of missing data.
+The consequence of this for the {class}`.FixedPedigree` simulations
+is that genetic material will often be at a "dead end" at many
+levels of the pedigree during a simulation.
+
+```{code-cell}
+ped_txt = """\
+# id parent0 parent1 time is_sample
+0   5   4   0.0 1
+1   3   3   0.0 1
+2   3   4   0.0 1
+3   7   7   1.0 0
+4   8   8   1.0 0
+5   8   6   1.0 0
+6   9   10  2.0 0
+7   NA  NA  2.0 0
+8   9   9   2.0 0
+9   11  12  3.0 0
+10  11  12  3.0 0
+11  NA  NA  4.0 0
+12  NA  NA  4.0 0
+"""
+pedigree = msprime.parse_pedigree(io.StringIO(ped_txt), sequence_length=100)
+
+draw_pedigree(pedigree.tree_sequence())
+```
+
+This pedigree has a dead end for individual 7, which is
+reflected in the trees that we output from the {class}`.FixedPedigree` simulation:
+
+```{code-cell}
+ped_ts = msprime.sim_ancestry(
+    initial_state=pedigree, model="fixed_pedigree", random_seed=41,
+    recombination_rate=0.001)
+node_labels = {node.id: f"{node.individual}({node.id})" for node in ped_ts.nodes()}
+SVG(ped_ts.draw_svg(y_axis=True,  node_labels=node_labels, size=(600,200)))
+```
+
+When we provide this tree sequence as the
+{ref}`initial state<sec_ancestry_initial_state>` of another simulation, we
+then effectively start the simulation at the earlier time point corresponding to
+individual 7. Ancestral material in older parts of the pedigree is then
+brought into the simulation in the same manner as
+{ref}`ancient samples<sec_ancestry_samples_sampling_time>` as we proceed backwards
+in time.
+
+:::{warning}
+It is important to note that although ancestral material for pedigree
+founders is introduced into the recapitation simulation at the correct time,
+no account is taken of the number of lineages present in the pedigree
+when calculating population sizes. Thus, the pedigree must be seen as
+entirely **external** to the population model simulated during recapitation.
+:::
+
+(sec_ancestry_models_fixed_pedigree_demography)=
+
+#### Pedigrees and demography
+
+For complex simulations in which our pedigree individuals (and
+founders, in particular) are drawn from different populations,
+we will want our follow-up recapitation simulations to use
+a {ref}`demographic model<sec_demography>`. In order to do this,
+we must define the demographic model *before* creating
+the input pedigree. For example, consider the following
+simple demography:
+
+```{code-cell}
+demography = msprime.Demography()
+demography.add_population(name="A", initial_size=10)
+demography.add_population(name="B", initial_size=20)
+demography.add_population(name="C", initial_size=100)
+demography.add_population_split(time=10, derived=["A", "B"], ancestral="C");
+```
+
+We when define a pedigree in which the individuals belong
+to the two "leaf" populations:
+
+
+```{code-cell}
+ped_txt = """\
+# id    parent0 parent1 time    is_sample   population
+0   9   11  0.0 1   A
+1   10  10  0.0 1   A
+2   6   7   0.0 1   A
+3   7   7   0.0 1   A
+4   8   9   0.0 1   A
+5   10  10  0.0 1   A
+6   14  14  1.0 0   A
+7   12  12  1.0 0   A
+8   16  17  1.0 0   A
+9   12  13  1.0 0   A
+10  14  12  1.0 0   A
+11  15  16  1.0 0   A
+12  NA  NA  2.0 0   A
+13  NA  NA  2.0 0   A
+14  NA  NA  2.0 0   A
+15  NA  NA  2.0 0   A
+16  NA  NA  2.0 0   A
+17  NA  NA  2.0 0   A
+18  23  22  0.0 1   B
+19  22  22  0.0 1   B
+20  21  22  0.0 1   B
+21  25  25  1.0 0   B
+22  24  25  1.0 0   B
+23  24  24  1.0 0   B
+24  27  26  2.0 0   B
+25  27  26  2.0 0   B
+26  30  28  3.0 0   B
+27  30  29  3.0 0   B
+28  31  32  4.0 0   B
+29  31  32  4.0 0   B
+30  32  31  4.0 0   B
+31  NA  NA  5.0 0   B
+32  NA  NA  5.0 0   B
+"""
+pedigree = msprime.parse_pedigree(
+    io.StringIO(ped_txt), demography=demography, sequence_length=100)
+
+draw_pedigree(pedigree.tree_sequence())
+```
+
+
+
+```{code-cell}
+ped_ts = msprime.sim_ancestry(
+    initial_state=pedigree, model="fixed_pedigree", random_seed=41)
+node_labels = {node.id: f"{node.individual}({node.id})" for node in ped_ts.nodes()}
+SVG(ped_ts.draw_svg(y_axis=True,  node_labels=node_labels, size=(600,200)))
+```
+
+:::{todo}
+This section is incomplete. https://github.com/tskit-dev/msprime/issues/2009
+:::
+
+
 (sec_ancestry_missing_data)=
 
 ## Missing data
diff --git a/docs/api.md b/docs/api.md
index 70252e159..1287dd8c8 100644
--- a/docs/api.md
+++ b/docs/api.md
@@ -24,6 +24,7 @@ for discussion and examples of individual features.
   SmcApproxCoalescent
   SmcPrimeApproxCoalescent
   DiscreteTimeWrightFisher
+  FixedPedigree
   BetaCoalescent
   DiracCoalescent
   SweepGenicSelection
@@ -63,6 +64,9 @@ for discussion and examples of individual features.
 .. autosummary::
 
   RateMap
+  PedigreeBuilder
+  parse_pedigree
+
 
 ```
 
@@ -103,6 +107,11 @@ for discussion and examples of individual features.
 .. autoclass:: msprime.DiscreteTimeWrightFisher
 ```
 
+
+```{eval-rst}
+.. autoclass:: msprime.FixedPedigree
+```
+
 ```{eval-rst}
 .. autoclass:: msprime.BetaCoalescent
 ```
@@ -253,3 +262,16 @@ The node was created by a census event. Please see the
 .. autofunction:: msprime.log_mutation_likelihood
 ```
 
+### Pedigrees
+
+```{eval-rst}
+.. autoclass:: msprime.PedigreeBuilder
+    :members:
+```
+
+```{eval-rst}
+.. autofunction:: msprime.parse_pedigree
+```
+
+```{eval-rst}
+.. autofunction:: msprime.write_pedigree
diff --git a/docs/pedigrees.md b/docs/pedigrees.md
new file mode 100644
index 000000000..a87e7f95c
--- /dev/null
+++ b/docs/pedigrees.md
@@ -0,0 +1,171 @@
+---
+jupytext:
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.12
+    jupytext_version: 1.9.1
+kernelspec:
+  display_name: Python 3
+  language: python
+  name: python3
+---
+
+
+```{code-cell}
+:tags: [remove-cell]
+
+import io
+import msprime
+import tskit
+from IPython.display import SVG, set_matplotlib_formats
+from matplotlib import pyplot as plt
+import numpy as np
+import networkx as nx
+
+set_matplotlib_formats("svg")
+```
+
+(sec_pedigrees)=
+
+# Pedigrees
+
+Pedigrees describe parent-offspring relationships between individual organisms,
+and can be provided as input to contrain simulations of genetic ancestry
+(see the {ref}`sec_ancestry_models_fixed_pedigree` section for details).
+In this section we describe the data structures used to encode pedigree
+data in msprime, and the utilities used to create input pedigrees.
+
+:::{todo}
+This page is incomplete. Needs another pass through.
+https://github.com/tskit-dev/msprime/issues/2008
+:::
+
+(sec_pedigrees_encoding)=
+
+## Pedigree encoding
+
+Msprime uses the tskit
+{ref}`node<tskit:sec_node_table_definition>` and
+{ref}`individual<tskit:sec_individual_table_definition>` tables to encode pedigree
+information.
+Each {ref}`individual<sec_data_model_definitions_individual>`
+is defined by a row in the individual table, and
+this row index is the individual's **id**. Pedigree relationships
+are defined using the ``parents`` column, which contains the IDs
+of an individual's parents. Further information is stored in the
+{ref}`nodes<sec_data_model_definitions_node>` associated with the individual.
+
+:::{note}
+It is important to note that sample status, time and population information
+is associated with an individual's nodes, which is an artefact of
+tskit's node-centric design.
+:::
+
+It is possible to specify a pedigree model directly using tskit APIs,
+but it is simpler to use the {class}`.PedigreeBuilder` utility
+class (or the simplified text {ref}`sec_pedigrees_file_format`).
+In the following example we build a simple trio:
+
+```{code-cell}
+pb = msprime.PedigreeBuilder()
+mom_id  = pb.add_individual(time=1)
+dad_id  = pb.add_individual(time=1)
+pb.add_individual(time=0, parents=[mom_id, dad_id], is_sample=True)
+pedigree = pb.finalise()
+# TODO replace with display(pedigree) when its implemented in tskit
+# https://github.com/tskit-dev/tskit/issues/2093
+print(pedigree)
+```
+
+The pedigree returned by the {meth}`~.PedigreeBuilder.finalise` method
+contains the pedigree information defined by the calls to
+{meth}`~.PedigreeBuilder.add_individual`. For this trio, we began by adding
+the parents, and because they are founders, we don't provide any information
+about their parents. Each call to returns the integer ID of
+newly added individual.
+
+### Requirements
+
+
+
+(sec_pedigrees_file_format)=
+
+## Simplified file format
+
+The methods described in {ref}`sec_pedigrees_encoding` are general and
+allow for arbitrary metadata to be associated with individuals. It
+is often convenient to work with a text based representation of
+the pedigree, which is supported by the {func}`.parse_pedigree`.
+
+### Basic structure
+
+```{code-cell}
+txt = """\
+# id parent0 parent1 time
+child mom dad 0
+mom NA NA 1
+dad NA NA 1
+"""
+pedigree = msprime.parse_pedigree(io.StringIO(txt))
+display(pedigree.individuals)
+display(pedigree.nodes)
+```
+
+
+(sec_pedigrees_visualisation)=
+
+## Visualising pedigrees
+
+```{code-cell}
+import io
+txt = """\
+# id parent0 parent1 time
+child mom dad 0
+mom NA NA 1
+dad NA NA 1
+"""
+pedigree = msprime.parse_pedigree(io.StringIO(txt))
+display(pedigree.individuals)
+
+display(pedigree.nodes)
+```
+
+```{code-cell}
+pedigree.sequence_length = 100
+ts = msprime.sim_ancestry(initial_state=pedigree, model="fixed_pedigree",
+    random_seed=1)
+
+SVG(ts.draw_svg())
+```
+
+
+
+## Visualising pedigrees
+
+```{code-cell}
+
+
+def draw_pedigree(ped_ts):
+
+    G = nx.DiGraph()
+    for ind in ped_ts.individuals():
+        time = ped_ts.node(ind.nodes[0]).time
+        G.add_node(ind.id, time=time)
+        for p in ind.parents:
+            if p != tskit.NULL:
+                G.add_edge(ind.id, p)
+    pos = nx.multipartite_layout(G, subset_key="time", align="horizontal")
+    nx.draw_networkx(G, pos, with_labels=True)
+
+tables = msprime.pedigrees.sim_pedigree(
+    population_size=20, end_time=3, num_samples=5, direction="backward", random_seed=1)
+tables.individuals
+tables.sequence_length = 1
+
+
+draw_pedigree(tables.tree_sequence())
+
+```
+
+
diff --git a/msprime/__init__.py b/msprime/__init__.py
index cebdc4b38..eeeb60037 100644
--- a/msprime/__init__.py
+++ b/msprime/__init__.py
@@ -57,7 +57,7 @@
     Population,
 )
 
-from msprime.pedigrees import PedigreeBuilder
+from msprime.pedigrees import PedigreeBuilder, parse_pedigree, write_pedigree
 from msprime.intervals import RateMap
 from msprime.likelihood import log_arg_likelihood, log_mutation_likelihood
 
@@ -151,9 +151,11 @@
     "log_arg_likelihood",
     "log_mutation_likelihood",
     "mutate",
+    "parse_pedigree",
     "sim_ancestry",
     "sim_mutations",
     "simulate",
+    "write_pedigree",
 ]
 
 
diff --git a/msprime/pedigrees.py b/msprime/pedigrees.py
index febb334ee..3a019b8d0 100644
--- a/msprime/pedigrees.py
+++ b/msprime/pedigrees.py
@@ -1,5 +1,5 @@
 #
-# Copyright (C) 2019-2021 University of Oxford
+# Copyright (C) 2019-2022 University of Oxford
 #
 # This file is part of msprime.
 #
@@ -19,12 +19,190 @@
 """
 Pedigree utilities.
 """
+from __future__ import annotations
+
+import dataclasses
+
 import numpy as np
 import tskit
 
 from msprime import demography as demog_mod
 
 
+@dataclasses.dataclass
+class Column:
+    name: str
+
+    def parse_value(self, value: str) -> str:
+        return value
+
+
+class IdColumn(Column):
+    pass
+
+
+class NumericColumn(Column):
+    def parse_value(self, value: str) -> float:
+        return float(value)
+
+
+class BooleanColumn(Column):
+    def parse_value(self, value: str) -> bool:
+        return bool(int(value))
+
+
+class StringOrIntColumn(Column):
+    def parse_value(self, value: str) -> str | int:
+        try:
+            return int(value)
+        except ValueError:
+            pass
+        return value
+
+
+@dataclasses.dataclass
+class TextFileIndividual:
+    id: str  # noqa: A003
+    parent0: str | None
+    parent1: str | None
+    time: float
+    is_sample: bool | None = None
+    population: int = 0
+
+
+def _parse_pedigree_header(
+    header: str, demography: None | demog_mod.Demography
+) -> list[Column]:
+    tokens = header.split()
+    if tokens[0] != "#":
+        raise ValueError("First line must be the header and start with #")
+    required_cols = [
+        IdColumn("id"),
+        IdColumn("parent0"),
+        IdColumn("parent1"),
+        NumericColumn("time"),
+    ]
+    demography = (
+        demog_mod.Demography.isolated_model([1]) if demography is None else demography
+    )
+    optional_cols = [
+        BooleanColumn("is_sample"),
+        # There is the potential for ambiguity here in the population column
+        # if the string name of the population was an integer not equal to
+        # its ID. However, the Demography class requires that population names
+        # are valid Python identifiers, so this can't happen in practice.
+        StringOrIntColumn("population"),
+    ]
+    col_names = [col.name for col in required_cols]
+    tokens = tokens[1:]
+    if tokens[: len(required_cols)] != col_names:
+        raise ValueError(f"The {col_names} columns are required")
+
+    columns = list(required_cols)
+    optional_col_map = {col.name: col for col in optional_cols}
+    for colname in tokens[len(required_cols) :]:
+        if colname not in optional_col_map:
+            raise ValueError(f"Column '{colname}' not supported")
+        columns.append(optional_col_map[colname])
+    return columns
+
+
+def parse_pedigree(text_file, demography=None, sequence_length=None):
+    """
+    Parse a text describing a pedigree used for input to the
+    :class:`.FixedPedigree` ancestry model.
+    """
+    # First line must be the header and contain the column headers
+    header = next(text_file, None)
+    if header is None:
+        raise ValueError("Pedigree file must contain at least a header")
+    columns = _parse_pedigree_header(header, demography=demography)
+    individuals = []
+    tsk_id_map = {"NA": -1}
+    for tsk_id, row in enumerate(text_file):
+        line = tsk_id + 2
+        tokens = row.split()
+        if len(tokens) != len(columns):
+            raise ValueError(
+                "Incorrect number of columns "
+                f"(found={len(tokens)} need={len(columns)}) at line {line}: {row}"
+            )
+        kwargs = {}
+        for value, column in zip(tokens, columns):
+            try:
+                kwargs[column.name] = column.parse_value(value)
+            except ValueError as ve:
+                raise ValueError(
+                    f"Error parsing value for column '{column.name}' "
+                    f"on line {line}: {ve}"
+                )
+        ind = TextFileIndividual(**kwargs)
+        if ind.id == "NA":
+            raise ValueError(f"'NA' cannot be used as an individual ID (line {line})")
+        if ind.id in tsk_id_map:
+            raise ValueError(f"Duplicate ID at line {line}")
+        tsk_id_map[ind.id] = len(individuals)
+        individuals.append(ind)
+
+    # Convert the list of individuals to a tskit pedigree.
+    builder = PedigreeBuilder(
+        individuals_metadata_schema=tskit.MetadataSchema.permissive_json(),
+        demography=demography,
+    )
+    for ind in individuals:
+        parents = []
+        for parent in [ind.parent0, ind.parent1]:
+            if parent not in tsk_id_map:
+                raise ValueError(
+                    f"Parent ID '{parent}' not defined (line {tsk_id_map[ind.id] + 2})"
+                )
+            parents.append(tsk_id_map[parent])
+        builder.add_individual(
+            parents=parents,
+            time=ind.time,
+            is_sample=ind.is_sample,
+            metadata={"file_id": ind.id},
+            population=ind.population,
+        )
+    return builder.finalise(sequence_length)
+
+
+def write_pedigree(ts, out):
+
+    print("# id\tparent0\tparent1\ttime\tis_sample\tpopulation", file=out)
+    for ind in ts.individuals():
+        if len(ind.nodes) != 2:
+            raise ValueError(
+                "Invalid pedigree format: each individual must be associated with "
+                "two nodes"
+            )
+        nodes = [ts.node(node) for node in ind.nodes]
+        time = {node.time for node in nodes}
+        is_sample = {node.is_sample() for node in nodes}
+        population = {node.population for node in nodes}
+        if len(time) != 1:
+            raise ValueError("Pedigree individuals must have the same node time")
+        if len(is_sample) != 1:
+            raise ValueError("Pedigree individuals must have the same sample status")
+        if len(population) != 1:
+            raise ValueError("Pedigree individuals must have the same node population")
+
+        time = time.pop()
+        is_sample = is_sample.pop()
+        population = population.pop()
+        parents = []
+        for parent in ind.parents:
+            if parent == tskit.NULL:
+                parents.append("NA")
+            else:
+                parents.append(parent)
+        print(
+            f"{ind.id}\t{parents[0]}\t{parents[1]}\t{time}\t{is_sample}\t"
+            f"{population}",
+            file=out,
+        )
+
+
 class PedigreeBuilder:
     def __init__(self, demography=None, individuals_metadata_schema=None):
         if demography is None:
@@ -41,8 +219,17 @@ def __init__(self, demography=None, individuals_metadata_schema=None):
             self.individuals.metadata_schema = individuals_metadata_schema
 
     def add_individual(
-        self, *, time, is_sample=None, parents=None, population=None, metadata=None
+        self,
+        *,
+        time,
+        is_sample=None,
+        parents=None,
+        population=None,
+        metadata=None,
     ):
+        """
+        Adds the specified individual, returning its ID.
+        """
         if is_sample is None:
             # By default, individuals at time 0 are marked as samples
             is_sample = time == 0
@@ -91,12 +278,16 @@ def add_individuals(self, *, parents, time):
         )
         return ind_ids
 
-    def finalise(self, sequence_length):
+    def finalise(self, sequence_length=None):
+        """
+        Returns the table collection representing the defined pedigree.
+        """
         copy = self.tables.copy()
-        copy.sequence_length = sequence_length
+        copy.sequence_length = 1
         # Not strictly necessary, but it's useful when the tables are written
         # to file and are loaded as as a tree sequence.
         copy.build_index()
+        copy.sequence_length = -1 if sequence_length is None else sequence_length
         return copy
 
 
@@ -189,61 +380,3 @@ def sim_pedigree(
 
     tables.sequence_length = -1 if sequence_length is None else sequence_length
     return tables
-
-
-def parse_fam(fam_file):
-    """
-    Parse PLINK .fam file and convert to tskit IndividualTable.
-    Assumes fam file contains five columns: FID, IID, PAT, MAT, SEX
-    :param fam_file: PLINK .fam file object
-    :param tskit.TableCollection tc: TableCollection with IndividualTable to
-        which the individuals will be added
-    """
-    individuals = np.loadtxt(
-        fname=fam_file,
-        dtype=str,
-        ndmin=2,  # read file as 2-D table
-        usecols=(0, 1, 2, 3, 4),  # only keep FID, IID, PAT, MAT, SEX columns
-    )  # requires same number of columns in each row, i.e. not ragged
-
-    id_map = {}  # dict for translating PLINK ID to tskit IndividualTable ID
-    for tskit_id, (plink_fid, plink_iid, _pat, _mat, _sex) in enumerate(individuals):
-        # include space between strings to ensure uniqueness
-        plink_id = f"{plink_fid} {plink_iid}"
-        if plink_id in id_map:
-            raise ValueError("Duplicate PLINK ID: {plink_id}")
-        id_map[plink_id] = tskit_id
-    id_map["0"] = -1  # -1 is used in tskit to denote "missing"
-
-    tc = tskit.TableCollection(1)
-    tb = tc.individuals
-    tb.metadata_schema = tskit.MetadataSchema(
-        {
-            "codec": "json",
-            "type": "object",
-            "properties": {
-                "plink_fid": {"type": "string"},
-                "plink_iid": {"type": "string"},
-                "sex": {"type": "integer"},
-            },
-            "required": ["plink_fid", "plink_iid", "sex"],
-            "additionalProperties": True,
-        }
-    )
-    for plink_fid, plink_iid, pat, mat, sex in individuals:
-        sex = int(sex)
-        if not (sex in range(3)):
-            raise ValueError(
-                "Sex must be one of the following: 0 (unknown), 1 (male), 2 (female)"
-            )
-        metadata_dict = {"plink_fid": plink_fid, "plink_iid": plink_iid, "sex": sex}
-        pat_id = f"{plink_fid} {pat}" if pat != "0" else pat
-        mat_id = f"{plink_fid} {mat}" if mat != "0" else mat
-        tb.add_row(
-            parents=[
-                id_map[pat_id],
-                id_map[mat_id],
-            ],
-            metadata=metadata_dict,
-        )
-    return tb
diff --git a/requirements/CI-docs/requirements.txt b/requirements/CI-docs/requirements.txt
index 32b36ebf9..7a581416f 100644
--- a/requirements/CI-docs/requirements.txt
+++ b/requirements/CI-docs/requirements.txt
@@ -2,6 +2,7 @@ daiquiri==3.0.1
 demes==0.2.1
 demesdraw==0.2.0
 jupyter-book==0.12.1
+networkx==2.6.3
 matplotlib==3.5.0
 newick==1.3.1
 sphinx==4.4.0
diff --git a/requirements/development.txt b/requirements/development.txt
index 61d4111af..f8ee2c748 100644
--- a/requirements/development.txt
+++ b/requirements/development.txt
@@ -25,6 +25,7 @@ sphinx>=4.4
 sphinx-argparse
 sphinx-issues
 jupyter-book
+networkx
 
 newick>=1.3.0
 # https://github.com/tskit-dev/msprime/issues/1860


Test patch:
diff --git a/tests/test_pedigree.py b/tests/test_pedigree.py
index d6b2710a9..1db273683 100644
--- a/tests/test_pedigree.py
+++ b/tests/test_pedigree.py
@@ -1,7 +1,7 @@
 import collections
-import dataclasses
+import io
 import sys
-import tempfile
+import textwrap
 
 import numpy as np
 import pytest
@@ -1087,225 +1087,6 @@ def test_trio_parents_different_pops_with_migrations(self):
         assert ts.node(6).time > 2
 
 
-@dataclasses.dataclass
-class FamEntry:
-    fid: str = "0"
-    iid: str = "0"
-    pat: str = "0"
-    mat: str = "0"
-    sex: str = "0"
-    phen: str = None
-
-    def get_row(self, delimiter="\t"):
-        return delimiter.join(
-            [x for x in dataclasses.asdict(self).values() if x is not None]
-        )
-
-
-class TestParseFam:
-    """
-    Tests for the parse_fam function.
-    """
-
-    def get_parsed_fam(self, entries, delimiter="\t"):
-        content = "\n".join([entry.get_row(delimiter=delimiter) for entry in entries])
-        with tempfile.TemporaryFile() as f:
-            f.write(bytes(content, "utf-8"))
-            f.seek(0)
-            return pedigrees.parse_fam(f)
-
-    def test_empty_file(self):
-        entries = []
-        with pytest.warns(UserWarning):
-            tb = self.get_parsed_fam(entries=entries)
-        assert len(tb) == 0
-
-    @pytest.mark.parametrize("iid", ["1", "a", "100", "abc"])
-    def test_single_line(self, iid):
-        entries = [FamEntry(iid=iid)]
-        tb = self.get_parsed_fam(entries=entries)
-        assert len(tb) == 1
-        assert np.array_equal(tb[0].parents, [-1, -1])
-        assert tb[0].metadata["plink_fid"] == "0"
-        assert tb[0].metadata["plink_iid"] == str(iid)
-        assert tb[0].metadata["sex"] == 0
-
-    @pytest.mark.parametrize("iids", [("1", "2"), ("a", "b")])
-    def test_multiple_line_file(self, iids):
-        # test both integer and string IIDs
-        iid1, iid2 = iids
-        entries = [FamEntry(iid=iid1), FamEntry(iid=iid2)]
-        tb = self.get_parsed_fam(entries=entries)
-        assert len(tb) == 2
-        for idx in range(2):
-            assert np.array_equal(tb[idx].parents, [-1, -1])
-            assert tb[idx].metadata["plink_fid"] == "0"
-            assert tb[idx].metadata["plink_iid"] == str(entries[idx].iid)
-            assert tb[idx].metadata["sex"] == 0
-
-    @pytest.mark.parametrize("n_cols", range(1, 5))
-    def test_insufficient_cols(self, n_cols):
-        fields = list(dataclasses.asdict(FamEntry()))
-        entry = FamEntry(iid="1")
-        for field in fields[n_cols:]:
-            entry.__setattr__(field, None)
-        with pytest.raises(IndexError):
-            self.get_parsed_fam(entries=[entry])
-
-    def test_unrelated_duplicate_iids(self):
-        # Individuals have the same IID, but are in different families
-        entries = [FamEntry(iid="1"), FamEntry(fid="1", iid="1")]
-        tb = self.get_parsed_fam(entries=entries)
-        assert len(tb) == 2
-        assert tb[0].metadata["plink_fid"] == "0"
-        assert tb[0].metadata["plink_iid"] == "1"
-        assert tb[1].metadata["plink_fid"] == "1"
-        assert tb[1].metadata["plink_iid"] == "1"
-
-    def test_duplicate_rows(self):
-        entries = [FamEntry(iid="1"), FamEntry(iid="1")]
-        with pytest.raises(ValueError):
-            self.get_parsed_fam(entries=entries)
-
-    def test_space_delimited(self):
-        entries = [FamEntry(iid="1")]
-        tb = self.get_parsed_fam(entries=entries, delimiter=" ")
-        assert np.array_equal(tb[0].parents, [-1, -1])
-        assert tb[0].metadata["plink_fid"] == "0"
-        assert tb[0].metadata["plink_iid"] == "1"
-        assert tb[0].metadata["sex"] == 0
-
-    def test_missing_phen_col(self):
-        entries = [FamEntry(iid="1", phen="1")]
-        tb = self.get_parsed_fam(entries=entries)
-
-        entries = [FamEntry(iid="1")]  # remove last column (PHEN column)
-        tb_missing = self.get_parsed_fam(entries=entries)
-
-        assert tb == tb_missing
-
-    @pytest.mark.parametrize("sex", [-2, 3, "F"])
-    def test_bad_sex_value(self, sex):
-        entries = [FamEntry(iid="1", sex=str(sex))]
-        with pytest.raises(ValueError):
-            self.get_parsed_fam(entries=entries)
-
-    def test_empty_sex_value(self):
-        entries = [FamEntry(iid="1", sex="")]
-        with pytest.raises(IndexError):
-            self.get_parsed_fam(entries=entries)
-
-    def test_single_family_map_parent_ids(self):
-        # PAT is mapped if the individual exists in the dataset
-        entries = [FamEntry(iid="1"), FamEntry(iid="2", pat="1")]
-        tb = self.get_parsed_fam(entries=entries)
-        assert np.array_equal(tb[1].parents, [0, -1])
-
-        # MAT is mapped if the individual exists in the dataset
-        entries = [FamEntry(iid="1"), FamEntry(iid="2", mat="1")]
-        tb = self.get_parsed_fam(entries=entries)
-        assert np.array_equal(tb[1].parents, [-1, 0])
-
-        # both parent IDs are remapped if the both parents exist in the dataset
-        entries = [
-            FamEntry(iid="1"),
-            FamEntry(iid="2"),
-            FamEntry(iid="3", pat="1", mat="2"),
-        ]
-        tb = self.get_parsed_fam(entries=entries)
-        assert np.array_equal(tb[2].parents, [0, 1])
-
-    def test_missing_parent_id(self):
-        # KeyError raised if at least one parent (PAT) does not exist in dataset
-        entries = [
-            FamEntry(iid="2"),
-            FamEntry(iid="3", pat="1", mat="2"),
-        ]
-        with pytest.raises(KeyError):
-            self.get_parsed_fam(entries=entries)
-
-        # KeyError raised if at least one parent (MAT) does not exist in dataset
-        entries = [
-            FamEntry(iid="1"),
-            FamEntry(iid="3", pat="1", mat="2"),
-        ]
-        with pytest.raises(KeyError):
-            self.get_parsed_fam(entries=entries)
-
-        # KeyError raised if both parents do not exist in dataset
-        entries = [FamEntry(iid="1", pat="2", mat="3")]
-        with pytest.raises(KeyError):
-            self.get_parsed_fam(entries=entries)
-
-    def test_multiple_family_map_parent_ids(self):
-        # parents mapped correctly when the same parent ID is used in different families
-        entries = [
-            FamEntry(iid="2"),
-            FamEntry(iid="1"),
-            FamEntry(fid="1", iid="2"),
-            FamEntry(fid="1", iid="1"),
-            FamEntry(iid="3", pat="1", mat="2"),
-            FamEntry(fid="1", iid="3", pat="1", mat="2"),
-        ]
-        tb = self.get_parsed_fam(entries=entries)
-        for idx in range(4):
-            assert np.array_equal(tb[idx].parents, [-1, -1])
-        assert np.array_equal(tb[4].parents, [1, 0])
-        assert np.array_equal(tb[5].parents, [3, 2])
-
-        # KeyError raised when FID does not match, even if parent ID matches
-        entries = [
-            FamEntry(iid="2"),
-            FamEntry(iid="1"),
-            FamEntry(iid="3", pat="1", mat="2"),
-            FamEntry(
-                fid="1", iid="1", pat="2", mat="3"
-            ),  # there is no parent with FID=1, IID=3
-            FamEntry(fid="1", iid="2"),
-        ]
-        with pytest.raises(KeyError):
-            self.get_parsed_fam(entries)
-
-    def test_grandparents(self):
-        entries = [
-            FamEntry(iid="4"),
-            FamEntry(iid="3"),
-            FamEntry(iid="2"),
-            FamEntry(iid="1"),
-            FamEntry(iid="6", pat="3", mat="4"),
-            FamEntry(iid="5", pat="1", mat="2"),
-            FamEntry(iid="7", pat="5", mat="6"),
-        ]
-        tb = self.get_parsed_fam(entries=entries)
-        assert np.array_equal(tb[4].parents, [1, 0])
-        assert np.array_equal(tb[5].parents, [3, 2])
-        assert np.array_equal(tb[6].parents, [5, 4])
-
-    def test_children_before_parents(self, tmp_path):
-        entries = [
-            FamEntry(iid="1", pat="2", mat="3"),
-            FamEntry(iid="2"),
-            FamEntry(iid="3"),
-        ]
-        content = "\n".join([entry.get_row() for entry in entries])
-        fam_path = f"{tmp_path}/test.fam"
-        with open(fam_path, "w+") as f:
-            f.write(content)
-            f.seek(0)
-            tb = pedigrees.parse_fam(f)
-
-        tc = tskit.TableCollection(1)
-        # Issue 1489 will make this better
-        tc.individuals.metadata_schema = tb.metadata_schema
-        for row in tb:
-            tc.individuals.append(row)
-        tc.tree_sequence()  # creating tree sequence should succeed
-
-        assert np.array_equal(tb[0].parents, [1, 2])
-        assert np.array_equal(tb[1].parents, [-1, -1])
-        assert np.array_equal(tb[2].parents, [-1, -1])
-
-
 class TestPedigreeBuilder:
     def test_is_sample_default(self):
         pb = pedigrees.PedigreeBuilder()
@@ -1471,6 +1252,281 @@ def test_indexed(self):
         assert t1.has_index()
 
 
+def parse_indented(text, demography=None):
+    return msprime.parse_pedigree(
+        io.StringIO(textwrap.dedent(text)), demography=demography
+    )
+
+
+class TestParsePedigreeExamples:
+    def test_single_all_columns(self):
+        text = """\
+        # id parent0 parent1 time is_sample
+        X NA NA 0 1
+        """
+        tables = parse_indented(text)
+        assert len(tables.individuals) == 1
+        assert len(tables.nodes) == 2
+        assert list(tables.individuals[0].parents) == [-1, -1]
+        assert tables.nodes[0].time == 0
+        assert tables.nodes[0].flags == tskit.NODE_IS_SAMPLE
+        assert tables.individuals[0].metadata == {"file_id": "X"}
+
+    def test_single_row_min_columns(self):
+        text = """\
+        # id parent0 parent1 time
+        X NA NA 0
+        """
+        tables = parse_indented(text)
+        assert len(tables.individuals) == 1
+        assert len(tables.nodes) == 2
+        assert list(tables.individuals[0].parents) == [-1, -1]
+        assert tables.nodes[0].time == 0
+        assert tables.nodes[0].flags == tskit.NODE_IS_SAMPLE
+        assert tables.individuals[0].metadata == {"file_id": "X"}
+
+    def test_no_rows(self):
+        text = """\
+        # id parent0 parent1 time is_sample
+        """
+        tables = parse_indented(text)
+        assert len(tables.individuals) == 0
+        assert len(tables.nodes) == 0
+
+    def test_trio(self):
+        text = """\
+        # id parent0 parent1 time
+        child mom dad 0
+        mom NA NA 1
+        dad NA NA 1
+        """
+        tables = parse_indented(text)
+        assert len(tables.individuals) == 3
+        assert len(tables.nodes) == 6
+        assert tables.individuals[0].metadata["file_id"] == "child"
+        assert list(tables.individuals[0].parents) == [1, 2]
+        assert tables.individuals[1].metadata["file_id"] == "mom"
+        assert list(tables.individuals[1].parents) == [-1, -1]
+        assert tables.individuals[2].metadata["file_id"] == "dad"
+        assert list(tables.individuals[2].parents) == [-1, -1]
+
+    def test_single_population_name(self):
+        text = """\
+        # id parent0 parent1 time population
+        A NA NA 0 pop_0
+        """
+        demography = msprime.Demography.isolated_model([1])
+        tables = parse_indented(text, demography=demography)
+        assert len(tables.individuals) == 1
+        assert tables.nodes[0].population == 0
+        assert tables.nodes[1].population == 0
+
+    def test_single_population_id(self):
+        text = """\
+        # id parent0 parent1 time population
+        A NA NA 0 0
+        """
+        demography = msprime.Demography.isolated_model([1])
+        tables = parse_indented(text, demography=demography)
+        assert len(tables.individuals) == 1
+        assert tables.nodes[0].population == 0
+        assert tables.nodes[1].population == 0
+
+    def test_two_populations_name(self):
+        text = """\
+        # id parent0 parent1 time population
+        A NA NA 0 pop_1
+        """
+        demography = msprime.Demography.isolated_model([1, 1])
+        tables = parse_indented(text, demography=demography)
+        assert len(tables.individuals) == 1
+        assert tables.nodes[0].population == 1
+        assert tables.nodes[1].population == 1
+
+    def test_two_populations_id(self):
+        text = """\
+        # id parent0 parent1 time population
+        A NA NA 0 1
+        """
+        demography = msprime.Demography.isolated_model([1, 1])
+        tables = parse_indented(text, demography=demography)
+        assert len(tables.individuals) == 1
+        assert tables.nodes[0].population == 1
+        assert tables.nodes[1].population == 1
+
+
+class TestParsePedigreeErrors:
+    def test_empty_file(self):
+        with pytest.raises(ValueError, match="at least a header"):
+            parse_indented("")
+
+    def test_bad_header(self):
+        with pytest.raises(ValueError, match="start with #"):
+            parse_indented("id")
+
+    @pytest.mark.parametrize("column", ["id", "parent0", "parent1", "time"])
+    def test_missing_required(self, column):
+        columns = [col for col in ["id", "parent0", "parent1", "time"] if col != column]
+        header = "# " + " ".join(columns)
+        with pytest.raises(ValueError, match="columns are required"):
+            parse_indented(header)
+
+    def test_unknown_column(self):
+        text = """\
+        # id parent0 parent1 time unknown_col
+        """
+        with pytest.raises(ValueError, match="'unknown_col' not supported"):
+            parse_indented(text)
+
+    def test_parital_row_min(self):
+        text = """\
+        # id parent0 parent1 time
+        child
+        """
+        with pytest.raises(ValueError, match="Incorrect number of columns"):
+            parse_indented(text)
+
+    def test_parital_row_missing_optional(self):
+        text = """\
+        # id parent0 parent1 time is_sample
+        child NA NA 0
+        """
+        with pytest.raises(ValueError, match="Incorrect number of columns"):
+            parse_indented(text)
+
+    def test_duplicate_id(self):
+        text = """\
+        # id parent0 parent1 time
+        child NA NA 0
+        child NA NA 1
+        """
+        with pytest.raises(ValueError, match="Duplicate ID"):
+            parse_indented(text)
+
+    def test_missing_id(self):
+        text = """\
+        # id parent0 parent1 time
+        child missing NA 0
+        """
+        with pytest.raises(ValueError, match="Parent ID 'missing' not defined"):
+            parse_indented(text)
+
+    def test_bad_time(self):
+        text = """\
+        # id parent0 parent1 time
+        child NA NA bad_time
+        """
+        with pytest.raises(ValueError, match="value for column 'time' on line 2"):
+            parse_indented(text)
+
+    def test_bad_is_sample(self):
+        text = """\
+        # id parent0 parent1 time is_sample
+        child NA NA 0 abcd
+        """
+        with pytest.raises(ValueError, match="value for column 'is_sample' on line 2"):
+            parse_indented(text)
+
+    def test_id_is_NA(self):
+        text = """\
+        # id parent0 parent1 time
+        NA NA NA 0
+        """
+        with pytest.raises(ValueError, match="'NA' cannot be used"):
+            parse_indented(text)
+
+    def test_line_no(self):
+        text = """\
+        # id parent0 parent1 time
+        me NA NA 0
+        NA NA NA 0
+        """
+        # People (and vim) usually count file lines from 1, so error on line 3
+        with pytest.raises(ValueError, match="line 3"):
+            parse_indented(text)
+
+
+class TestPedigreeTextRoundTrip:
+    def verify(self, pedigree, demography=None):
+        buff = io.StringIO()
+        pedigrees.write_pedigree(pedigree.tree_sequence(), buff)
+        buff.seek(0)
+        parsed = pedigrees.parse_pedigree(
+            buff, demography, sequence_length=pedigree.sequence_length
+        )
+        pedigree.assert_equals(parsed, ignore_metadata=True)
+
+    @pytest.mark.parametrize("direction", ["backward", "forward"])
+    def test_sim_pedigree(self, direction):
+        ped_tables = pedigrees.sim_pedigree(
+            population_size=5,
+            end_time=5,
+            random_seed=1234,
+            direction=direction,
+            sequence_length=1,
+        )
+        self.verify(ped_tables)
+
+    def test_simulate_pedigree(self):
+        ped_tables = simulate_pedigree(
+            num_founders=20,
+            num_generations=5,
+            sequence_length=1,
+        )
+        self.verify(ped_tables)
+
+    def test_two_populations(self):
+        t1 = pedigrees.sim_pedigree(
+            population_size=5,
+            end_time=3,
+            random_seed=1234,
+            sequence_length=1,
+        )
+        t1.populations.add_row({"name": "pop_1", "description": ""})
+        t2 = t1.copy()
+        t2.nodes.population = t2.nodes.population + 1
+        ped = join_pedigrees([t1, t2])
+        demography = msprime.Demography.isolated_model([1, 1])
+        self.verify(ped.tables, demography=demography)
+
+
+class TestWritePedigreeErrors:
+    @pytest.mark.parametrize("num_nodes", [0, 1, 3])
+    def test_incorrect_num_nodes(self, num_nodes):
+        tables = tskit.TableCollection(1)
+        tables.individuals.add_row(parents=[-1, -1])
+        for _ in range(num_nodes):
+            tables.nodes.add_row(time=0, individual=0)
+        with pytest.raises(ValueError, match="two nodes"):
+            pedigrees.write_pedigree(tables.tree_sequence(), io.StringIO())
+
+    def test_mismatched_node_times(self):
+        tables = tskit.TableCollection(1)
+        tables.individuals.add_row(parents=[-1, -1])
+        tables.nodes.add_row(flags=0, time=0, individual=0)
+        tables.nodes.add_row(flags=1, time=0, individual=0)
+        with pytest.raises(ValueError, match="same sample status"):
+            pedigrees.write_pedigree(tables.tree_sequence(), io.StringIO())
+
+    def test_mismatched_node_flags(self):
+        tables = tskit.TableCollection(1)
+        tables.individuals.add_row(parents=[-1, -1])
+        tables.nodes.add_row(time=0, individual=0)
+        tables.nodes.add_row(time=1, individual=0)
+        with pytest.raises(ValueError, match="same node time"):
+            pedigrees.write_pedigree(tables.tree_sequence(), io.StringIO())
+
+    def test_mismatched_node_populations(self):
+        tables = tskit.TableCollection(1)
+        tables.individuals.add_row(parents=[-1, -1])
+        tables.populations.add_row()
+        tables.populations.add_row()
+        tables.nodes.add_row(time=0, individual=0, population=0)
+        tables.nodes.add_row(time=0, individual=0, population=1)
+        with pytest.raises(ValueError, match="same node population"):
+            pedigrees.write_pedigree(tables.tree_sequence(), io.StringIO())
+
+
 if __name__ == "__main__":
     # Easy way to output a simulated pedigree for command line testing.
     # Run with python3 tests/test_pedigree.py NUM_FOUNDERS NUM_GENERATIONS > out.trees