Structure and group labels¶
Principle¶
Many statistics, including \(F_{ST}\), require that a structure is defined.
In EggLib, the structure is controlled by an instance of a specific
class (Structure
) that defines the groups and identifies the samples
that belong to them. These structure objects define which samples belong
to the outgroup and, how the other samples are organized in
individuals and/or populations and/or clusters of populations. In this
way, a three-level hierarchical structure can be defined. It is not required
to provide information for either of the three levels.
The rationale of this system is to allow processing the same data set
assuming different structurations (in order to either compare the effect of different
ways to group samples, or analyze different subsets of the data separately).
The approach is to keep the data set static, and provide a separate Structure
instance holding the description of the structure for each analysis.
There are threes ways of creating a Structure
instance. First, based on
group labels of an Align
instance. Second, by providing the sample
size of each population. Third, from a more flexible explicit description
of the structure as a dict
instance. Those approaches
are described in the next paragraphs.
Group labels¶
To account for structure, EggLib uses a specific system that consists in attaching labels to each samples. All samples bearing the same label are supposed to belong to the same group. There can be several labels per samples. We sometimes refer to an index among group labels as a level. Different levels of structure are aimed to represent either nested levels of structure (clusters of populations and/or populations and/or individuals) or alternative levels of structure. There are little restrictions on group labels besides being strings.
All Align
and Container
instances have a list of labels attached to each sample. They
can be set or edited using either LabelView
which is a part of
SampleView
(see Principle of proxy types) or direct class-level methods
get_label()
and set_label()
. By default, this list of labels is empty.
The aim of these group labels is essentially to be interpreted by the
class Structure
(actually through the method from_labels()
)
to be used for computing diversity statistics or restrict an operation
to a given subgroup.
In the following
we explain how to specify group labels in Fasta files such as EggLib will properly
interpret them and store them within an Align
or Container
instance.
Group labels in Fasta files¶
Single group label¶
Let align2.fas
be a Fasta file with six samples, the first three belonging to
population “pop1” and the other three belonging to population “pop2”. EggLib
supports a specific system of tags within sequence headers in the Fasta format
to indicate group labels. The tags must appear as suffix starting with a @
followed by a string, as in the following example:
>sample1@pop1
ACCGTGGAGAGCGCGTTGCA
>sample2@pop1
ACCGTGGAGAGCGCGTTGCA
>sample3@pop1
ACCGTGGAGAGCGCGTTGCA
>sample4@pop2
ACCGTGGAGAGCGCGTTGCA
>sample5@pop2
ACCGTGGAGAGCGCGTTGCA
>sample6@pop2
ACCGTGGAGAGCGCGTTGCA
To import group labels, one is required to set the labels option of
from_fasta()
to True
:
aln2 = egglib.io.from_fasta('align2.fas', alphabet=egglib.alphabets.DNA, labels=True)
print(aln2.get_name(0), aln2.get_label(0, 0))
sample1 1
If labels have not been imported like in the following:
aln2 = egglib.io.from_fasta('align2.fas', alphabet=egglib.alphabets.DNA)
print(aln2.get_label(0, 0))
In that case, print(aln2.get_label(0, 0))
would cause an error because,
by default, there is no group labels at all included in Align
instances.
Note
io.from_fasta()
has an option (label_marker) to use a different character than
@
to separate the name and the tag.
Multiple group labels¶
There can be any number of group levels, either nested or not. To specify several labels for a sample, one can write several strings separated by commas, as in the following example:
>sam1@c1,i1
ACCGTGGAGAGCGCGTTGCA
>sam2@c1,i1
ACCGTGGAGAGCGCGTTGCA
>sam3@c1,i2
ACCGTGGAGAGCGCGTTGCA
>sam1@c1,pop1,i1
ACCGTGGAGAGCGCGTTGCA
>sam5@c2,pop1,i1
ACCGTGGAGAGCGCGTTGCA
>sam6@c2,pop1,i2
ACCGTGGAGAGCGCGTTGCA
The example above also demonstrate that it is possible to omit group labels for part of the samples,
although it is probably better to avoid it (because it is error-prone). Labels
absent for a given level are not added or initialised in any way. As a result, if the
file shown above is saved as align3.fas
, the code:
aln3 = egglib.io.from_fasta('align3.fas', alphabet = egglib.alphabets.DNA, labels=True)
print(aln3.get_name(0))
print(aln3.get_label(0, 0))
print(aln3.get_label(0, 1))
will yield:
sample1
c1
i1
but the code:
print(aln3.get_label(0, 2))
will result in an error.
So, at this point, one should understand labels as list of 0, 1 or more
arbitrary identifiers attached to each sample. How this labels will be
used to group samples in individuals, populations or possibly multi-level
hierachical structure is up to the Structure
class.
Note
The separator can also be changed. This can be done using the
label_separator of the io.from_fasta()
method.
Outgroup specification¶
Tools analysing diversity in EggLib can account for one or more outgroup samples. If individuals are defined in the main group (ingroup), it is required that outgroup samples also come as one or more individuals sharing the same ploidy.
Individuals must be specified by the label #
(obligatory as the only
or first label), followed by maximum one other label (if the individual
level is to be considered). This, outgroup samples are denoted by the
tag @#
or @#,IND
where IND
is an individual label. Note that
if there are individuals in the ingroup then there must be also individuals
in the outgroup (with the same ploidy).
Very importantly, one must keep in mind that neither io.from_fasta()
nor Align
have a notion of the outgroup. They don’t interpret
the #
label as special and don’t process outgroup samples differently
of other samples. It will be the job of Structure
to separate the
outgroup from the rest of the samples. This means that, if you have outgroup
samples including in your data, you must use a Structure
instance.
Also, if you want to use another label than #
to identify outgroup samples,
you need to tell it to Structure
(or to the method struct_from_labels()
).
Creating a structure from an alignment¶
To present the usage of the Structure
class, we will use a
complete, albeit over-simplified example. Consider the Fasta
alignment below:
>sample01@c1,p1,i01
CTTCCGGGAAGCGCCAGCAGAAGGTTGCTGCTAAGGCCCGCACACGTCTGCAGCACTTCG
>sample02@c1,p1,i01
CTTCCGCGCAGGGCCAGGAGCATGTAGCTTCTAAGGCTTGCACAGGTCTTCAGCACTACG
>sample03@c1,p1,i02
CTTCCGCGCAGGGCCAGGAGCATGTAGCTTCTAAGGCTTGCACAGGTCTTCAGCACTACG
>sample04@c1,p1,i02
CTTCCGCGCAGGGCCAGGAGCATGTAGCTTCTAAGGCTTGCACAGGTCTTCAGCACTACG
>sample05@c1,p1,i03
CTACCGTGACGAGCTAGCCGAGCCTGACGCAGGGGGCGAGTAAGGGAGATTACGACTTGG
>sample06@c1,p1,i03
CTTCCGCGCAGGGCCAGGAGAATGTTGCTTCTAAGGCTTGCACAGGTCTTCAGCACTAAG
>sample07@c1,p2,i04
CTGCTATGACGAACTACCCGAGCCTGGGGCATGGGGCGTGTATGGGAGCTTACAACTTGG
>sample08@c1,p2,i04
CTTACGCGACGTGCCAGCATAGGGAAGCTGCTAAGGCCTGCACACGTCCGCAGCACTACG
>sample09@c1,p2,i05
CTGCTATGACGAACTACCCGAGGCTGGGGCATGGGGCGTGTATGGGAGCTTACAACTTGG
>sample10@c1,p2,i05
CTGCTATGACGAACTACCCGAGGCTGGGGCATGGGGCGTGTATGGGAGCTTACAACTTGG
>sample11@c1,p2,i06
CTTACGCGACGCGCCAGCAGAGGGATGCTGCTAAGGCCTGCACACGTCCGCAGCACTACG
>sample12@c1,p2,i06
CTTACGCGACGCGCCAGCAGAGGGATGCTGCTAAGGCCTGCACACGTCCGCAGCACTACG
>sample13@c1,p2,i07
CTGCTATGACGAACTACCCGAGCCTGGGGCATGGGGCGTGTATGGGAGCTTACAACTTGG
>sample14@c1,p2,i07
CTGCTATGACGAACTACCCGAGGCTGGGGCATGGGGCGTGTATGGGAGCTTACAACTTGG
>sample15@c2,p3,i08
CTCCGGGGCCGGTTTCGCATAACGTCGCGCAGGGGACGTGTAGGGGCGCATACACCTGGG
>sample16@c2,p3,i08
CTCCGGGGCCGGTTTCGCATAACGTCGCGCAGGGGACGTGTAGGGGCGCATACACCTGGG
>sample17@c2,p3,i09
TTGCCGGGTCGAACTAGCCGACCTTGGCGCAGGGGTCGTTTAAGGGTCCTTACAACTTGG
>sample18@c2,p3,i09
TTGCCGGGTCGAACTAGCCGACCTTGGCGCAGGGGTCGTTTAAGGGTCCTTACAACTTGG
>sample19@c2,p3,i10
CTCCGGCGCCGGTTTCGCATAACGTCGCGCAGGGGACGTGTAGGGGCGCATACACCTGGG
>sample20@c2,p3,i10
TTGCCGGGTCGAACTAGCCGACCTTGGCGAAGGGGTCGTTTAAGGGACCTTACAACTTGG
>sample21@c2,p4,i11
TTGCCAGGACGAACTAGCCGCGCCTGGCGCAGGGGTCGTTTAAGGGAGCTTACAACTTGG
>sample22@c2,p4,i11
TTGCCAGGACGAACTAGCCGCGCCTGGCGCAGGGGTCGTTTAAGGGAGCTTACAACTTGG
>sample23@c2,p4,i12
TTGCCGGGACGAACTAGCCGAGCCTGGCGCAGGGGTCGTTTAAGGGAGCTTACAACTTGG
>sample24@c2,p4,i12
TTGCCGGGACGAACTAGCCGAGCCTGGCGCAGGGGTCGTTTAAGGGAGCTTACAACTTGG
>sample25@c2,p5,i13
CTACCGTGACGAACTAGCCGAGCCTGGCGCAGGGGGCGAGTAAGGGAGAGTACAACTTGG
>sample26@c2,p5,i13
CTACCGTGACGAACTAGCCGAGCCTGGCGCAGGGGGCGAGTAAGGGAGAGTACAACTTGG
>sample27@c2,p5,i14
CTACCGTGACGAACTAGCCGAGCCTGGCGCAGGGGGCGAGTAAGGGAGAGTACAACTTGG
>sample28@c2,p5,i14
CTACCGTGACGAACTAGCCGAGCCTGGCGCAGGGGGCGAGTAAGGGAGAGTACAACTTGG
>sample29@c2,p5,i15
TTGCCGCGACGAACTAGCCGAGCCTGGCGCAGGGGTCGTTTAAGGGAGCTAACAACTTGG
>sample30@c2,p5,i15
CTACCGTGACGAACTAGCCGAGCCTGGCGCAGGGGGCGAGTAAGGGAGAGTACAACTTGG
>sample31@#,i98
CATACCACCTTGGCCCGGAGAGTGCGGAGTACCGGGCGTGGAAGGCTGCATGCAAATGGA
>sample32@#,i98
CATACCACCTTGGCCCGGAGAGTGCGGAGTACCGGGCGTGGAAGGCTGCATGCAAATGGA
>sample33@#,i99
CATACCACCTTGGCCCGGAGAGAGCGCAGTGCCGGGCGTGGAAGGCTGCATTCAAATGCG
>sample34@#,i99
CATACCACCTTGGCCCGGAGAGAGCGCAGTGCCGGGCGTGGAAGGCTGCATTCAAATGCG
It has a total of 30 ingroup samples and 4 outgroup samples. These are
actually respectively 15 and 2 individuals, and the ingroup is organized
in two clusters of respectively two and three populations, themselves composed of two,
three, or four individuals each. Remember the labels are arbitrary. In this case,
cluster labels are c1
and c2
, population labels p1
to p5
and individual labels i01
to i15
(i98
and i99
for the
two outgroup individuals).
Let use name this file align5.fas
and import it with group labels:
>>> aln = egglib.io.from_fasta('align5.fas', alphabet=egglib.alphabets.DNA, labels=True)
Now, we will directly show a Structure
instance incorporating
all structure information (all three levels) can be created:
>>> struct = egglib.struct_from_labels(aln, lvl_clust=0, lvl_pop=1, lvl_indiv=2)
>>> print(struct.as_dict())
({'c2': {'p3': {'i09': [16, 17], 'i08': [14, 15], 'i10': [18, 19]},
'p5': {'i13': [24, 25], 'i15': [28, 29], 'i14': [26, 27]},
'p4': {'i11': [20, 21], 'i12': [22, 23]}},
'c1': {'p1': {'i01': [0, 1], 'i03': [4, 5], 'i02': [2, 3]},
'p2': {'i05': [8, 9], 'i04': [6, 7], 'i07': [12, 13], 'i06': [10, 11]}}},
{'i99': [32, 33], 'i98': [30, 31]})
We used the function struct_from_labels()
that generates a new
Structure
instance based on group labels of an Align
(or Container
). To use this method, it is necessary to tell
which group level corresponds to the clusters, populations, and individuals
in such a way that they are properly hierarchical. It is possible to skip
any of these three levels of structure, simply by dropping the corresponding
option parameter(s).
The method as_dict()
is aimed to provide an intuitive
representation of the structure held by the instance. In practice, it is as
intuitive as possible while being flexible enough to represent all possible cases.
Dictionary representation of Structure
instances¶
It is a tuple
containing two items,
each being a dict
. The first one represents the ingroup and the second
represents the outgroup.
The ingroup dictionary is itself a dictionary holding more dictionaries, one for each cluster of populations. Each cluster dictionary is a dictionary of populations, populations being themselves represented by a dictionary. A population dictionary is, again, a dictionary of individuals. Finally individuals are represented by lists or integers.
An individual list contains the index
of all samples belonging to this individual. For haploid data, individuals
will be one-item lists. In other cases, all individual lists are required to have
the same number of items (consistent ploidy). Note that, if the ploidy is more
than one, nothing enforces that samples of a given individual are grouped within
the original data, meaning that you can shuffle labels in Align
instances
(or in your Fasta file) if you need to.
The keys of the ingroup dictionary are the labels identifying each cluster. Within a cluster dictionary, the keys are population labels. Finally, within a population dictionary, the keys are individual labels.
The second dictionary represents the outgroup. Its structure is simpler: it has individual labels as keys, and lists of corresponding sample indexes as values. The outgroup dictionary is similar to any ingroup population dictionary. The ploidy is required to match over all ingroup and outgroup individuals.
If we go back to our example, we see that the returned dictionary for the
ingroup has two items, with keys c1
and c2
, respectively, and that
the correct structure appears at lower levels, with two-item (diploid) individuals
within populations withing clusters. Similarly, the two outgroup individuals,
labelled i98
and i99
, appear as expected in the second dictionary returned by
the as_dict()
method. Ultimately, the values contained by the
lists are the lowest levels are the index referring to the original
Align
instance (from 0 to 29 in the ingroup, 30 to 33 in the
outgroup).
Alternative structure¶
Occasionally, one will want to generate different Structure
instances
based on different levels of structure in group labels (for example if there are
alternative structurations of the data). It is not required that all levels of
a Structure
instances are populated, and it is not necessary to import
all structure levels of an Align
. The example below demonstrates all
this by importing the first level (previously, clusters) as populations in a new
instance, skipping all other information:
>>> struct2 = egglib.struct_from_labels(aln, lvl_pop=0)
>>> print(struct2.as_dict())
({None: {'c2': {'24': [24], '25': [25], '23': [23], '27': [27],
'15': [15],'14': [14], '17': [17], '16': [16],
'19': [19], '18': [18], '22': [22], '28': [28],
'26': [26], '29': [29], '20': [20], '21': [21]},
'c1': {'11': [11], '10': [10], '13': [13], '12': [12],
'1': [1], '0': [0], '3': [3], '2': [2], '5': [5],
'4': [4], '7': [7], '6': [6], '9': [9], '8': [8]}}},
{'33': [33], '32': [32], '31': [31], '30': [30]})
Note that it is also possible to recycle an already existing Structure
instance instead creating a new one (with the method from_labels()
of Structure
instances.).
Since we did not specify any group label index for the cluster level, there is
no information regarding clusters, and all populations are placed in a single
cluster. The default label is None
in that case. The two labels c1
and c2
are
now considered as population labels. At the lowest level (also in the outgroup),
all samples are placed in a single-item individuals because, likewise, no index
has been provided for the individual level. Then, haploidy is assumed, and the
sample index is used as default value for individual labels (incremented in the
outgroup).
This example demonstrates that the group labels in Align
instances have no particular meaning per se until they are interpreted
while configuring a Structure
instance.
Creating a structure manually¶
Simple structure¶
If your data are organized in an intuitive way (that is, samples organized
per individual and individuals grouped per population), and if the cluster
level is not needed, you can use the function struct_from_samplesizes()
.
This function takes a list of sample sizes (one item per population).
For example, if your dataset contains two populations of 20 haploid
individuals, you can enter simply:
>>> struct = egglib.struct_from_samplesizes([20, 20])
>>> print(struct.as_dict())
({None: {'pop1': {'idv1': [0], 'idv2': [1], 'idv3': [2], 'idv4': [3],
'idv5': [4], 'idv6': [5], 'idv7': [6], 'idv8': [7],
'idv9': [8], 'idv10': [9], 'idv11': [10],'idv12': [11],
'idv13': [12], 'idv14': [13], 'idv15': [14], 'idv16': [15],
'idv17': [16], 'idv18': [17], 'idv19': [18], 'idv20': [19]},
'pop2': {'idv21': [20], 'idv22': [21], 'idv23': [22], 'idv24': [23],
'idv25': [24], 'idv26': [25], 'idv27': [26], 'idv28': [27],
'idv29': [28], 'idv30': [29], 'idv31': [30], 'idv32': [31],
'idv33': [32], 'idv34': [33], 'idv35': [34], 'idv36': [35],
'idv37': [36], 'idv38': [37], 'idv39': [38], 'idv40': [39]}}}, {})
This function supports ploidy and outgroup individuals, so you can also declare, for example, two populations of 10 diploid individuals plus one outgroup individual:
>>> struct = egglib.struct_from_samplesizes([10, 10], ploidy=2, outgroup=1)
>>> print(struct.as_dict())
({None: {'pop1': {'idv1': [0, 1], 'idv2': [2, 3], 'idv3': [4, 5],
'idv4': [6, 7], 'idv5': [8, 9], 'idv6': [10, 11],
'idv7': [12, 13], 'idv8': [14, 15], 'idv9': [16, 17],
'idv10': [18, 19]},
'pop2': {'idv11': [20, 21], 'idv12': [22, 23], 'idv13': [24, 25],
'idv14': [26, 27], 'idv15': [28, 29], 'idv16': [30, 31],
'idv17': [32, 33], 'idv18': [34, 35], 'idv19': [36, 37],
'idv20': [38, 39]}}}, {'idv21': [40, 41]})
Be careful that the order of samples in the Align
or
Site
you’ll be analyzing with the resulting Structure
instance must be consistent. Populations must be grouped together in the
order indicated (if sample sizes differ), as well as individuals in
populations, and the outgroup must be at the end.
Fully flexible dictionary¶
It is possible to create a Structure
instance from user-provided
data formatted as dictionaries, using either the function struct_from_dict()
or the equivalent method
of Structure
instances
to recycle an existing instance.
This approach allows maximal flexibility but requires that you create a
properly formatted dictionary. Both methods take an ingroup and an outgroup
arguments, which are formatted exactly as the output of as_dict()
(see Dictionary representation of Structure instances). This feature can be used to import structure
information without relying on EggLib’s tag system of Fasta sequence names.
Using the structure¶
Once a Structure
has been configured to represent the structuration
of the data set, it can be used as a descriptor while computing diversity
statistics. This will make available a wide array of statistics requiring
this type of information. For example, the statistics with codes
Fis
, Snn
, WCist
, and WCisct
require individual and/or
population structure information and won’t be computed if no structure
is provided:
>>> cs = egglib.stats.ComputeStats()
>>> cs.add_stats('Fis', 'FistWC', 'FisctWC', 'Snn')
>>> print(cs.process_align(aln))
{'Snn': None, 'Fis': None, 'WCisct': None, 'WCist': None}
To provide the Structure
to stats.ComputeStats
, one just
needs to pass the instance as a value for the struct argument of the
class constructor (or the configure()
method of ComputeStats
instances):
>>> cs.configure(struct=struct)
>>> print(cs.process_align(aln))
{'Snn': 0.5111111111111112, 'Fis': 0.6442024052643518,
'FistWC': (0.39885944313988586, 0.4485871173702674, 0.6685233526761218),
'FisctWC': (0.39885944313988575, 0.4964142771064885, 0.2339234438730581, 0.6972741981129913)}
The code above shows that, with proper structure, we can compute statistics
taking into account the individual, population, and cluster levels. In particular,
FisctWC
takes all levels into account. In comparison, FistWC
ignores
the cluster level, but nothing prevents you from computing it at this point.
The code below shows that we can analyse the same data with a different structure
(using the second instance we created before, using the clusters as populations
and ignoring other levels):
>>> cs.configure(struct=struct2)
>>> print(cs.process_align(aln))
{'Snn': 0.9666666666666667, 'Fis': None, 'FistWC': None, 'FisctWC': None}
Since the individual level is not available, the statistics Fis
,
WCist
, and WCisct
(which also requires the cluster level) cannot
be computed. Only Snn
can. It is still possible to call for statistics
that cannot be computed, but their value will be set to None
.