<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Genotype data in FAPS
Step1:  Tom Ellis, March 2017
In most cases, researchers will have a sample of offspring, maternal and candidate paternal individuals typed at a set of markers. In this section we'll look in more detail at how FAPS deals with genotype data to build a matrix we can use for sibship inference.
This notebook will examine how to
Step2:  The object we just created contains information about the genotypes of each of the ten parent individuals. Genotypes are stored as NxLx2-dimensional arrays, where N is the number of individuals and L is the number of loci. We can view the genotype for the first parent like so (recall that Python starts counting from zero, not one)
Step3:  You could subset the array by indexes the genotypes, for example by taking only the first two individuals and the first five loci
Step4:  For realistic examples with many more loci, this obviously gets unwieldy pretty soon. It's cleaner to supply a list of individuals to keep or remove to the subset and drop functions. These return return a new genotypeArray for the individuals of interest.
Step5:  Information on individuals
A genotypeArray contains other useful information about the individuals
Step6:  make_sibships is a convenient way to generate a single half-sibling array from individuals in mypop. This code mates makes a half-sib array with individual 0 as the mothers, with individuals 1, 2 and 3 contributing male gametes. Each father has four offspring each.
Step7:  With this generation we can extract a little extra information from the genotypeArray than we could from the parents about their parents and family structure.
Step8:  Of course with real data we would not normally know the identity of the father or the number of families, but this is useful for checking accuracy in simulations. It can also be useful to look up the positions of the parents in another list of names. This code finds the indices of the mothers and fathers of the offspring in the names listed in mypop.
Step9:  Information on markers
Pull out marker names with marker. The names here are boring because they are simulated, but your data can have as exciting names as you'd like.
Step10:  Check whether the locus names for parents and offspring match. This is obvious vital for determining who shares alleles with whom, but easy to overlook! If they don't match, the most likely explanation is that you have imported genotype data and misspecified where the genotype data start (the genotype_col argument).
Step11:  FAPS uses population allele frequencies to calculate the likelihood that paternal alleles are drawn at random.
They are are useful to check the markers are doing what you think they are.
Pull out the population allele frequencies for each locus
Step12:  We can also check for missing data and heterozygosity for each marker and individual. By default, data for each marker are returned
Step13:  To get summaries for each individual
Step14:  In this instance there is no missing data, because data are simulated to be error-free. See the next section on an empircal example where this is not true.
Importing genotype data
You can import genotype data from a text or CSV (comma-separated text) file. Both can be easily exported from a spreadsheet program. Rows index individuals, and columns index each typed locus. More specifically
Step15:  Again, Python starts counting from zero rather than one, so the first column is really column zero, and so on. Because these are CSV, there was no need to specify that data are delimited by commas, but this is included for illustration.
Offspring are divided into 60 maternal families of different sizes. You can call the name of the mother of each offspring. You can also call the names of the fathers, with offspring.fathers, but since these are unknown this is not informative.
Step16:  Offspring names are a combination of maternal family and a unique ID for ecah offspring.
Step17:  You can call summaries of genotype data to help in data cleaning. For example, this code shows the proportion of loci with missing genotype data for the first ten offspring individuals.
Step18:  This snippet shows the proportion of missing data points and heterozygosity for the first ten loci. These can be helpful in identifying dubious loci.
Step19:  Multiple families
In real data set we generally work with multplie half-sibling arrays at once. For downstream analyses we need to split up the genotype data into families to reflect this. This is easy to do with split and a vector of labels to group offspring by. This returns a dictionary of genotypeArray objects labelled by maternal family. These snippet splits up the data and prints the maternal family names.
Step20:  Each entry is an individual genotypeArray. You can pull out individual families by indexing the dictionary by name. For example, here are the names of the offspring in family J1246
Step21:  To perform operations on each genotypeArray we now have to iterate over each element. A convenient way to do this is with dictionary comprehensions by separating out the labels from the genotypeArray objects using items.
As an example, here's how you call the number of offspring in each family. It splits up the dictionary into keys for each family, and calls size on each genotypeArray (labelled genArray in the comprehension).
Step22:  You can achieve the same thing with a list comprehension, but you lose information about family ID. It is also more difficult to pass a list on to downstream functions. This snippet shows the first ten items.
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import faps as fp
print("Created using FAPS version {}.".format(fp.__version__))
Explanation: Genotype data in FAPS
End of explanation
allele_freqs = np.random.uniform(0.3,0.5,10)
mypop = fp.make_parents(5, allele_freqs, family_name='my_population')
Explanation: Tom Ellis, March 2017
In most cases, researchers will have a sample of offspring, maternal and candidate paternal individuals typed at a set of markers. In this section we'll look in more detail at how FAPS deals with genotype data to build a matrix we can use for sibship inference.
This notebook will examine how to:
Generate simple genotypeArray objects and explore what information is contained in them.
Import external genotype data.
Work with genotype data from multiple half sib families.
Checking genotype data is an important step before committing to a full analysis. A case study of data checking and cleaning using an empirical dataset is given in here.
In the next section we'll see how to combine genotype information on offspring and a set of candidate parents to create an array of likelihoods of paternity for dyads of offspring and candidate fathers.
Note that the first half of this tutorial only deals with the case where you have a genotypeArray object for a single maternal family. If you have multiple families, you can apply what is here to each one, but at some point you'll have to iterate over those families. See below and the specific tutorial on that.
Currently, FAPS genotypeArray objects assume you are using biallelic, unlinked SNPs for a diploid. If your system deviates from these criteria in some way you can also skip this stage by creating your own array of paternity likelihoods using an appropriate likelihood function, and importing this directly as a paternityArrays. See the next section for more on paternityArray objects and how they should look.
genotypeArray objects
Basic genotype information
Genotype data are stored in a class of objects called a genotypeArray. We'll illustrate how these work with simulated data, since not all information is available for real-world data sets. We first generate a vector of population allele frequencies for 10 unlinked SNP markers, and use these to create a population of five adult individuals. This is obviously an unrealisticaly small dataset, but serves for illustration. The optional argument family_names allows you to name this generation.
End of explanation
mypop.geno[0]
Explanation: The object we just created contains information about the genotypes of each of the ten parent individuals. Genotypes are stored as NxLx2-dimensional arrays, where N is the number of individuals and L is the number of loci. We can view the genotype for the first parent like so (recall that Python starts counting from zero, not one):
End of explanation
mypop.geno[:2, :5]
Explanation: You could subset the array by indexes the genotypes, for example by taking only the first two individuals and the first five loci:
End of explanation
print(mypop.subset([0,2]).names)
print(mypop.drop([0,2]).names)
Explanation: For realistic examples with many more loci, this obviously gets unwieldy pretty soon. It's cleaner to supply a list of individuals to keep or remove to the subset and drop functions. These return return a new genotypeArray for the individuals of interest.
End of explanation
print(mypop.names) # individual names
print(mypop.size)  # number of individuals
print(mypop.nloci) # numbe of loci typed.
Explanation: Information on individuals
A genotypeArray contains other useful information about the individuals:
End of explanation
progeny = fp.make_sibships(mypop, 0, [1,2,3], 4, 'myprogeny')
Explanation: make_sibships is a convenient way to generate a single half-sibling array from individuals in mypop. This code mates makes a half-sib array with individual 0 as the mothers, with individuals 1, 2 and 3 contributing male gametes. Each father has four offspring each.
End of explanation
print(progeny.fathers)
print(progeny.mothers)
print(progeny.families)
print(progeny.nfamilies)
Explanation: With this generation we can extract a little extra information from the genotypeArray than we could from the parents about their parents and family structure.
End of explanation
print(progeny.parent_index('mother', mypop.names))
print(progeny.parent_index('father', mypop.names))
Explanation: Of course with real data we would not normally know the identity of the father or the number of families, but this is useful for checking accuracy in simulations. It can also be useful to look up the positions of the parents in another list of names. This code finds the indices of the mothers and fathers of the offspring in the names listed in mypop.
End of explanation
mypop.markers
Explanation: Information on markers
Pull out marker names with marker. The names here are boring because they are simulated, but your data can have as exciting names as you'd like.
End of explanation
mypop.markers == progeny.markers
Explanation: Check whether the locus names for parents and offspring match. This is obvious vital for determining who shares alleles with whom, but easy to overlook! If they don't match, the most likely explanation is that you have imported genotype data and misspecified where the genotype data start (the genotype_col argument).
End of explanation
mypop.allele_freqs()
Explanation: FAPS uses population allele frequencies to calculate the likelihood that paternal alleles are drawn at random.
They are are useful to check the markers are doing what you think they are.
Pull out the population allele frequencies for each locus:
End of explanation
print(mypop.missing_data())
print(mypop.heterozygosity())
Explanation: We can also check for missing data and heterozygosity for each marker and individual. By default, data for each marker are returned:
End of explanation
print(mypop.missing_data(by='individual'))
print(mypop.heterozygosity(by='individual'))
Explanation: To get summaries for each individual:
End of explanation
offspring = fp.read_genotypes(
    path = '../../data/offspring_2012_genotypes.csv',
    mothers_col=1,
    genotype_col=2)
Explanation: In this instance there is no missing data, because data are simulated to be error-free. See the next section on an empircal example where this is not true.
Importing genotype data
You can import genotype data from a text or CSV (comma-separated text) file. Both can be easily exported from a spreadsheet program. Rows index individuals, and columns index each typed locus. More specifically:
Offspring names should be given in the first column
If the data are offspring, names of the mothers are given in the second column.
If known for some reason, names of fathers can be given as well.
Genotype information should be given to the right of columns indicating individual or parental names, with locus names in the column headers.
SNP genotype data must be biallelic, that is they can only be homozygous for the first allele, heterozygous, or homozygous for the second allele. These should be given as 0, 1 and 2 respectively. If genotype data is missing this should be entered as NA.
The following code imports genotype information on real samples of offspring from half-sibling array of wild-pollinated snpadragon seedlings collected in the Spanish Pyrenees. The candidate parents are as many of the wild adult plants as we could find. You will find the data files on the IST Austria data repository (DOI:10.15479/AT:ISTA:95). Aside from the path to where the data file is stored, the two other arguments specify the column containing names of the mothers, and the first column containing genotype data of the offspring.
End of explanation
np.unique(offspring.mothers)
Explanation: Again, Python starts counting from zero rather than one, so the first column is really column zero, and so on. Because these are CSV, there was no need to specify that data are delimited by commas, but this is included for illustration.
Offspring are divided into 60 maternal families of different sizes. You can call the name of the mother of each offspring. You can also call the names of the fathers, with offspring.fathers, but since these are unknown this is not informative.
End of explanation
offspring.names
Explanation: Offspring names are a combination of maternal family and a unique ID for ecah offspring.
End of explanation
print(offspring.missing_data('individual')[:10])
Explanation: You can call summaries of genotype data to help in data cleaning. For example, this code shows the proportion of loci with missing genotype data for the first ten offspring individuals.
End of explanation
print(offspring.missing_data('marker')[:10])
print(offspring.heterozygosity()[:10])
Explanation: This snippet shows the proportion of missing data points and heterozygosity for the first ten loci. These can be helpful in identifying dubious loci.
End of explanation
offs_split = offspring.split(by = offspring.mothers)
offs_split.keys()
Explanation: Multiple families
In real data set we generally work with multplie half-sibling arrays at once. For downstream analyses we need to split up the genotype data into families to reflect this. This is easy to do with split and a vector of labels to group offspring by. This returns a dictionary of genotypeArray objects labelled by maternal family. These snippet splits up the data and prints the maternal family names.
End of explanation
offs_split["J1246"].names
Explanation: Each entry is an individual genotypeArray. You can pull out individual families by indexing the dictionary by name. For example, here are the names of the offspring in family J1246:
End of explanation
{family : genArray.size for family,genArray in offs_split.items()}
Explanation: To perform operations on each genotypeArray we now have to iterate over each element. A convenient way to do this is with dictionary comprehensions by separating out the labels from the genotypeArray objects using items.
As an example, here's how you call the number of offspring in each family. It splits up the dictionary into keys for each family, and calls size on each genotypeArray (labelled genArray in the comprehension).
End of explanation
[genArray.size for genArray in offs_split.values()][:10]
Explanation: You can achieve the same thing with a list comprehension, but you lose information about family ID. It is also more difficult to pass a list on to downstream functions. This snippet shows the first ten items.
End of explanation
<END_TASK>