<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Sibship clustering
Tom Ellis, March 2017, updated June 2020
Step1:  FAPS uses information in a paternityArray to generate plausible full-sibship configurations. This information is stored as a sibshipCluster object, and is the basis for almost all inference about biological processes in FAPS.
This notebook will examine how to
Step2:  We take the first individal as the mother and mate her to four males, to create three full sibships of five offspring. We then generate a paternityArray based on the genotype data.
Step3:  It is straightforward to cluster offspring into full sibships. For now we'll stick with the default number of Monte Carlo draws.
Step4:  The default number of Monte Carlo draws is 1000, which seems to work in most cases. I have dropped it to 100 in cases where I wanted to call sibship_clustering many times, such as in an MCMC loop, when finding every possible candidate wasn't a priority. You could also use more draws if you really wanted to be sure you had completely sampled the space of compatible candidate fathers. Speeds are unlikely to increase linearly with number of draws
Step5:  We discussed this in figure 5 of the FAPS paper should you be interested in more on this.
Sibling relationships
Sibship clustering calculates likelihoods that each pair of offspring are full siblings, then builds a dendrogram from this. We can visualise this dendrogram if we so wish, although the output is not pretty.
Step6:  Offspring individuals are labelled by their index in the array. Since full sibships are of size five we should expect to see clusters of {0,1,2,3,4}, {5,6,7,8,9} and {10,11,12,13,14}. This is indeed what we do see. What is difficult to see on the dendrogram are the branches connecting full siblings at the very bottom of the plot. If we bisect this dendrogram at different places on the y-axis we can infer different ways to partition the offspring into full siblings.
sc is an object of class sibshipCluster that contains various information about the array. Of primary interest are the set of partition structures inferred from the dendrogram. There are sixteen partitions - one for each individual in the array (i.e. one for each bifurcation in the dendrogram).
Step7:  What is key about partition structures is that each symbol represents a unique but arbitrary family identifier. For example in the third row we see the true partition structure, with individuals grouped into three groups of five individuals.
Step8:  Beyond denoting who is in a family with whom, the labels are arbitrary, with no meaningful order. This partition would be identical to [0,0,0,0,0,1,1,1,1,1,2,2,2,2,2] or [10,10,10,10,10,7,7,7,7,7,22,22,22,22,22] for example.
Each partition is associated with a log likelihood and equivalent log probability. We can see from both scores that the third partition is most consistent with the data. This is of course the true partition.
Step9:  We also see that the first and second partitions have non-zero, but small likelihoods. Parititons 5-8 have negative infinity log likelihood - they are incompatible with the data. These partitions split up true full siblings, and there is no way to reconcile this with the data. In real world situations such partitions might have non-zero likelihoods if they were an unrelated candidate male compatible with one or more offspring through chance alone.
In some cases there can be rounding error when log probabilities are exponentiated and probabilities do not sum to one. This is classic machine error, and the reason it is good to work with log values wherever possible. We can check
Step10:  You can directly call the most likely partition. This is somewhat against the spirit of fractional analyses though...
Step11:  For information about fine scale relationships, sc.full_sib_matrix() returns an $n*n$ matrix, where $n$ is the number of offspring. Each element describes the (log) probability that a pair of individuals are full siblings, averaged over partition structures and paternity configurations. If we plot this using a heatmap you can clearly see the five full sibships jump out as blocks of yellow (>90% probability of being full siblings) against a sea of purple (near zero probability of being full siblings).
Step12:  Note that real datasets seldom look this tidy!
Inferring family structure
For this section we will simulate a slightly more interesting family structure. This block of code creates a half-sib array of 15 offspring from five fathers, where each father contributes five, four, three, two and one offspring respectively. It then performs sibship clustering on the array. We use 1000 candidate males and 50 loci.
Step13:  Number of families
We saw before that we could call a list of valid partitions for sc using sc.partitions. The output is not terribly enlightening on its own, however. We could instead ask how probable it is that there are x full sibships in the array, integrating over all partition structures. Here each number is the probability that there are 1, 2, ..., 15 families.
Step14:  We could show the same information graphically. Its clear that almost all the probability denisty is around $x=5$ families.
Step15:  Family size
We can also get the distribution of family sizes within the array, averaged over all partitions. This returns a vector of the same length as the number of offspring in the array. family_size returns the posterior probability of observing one or more families of size 1, 2, ... , 15. It will be clear that we are unable to distinguish a single sibship with high probability from multiple families of the same size, each with low probability; this is the price we pay for being able to integrate out uncertainty in partition structure.
Step16:  Plotting this shows that we are roughly equally likely to observe a family of sizes one, two, three, four and five.
Step17:  Identifying fathers
In the section on paternityArray objects we used the method prob_array() to create an matrix of probabilities that each candidate is the father of each offspring, or that the father was missing from a paternityArray object. That was based only on comparing alleles shared between the mother, the candidate father and individual offspring.
We can get a similar matrix of probabilities of paternity after accounting for information about paternity that is shared between siblings by calling posterior_paternity_matrix on the sibshipCluster object sc. For example, father 'a_1' is the true sire of the furst five individuals in this toy dataset. We can compare his probabilities of paternity for his real offspring to the other unrelated progeny. We see that there is an increase in support for him as the father of the first five individuals after accounting for information shared between siblings, and a decrease for the other progeny.
Step18:  Since posterior_paternity_matrix returns a matrix, some wrangling is required to get the output into a format that is useful for a human to read. The following sections describe two helper functions that use posterior_paternity_matrix to create (1) a summary of probable mating events and (2) a summary of the paternity of each offspring.
Mating events
We very frequently want to know who the fathers of the offspring were to say something about mating events. There are several levels of complexity. Firstly, you can use the sires method to return a list of all the males who could possibly have mated with the mother. This is essentially identifying mating events, but doesn't say anything about the paternity of individual offspring. For many applications, that may be all you need because it's the mating events that are the unit of interest, not the number of offspring per se.
Once you have a sibshipCluster object, doing this is easy
Step19:  The columns in the output tell you several bits of information. The most interesting of these are
Step20:  The first five rows of the table above show that these fathers have posterior probabilities of paternity of one or close to one, and seem to have sired the correct numbers of offspring each. Of note is that although a_1 to a_4 have posterior probabilities of exactly one, the posterior probability for a_5 is slightly less than one. This is because the first four fathers sired multiple offspring, and there is shared information between siblings about the father, but this is not the case for father a_5.
After the true fathers there are a long list of extra candidates with very low posterior probabilities of paternity. In this case we know they are not true fathers, but in real life we would not, and we would like to account for this uncertainty. 
Paternity of individuals
If you are interested in the paternity of individual offspring we can look at the output of the paternity method of sibshipCluster objects. This returns a data frame with a row for each offspring, and the name and log posterior probability of paternity for the four most likely candidate fathers.
Step21:  In this case, the data are simulated, so we know who the real fathers are, and can print them with progeny.fathers. We see that the top candidate for all offspring is indeed the true father, with strong support (log probabilitites close to zero).
<ASSISTANT_TASK:>
Python Code:
import faps as fp
import numpy as np
import pandas as pd
print("Created using FAPS version {}.".format(fp.__version__))
Explanation: Sibship clustering
Tom Ellis, March 2017, updated June 2020
End of explanation
np.random.seed(867)
allele_freqs = np.random.uniform(0.3,0.5,50)
adults = fp.make_parents(100, allele_freqs, family_name='a')
Explanation: FAPS uses information in a paternityArray to generate plausible full-sibship configurations. This information is stored as a sibshipCluster object, and is the basis for almost all inference about biological processes in FAPS.
This notebook will examine how to:
Use a paternityArray to cluster offspring into plausible full-sibships.
Compare the relative support for different partitions structures
Make some basic inferences about the size and number of full sibships, and who is related to whom.
Note that this tutorial only deals with the case where you have a single maternal family. If you have multiple families, you can apply what is here to each one, but you'll have to iterate over those families. See the specific tutorial on that.
Generating a sibshipCluster object
We will begin by generating a population of 100 adults with 50 loci.
End of explanation
progeny = fp.make_sibships(adults, 0, [1,2,3], 5, 'x')
mothers = adults.subset(progeny.mothers)
patlik  = fp.paternity_array(progeny, mothers, adults, mu = 0.0015, missing_parents=0.01)
Explanation: We take the first individal as the mother and mate her to four males, to create three full sibships of five offspring. We then generate a paternityArray based on the genotype data.
End of explanation
sc = fp.sibship_clustering(patlik)
Explanation: It is straightforward to cluster offspring into full sibships. For now we'll stick with the default number of Monte Carlo draws.
End of explanation
%timeit fp.sibship_clustering(patlik, ndraws=100)
%timeit fp.sibship_clustering(patlik, ndraws=1000)
%timeit fp.sibship_clustering(patlik, ndraws=10000)
Explanation: The default number of Monte Carlo draws is 1000, which seems to work in most cases. I have dropped it to 100 in cases where I wanted to call sibship_clustering many times, such as in an MCMC loop, when finding every possible candidate wasn't a priority. You could also use more draws if you really wanted to be sure you had completely sampled the space of compatible candidate fathers. Speeds are unlikely to increase linearly with number of draws:
End of explanation
from scipy.cluster.hierarchy import dendrogram
import matplotlib.pyplot as plt
dendrogram(sc.linkage_matrix)
plt.show()
Explanation: We discussed this in figure 5 of the FAPS paper should you be interested in more on this.
Sibling relationships
Sibship clustering calculates likelihoods that each pair of offspring are full siblings, then builds a dendrogram from this. We can visualise this dendrogram if we so wish, although the output is not pretty.
End of explanation
sc.partitions
Explanation: Offspring individuals are labelled by their index in the array. Since full sibships are of size five we should expect to see clusters of {0,1,2,3,4}, {5,6,7,8,9} and {10,11,12,13,14}. This is indeed what we do see. What is difficult to see on the dendrogram are the branches connecting full siblings at the very bottom of the plot. If we bisect this dendrogram at different places on the y-axis we can infer different ways to partition the offspring into full siblings.
sc is an object of class sibshipCluster that contains various information about the array. Of primary interest are the set of partition structures inferred from the dendrogram. There are sixteen partitions - one for each individual in the array (i.e. one for each bifurcation in the dendrogram).
End of explanation
sc.partitions[2]
Explanation: What is key about partition structures is that each symbol represents a unique but arbitrary family identifier. For example in the third row we see the true partition structure, with individuals grouped into three groups of five individuals.
End of explanation
print(sc.lik_partitions) # log likelihood of each partition
print(np.exp(sc.prob_partitions)) # probabilities of each partition
Explanation: Beyond denoting who is in a family with whom, the labels are arbitrary, with no meaningful order. This partition would be identical to [0,0,0,0,0,1,1,1,1,1,2,2,2,2,2] or [10,10,10,10,10,7,7,7,7,7,22,22,22,22,22] for example.
Each partition is associated with a log likelihood and equivalent log probability. We can see from both scores that the third partition is most consistent with the data. This is of course the true partition.
End of explanation
np.exp(sc.prob_partitions).sum()
Explanation: We also see that the first and second partitions have non-zero, but small likelihoods. Parititons 5-8 have negative infinity log likelihood - they are incompatible with the data. These partitions split up true full siblings, and there is no way to reconcile this with the data. In real world situations such partitions might have non-zero likelihoods if they were an unrelated candidate male compatible with one or more offspring through chance alone.
In some cases there can be rounding error when log probabilities are exponentiated and probabilities do not sum to one. This is classic machine error, and the reason it is good to work with log values wherever possible. We can check:
End of explanation
sc.mlpartition
Explanation: You can directly call the most likely partition. This is somewhat against the spirit of fractional analyses though...
End of explanation
sibmat = sc.full_sib_matrix()
plt.pcolor(np.exp(sibmat))
plt.colorbar()
plt.show()
Explanation: For information about fine scale relationships, sc.full_sib_matrix() returns an $n*n$ matrix, where $n$ is the number of offspring. Each element describes the (log) probability that a pair of individuals are full siblings, averaged over partition structures and paternity configurations. If we plot this using a heatmap you can clearly see the five full sibships jump out as blocks of yellow (>90% probability of being full siblings) against a sea of purple (near zero probability of being full siblings).
End of explanation
# Lists indexing sires and dams
sires = [1]*5 + [2]*4 + [3]*3 + [4]*2 +[5]*1
dam   = [0] * len(sires)
np.random.seed(542)
allele_freqs = np.random.uniform(0.3,0.5,30)
adults  = fp.make_parents(1000, allele_freqs, family_name='a')
progeny = fp.make_offspring(adults, dam_list=dam, sire_list=sires)
mothers = adults.subset(progeny.mothers)
patlik  = fp.paternity_array(progeny, mothers, adults, mu= 0.0015, missing_parents=0.01)
sc = fp.sibship_clustering(patlik)
Explanation: Note that real datasets seldom look this tidy!
Inferring family structure
For this section we will simulate a slightly more interesting family structure. This block of code creates a half-sib array of 15 offspring from five fathers, where each father contributes five, four, three, two and one offspring respectively. It then performs sibship clustering on the array. We use 1000 candidate males and 50 loci.
End of explanation
sc.nfamilies()
Explanation: Number of families
We saw before that we could call a list of valid partitions for sc using sc.partitions. The output is not terribly enlightening on its own, however. We could instead ask how probable it is that there are x full sibships in the array, integrating over all partition structures. Here each number is the probability that there are 1, 2, ..., 15 families.
End of explanation
%matplotlib inline
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(111)
#ax.bar(np.arange(0.5, len(sc.nfamilies())+0.5), sc.nfamilies())
ax.bar(np.arange(1,16), sc.nfamilies())
ax.set_xlabel('Number of full sibships')
ax.set_ylabel('Posterior probability')
plt.show()
Explanation: We could show the same information graphically. Its clear that almost all the probability denisty is around $x=5$ families.
End of explanation
sc.family_size()
Explanation: Family size
We can also get the distribution of family sizes within the array, averaged over all partitions. This returns a vector of the same length as the number of offspring in the array. family_size returns the posterior probability of observing one or more families of size 1, 2, ... , 15. It will be clear that we are unable to distinguish a single sibship with high probability from multiple families of the same size, each with low probability; this is the price we pay for being able to integrate out uncertainty in partition structure.
End of explanation
fig = plt.figure()
ax = fig.add_subplot(111)
ax.bar(np.arange(len(sires))+0.5, sc.family_size())
plt.show()
Explanation: Plotting this shows that we are roughly equally likely to observe a family of sizes one, two, three, four and five.
End of explanation
np.exp(
    pd.DataFrame({
        "Before clustering" : patlik.prob_array()[:, 1],
        "After clustering"  : sc.posterior_paternity_matrix()[:, 1]
    })
)
Explanation: Identifying fathers
In the section on paternityArray objects we used the method prob_array() to create an matrix of probabilities that each candidate is the father of each offspring, or that the father was missing from a paternityArray object. That was based only on comparing alleles shared between the mother, the candidate father and individual offspring.
We can get a similar matrix of probabilities of paternity after accounting for information about paternity that is shared between siblings by calling posterior_paternity_matrix on the sibshipCluster object sc. For example, father 'a_1' is the true sire of the furst five individuals in this toy dataset. We can compare his probabilities of paternity for his real offspring to the other unrelated progeny. We see that there is an increase in support for him as the father of the first five individuals after accounting for information shared between siblings, and a decrease for the other progeny.
End of explanation
sc.sires()
Explanation: Since posterior_paternity_matrix returns a matrix, some wrangling is required to get the output into a format that is useful for a human to read. The following sections describe two helper functions that use posterior_paternity_matrix to create (1) a summary of probable mating events and (2) a summary of the paternity of each offspring.
Mating events
We very frequently want to know who the fathers of the offspring were to say something about mating events. There are several levels of complexity. Firstly, you can use the sires method to return a list of all the males who could possibly have mated with the mother. This is essentially identifying mating events, but doesn't say anything about the paternity of individual offspring. For many applications, that may be all you need because it's the mating events that are the unit of interest, not the number of offspring per se.
Once you have a sibshipCluster object, doing this is easy:
End of explanation
np.unique(patlik.fathers, return_counts=True)
Explanation: The columns in the output tell you several bits of information. The most interesting of these are:
label is the name of the candidate father
prob is the probability that the male sired at least one offspring with the mother, as a weighted average over partition structures. For example, 
offspring shows the expected number of offspring sired by the male, as a weighted average over partition structures. Specifically, it's the sum over rows from posterior_paternity_matrix; see below.
Note that if you have multiple maternal families the output will look a bit different. See the tutorial on multiple maternal families for more.
We can check this table makes sense by reviewing who the real fathers really are. This snippet gives a list of the names of the five true fathers, followed by the number of offspring sired by each.
End of explanation
sc.paternity()
Explanation: The first five rows of the table above show that these fathers have posterior probabilities of paternity of one or close to one, and seem to have sired the correct numbers of offspring each. Of note is that although a_1 to a_4 have posterior probabilities of exactly one, the posterior probability for a_5 is slightly less than one. This is because the first four fathers sired multiple offspring, and there is shared information between siblings about the father, but this is not the case for father a_5.
After the true fathers there are a long list of extra candidates with very low posterior probabilities of paternity. In this case we know they are not true fathers, but in real life we would not, and we would like to account for this uncertainty. 
Paternity of individuals
If you are interested in the paternity of individual offspring we can look at the output of the paternity method of sibshipCluster objects. This returns a data frame with a row for each offspring, and the name and log posterior probability of paternity for the four most likely candidate fathers.
End of explanation
progeny.fathers
Explanation: In this case, the data are simulated, so we know who the real fathers are, and can print them with progeny.fathers. We see that the top candidate for all offspring is indeed the true father, with strong support (log probabilitites close to zero).
End of explanation
<END_TASK>