import numpy as np
from IPython.display import HTML
from bokeh.plotting import output_notebook, show
import genomes_dnj.lct_interval.series_plots as dm
import genomes_dnj.lct_interval.anal_series as an
import genomes_dnj.lct_interval_snp_anal.lct_interval_snp_anal as snp_anal
import genomes_dnj.lct_interval.series_masks as sm
from genomes_dnj.lct_interval.bokeh_match_count_plot import snp_match_count_plot_cls as smcp
output_notebook(hide_banner=True)
Full details of all of the methods are available in the source code. This notebook describes them and provides some examples to illustrate their use.
SNPs were grouped into series through an algorithmic process. The data for that process was a two dimensional array. There was a row for each SNP expressed by at least 16 1000 genomes chromosome samples. There were 5008 columns, one for each 1000 genomes chromosome sample. Values were true if a sample expressed a SNP and false if it did not.
The process for identifying a series started with a specific SNP and an array of the column indexes for the samples that expressed that SNP. Those indexes were used to generate an array of counts of the number of the test samples that expressed each SNP farther along the chromosome. The counts were compared with the 90% match criterion to identify additional SNPs that were expressed by 90% of the test samples. A second array was generated with counts of all samples that expressed each of the SNPs selected by the initial match process. The second match check required the count of the test samples expressing the SNP to be at least 90% of the count of all samples expressing the SNP.
This process was carried out recursively by making the last identified SNP in a series the start of another attempt to find more SNPs in the series. The process was continued until a 1,000,000 base interval was scanned without finding any additional SNPs in the series. The cases where this process resulted in at least 4 SNPs were identified as a series. Samples that expressed at least 90% of the SNPs in a series were considered to express the series.
Rounding can result in inclusion samples that express less than 90% of the SNPs. The recursive process can favor identification of series that have some history of fragment generation of the lower part of the series. But, samples that express the fragments are not counted as expressing the series.
Plots were constructed by selecting chromosome samples that satisfied some filter criteria. The criteria are based on a yes list of series and a no list of series. To be selected a chromosome sample must express all of the series in the yes list and none of the series in the no list.
Series are identified by the number of SNPs in the series and the number of 1000 genomes chromosome samples that express it. The series 11_765 associated with lactase persistence includes 11 SNPs and is expressed by 765 of the 5008 1000 genomes samples of chromsome 2.
import genomes_dnj.lct_interval.plot_colors as clr
HTML(clr.pop_colors)
Series are represented in a plot by a colored rectangle. The hue of the rectangle represents the most overexpressed population for the chromosome samples that express the series. More saturated colors represent more overexpressed populations. The table above shows the fully saturated colors used to represent overexpression of series by different 1000 genomes populations. The plot x axis maps to positions on chromosome 2. The height of the rectangle represents the log of the number of 1000 genomes chromosome 2 samples that express the series. Each SNP in the series is represented with a vertical line. Height and color are properties of the series regardless of how many chromosome samples express the series in a particular plot.
Most plots in these notebooks are superset plots. The intent of a superset plot is to find all of the series in the studied region expressed by some set of chromosome samples. Some of those series may be expressed specifically by the samples that satisfy the filter test. For other series those samples are just a subset of the 1000 genomes chromosome samples that express the series. A yes no comparison is used to select input chromosome samples. The output series are identified by trying to match the input samples with the samples that express each of the series in the studied region. Those series that meet the match criterion are added to the result set. When the intent was to limit the series to those expressed by all of the selected chromosome samples, the match criterion was generally 90%. When the intent was to identify all of the series most commonly expressed by the selected chromosome samples, the match criterion was generally 50%. The use of a 50% match prevents the selection of overlapping series that are expressed by different subsets of the selected chromosome samples. When the intent was to display a complete hierarchy of series that are all expressed by all the chosen set of chromosome samples, the match criterion was generally 0.1%. The data associated with each plot shows the number of input samples that actually express each of the result series.
lct_plt_obj = dm.superset_yes_no([dm.di_11_765], min_match=0.9)
plt = lct_plt_obj.do_plot()
show(plt)
This plot is an example of a superset plot. The only criterion for selecting chromosome samples is the expression of the series 11_765. The series in the plot are ones that are expressed by at least 90% of the 1000 genomes phase 3 chromsome 2 samples that express the series 11_765.
HTML(lct_plt_obj.get_html())
This table shows the plot data for the 11_765 superset plot. Each matched series has a row in the table. Series are ordered by the number of 1000 genomes chromosome samples that express the series. Matches give the number of those samples selected for the plot that express the particular series.
For example the top line indicates that the 10 SNP series expressed by 2206 chromosome 2 samples is expressed by all 765 of the chromsomes samples that satisfied this plot's yes no criteria. Those samples were 35% of the total 2206 1000 genomes chromosome 2 samples that express the series 10_2206. Those samples were 100% of the 765 chromosome samples that satisfied this plot's yes no criteria.
Number of matches and relative population expression are broken out by population. Both the African and South Asian populations are divided between those that came from the native area and those who live externally. The data in this table was the reason for that choice. The 32 instances of 11_765 found in American Southwest and Caribbean African populations are very unlikely to have come from Africa and significantly over represent the frequency of lactase persistence in Afican populations.
For calculations of population overexpression chromosome samples that failed the no test were excluded from the population under consideration. The relative expression calculation is based on a ratio of the samples that matched to the population size minus the samples that failed the no test.
A subset plot starts with a test sample set and identifies series that are expressed by some subset of that sample set. The samples that express each of the series in the result set must be a subset of the samples in the test set up to some match criterion. Generally, a 90% match criterion was used for these plots.
plt_obj = dm.subset_yes_no([dm.di_193_843], min_match=0.8)
plt = plt_obj.do_plot()
show(plt)
The plot above shows series that are expressed by chromosome 2 samples that are subsets of the chromosome samples that express the series 193_843. Those series are all part of the hierarchy that has been named the SAS tree. The match criterion has been relaxed to 0.8 because some of the chromosome samples that express the series 8_267 have experienced a recombination event that cause loss of enough 193_843 SNPs for those samples not to match 193_843.
HTML(plt_obj.get_html())
The data in the table above shows 19 series where 100% of the chromsome samples that express the series also express 193_843.
Three other series, 4_815, 6_713, and 8_267, are of particular interest. All three are contained within the interval of chromosome 2 covered by 193_843. Recombination events that lose association with 193_843 must be somewhere within the DNA covered by the 193_843 series. More than 80% of the chromosome samples that express any of the three series also express 193_843. But a signficant number of those chromosomes samples do not.
The reasons are two recombination events that are analyszed in the part of this notebook focused on allele masks. One involved a chromosome that expressed the series 8_267 and 6_713. It resulted from a recombination event that lost the last 23 SNPs of 193_843. The result of that recombination event became overexpressed through a selection process associated with the appearance of the series 9_39.
The other recombination event involved a chromosome that expressed 193_843 and the other major series that form the root of the SAS tree with a chromosome that expressed the series 5_684 that is a major component of the EAS tree hierarchy. That recombination event resulted in a number of overexpressed variations on the association of 4_815 and 5_684 that appear to have conferred some kind of selecrtive advantage on the recombinant chromosomes.
Another class of cases include 15_59, 20_56, 5_28, and 5_25. All of these series cover parts of the studied region of chromosome 2 that are beyond the end of 193_843. As is commonly the case hierarchies selected by these series have experienced recombination events between the more stable lower part of the region and the less stable upper part.
The two other cases are the series 10_25 and 8_17. Both have experienced recombination events within the region covered by 193_843. But neither case has generated an overexpressed result.
A basis plot is used to match samples to the most specific series in an hierarchy that the samples express. Basis filtering can be used in combination with yes no comparisons to explore recombination events for samples that express particular parts of an hierarchy.
plt_obj = dm.basis_yes_no([dm.di_6_1503, dm.di_26_1414])
plt = plt_obj.do_basis_plot()
show(plt)
This plot shows the basis series for all of the chromosome 2 samples in the phase 3 1000 genomes data that express both the series 6_1503 and 26_1414. These criteria identify the chromosome 2 samples that express hierarchies that are part of the EUR tree.
Series are labeled in the top left corner by their id. The number at the end of the series is the number of chromosome 2 samples that do express it and do not express any other more specific series in the lactase persistence region of chromosome 2. For example the 135 chromosome 2 samples that express the series 4_911 as a basis series do not express the series 11_765. All but 3 of the chromosome samples that express 11_765 also express 4_911. But, since 911 chromosamples express 4_911 and only 765 chromosome samples express 11_765, 11_765 is the more specific series.
The color used to represent the basis series is the color used to represent the whole series. But the height of the basis series rectangles is a function of the log of the number of chromosome samples that actually express it as a basis series.
HTML(plt_obj.get_basis_html())
The data accompanying the basis plots needs to be consulted to understand the population distribution of the chromosome samples that express a particular series as a basis. The data in this table for the series 26_1414, 4_911, and 11_765 provide some insight in the evolution of lactase persistence.
The 202 chromosome samples that express 26_1414 as a basis suggest that the root series of the EUR tree already had some kind of selective advantage. Those samples are heavily under expressed in African populations, modestly under expressed in European populations, modestly overexpressed in American and South Asian populations, and most heavily overexpressed in East Asian populations. The root series of the EUR tree does seem to have been favored by the expansion out of Africa. But, it appears to have been distributed over all of the other 1000 genomes population regions.
The 135 chromosomes that express 4_911 as a basis do not include any African or East Asian chromosomes. They are modestly overexpressed for European populations, more overexpressed for South Asian populations and most overexpressed for American populations. It seems unlikely that the American samples got to America through the migration either of East Asian or European populations.
Another kind of plot is used to examine the number of SNPs from a series that are expressed by samples from the 1000 genomes data. This plot shows number of series SNPs on the x axis and number of 1000 genomes chromosome samples that express that number of the series SNPs on the y axis.
This example shows that kind of plot for the 117 SNPs from the series 117_1685
count_data = an.sa_117_1685.unique_snps_per_allele()
plt_obj = smcp(count_data)
plt = plt_obj.do_plot()
show(plt)
This plot shows the numbers of samples that express different numbers of SNPs from the series 117_1685
count_data
This data shows the SNP fragment size and the number of samples expressing that number of 117_1685 SNPs as an array of tuples. Each tuple associate a number of SNPs with the number of samples that express that number of SNPs
In this example the distribution is strongly peaked with 2680 samples expressing none of the 117 SNPs and 1182 samples expressing all 117 of them. This kind of result was observed for all of the 76 series documented in the detailed results of this study.
The distribution includes several cases of peaks for significant numbers of samples expressing the same size fragment of the 117 SNPs. Investigations documented in the detailed results showed that multiple fragments of the same size generally indicate expression of multiple instances of fragments with the same or similar sequences of SNPs.
sa_64_1575 = snp_anal.series_anal_cls(dm.di_64_1575, sm.series_data)
count_data = sa_64_1575.unique_snps_per_allele()
plt_obj = smcp(count_data)
plt = plt_obj.do_plot()
show(plt)
This plot shows the SNP distribution analysis for the series 64_1575. It is another case where substantial numbers of fragments exist that generally can be traced back to recombination events and some kind of process that selected the fragment for overexpression.
count_data
This data is the array of SNP counts and sample numbers for the 64_1575 analysis.
Detailed analysis of series fragments can be carried out with allele mask plots. An allele mask is an array with a boolean value for each of the 5008 chromosome samples in the 1000 genomes phase 3 data. The value is true for those chromosome samples that satisfy some logical test. Allele masks for samples that express fragments of a specific size from some series are useful for analyzing those series fragments.
from genomes_dnj.lct_interval.series_masks import am_193_843, nam_193_843
an.sa_193_843.unique_snps_per_allele(am_193_843)
The an.sa_193_43 object is an instance of the series_anal_cls. It is used for analysis of the samples that express any subset of the 193 SNPs in the series 193_843. The allele mask am_193_843 is a boolean mask that identifies the 1000 genomes chromosome 2 samples that express the series 193_843. The unique_snps_allele method returns a list of the number of SNPs and the number of chromosome samples that express that number of 193_843 SNPs. This data shows that 696 of the 843 chromosome samples that express 193_843 express all 193 of its SNPs. Another 120 of the chromosomes express 192 of them. Smaller numbers of chromosome samples do express smaller numbers of SNPs that are fragments of the whole series.
an.sa_193_843.unique_snps_per_allele(nam_193_843)
The nam_193_843 boolean mask identifies chromosome 2 samples that do not express 193_843. The results are a list of numbers of 193_843 SNPs and the number of chromsome 2 samples that express that number of 193_843 SNPs. The results show that 3458 samples express none of the 193 SNPs.
The data does show that a significant number of samples express 1 or 2 of the 193 SNPs. It also shows that 5, 6, 12, 28, 29, 55 and 170 SNP fragments are expressed by numbers of samples that indicate substantial overexpression. There also is a background of other fragments of the 193 SNPs that could be available for some kind of future selection process.
aps_29, am_29 = an.sa_193_843.snps_from_aps_value(29, nam_193_843)
aps_29
This data looks at the 29 SNP fragments that are expressed by 94 chromosome samples. The method snps_from_aps_value takes a count of the number of SNPs and a boolean mask of chromosome 2 samples. It returns an array of the number of samples that express each of the 193 SNPs and a boolean mask of samples that express 29 SNP fragments of the 193 SNPs. Only samples identified in the input mask are elgible for the output mask.
In this case 93 of the chromsomes samples express the first 29 SNPs. That pattern reflects the result of some recombination event. The other chromosome sample expresses a very irregular collection of 29 different SNPs generated by some kind of unknown process.
plt_obj = dm.superset_allele_mask(am_29, min_match=0.5)
plt = plt_obj.do_plot()
show(plt)
This plot shows the series expressed by at least 50% of the chromosome samples that express a 29 SNP fragment of the series 193_843 SNPs. It indentifies the results of recombination events that have associated the series 4_815 from the South Asian tree with 5_684 and 32_1361 from the East Asian tree and with 64_1575, 10_2206, and 7_1868 from the European tree.
HTML(plt_obj.get_html())
The data from the plot show that 93 out of the 94 chromsome samples that express 29 193_843 SNPs express the series 4_815.
plt_obj = dm.superset_allele_mask(am_29, [dm.di_4_815], min_match=0.5)
plt = plt_obj.do_plot()
am_4_815 = plt_obj.plot_context.yes_allele_mask
show(plt)
HTML(plt_obj.get_html())
The data from the plot show that 87 of the 93 chromsome samples that express 29 193_843 SNPs and the series 4_815 also express the series 5_864.
an.sa_193_843.alleles_per_snp(am_4_815)
This data confirms that the 93 chromosome samples that express 4_815 are the 93 chromsome samples that express the first 29 SNPs of the series 193_843.
plt_obj = dm.superset_allele_mask(am_29, [dm.di_4_815, dm.di_5_684, dm.di_32_1361, dm.di_64_1575], min_match=0.001)
plt = plt_obj.do_plot()
show(plt)
This plot limits the chromosme samples to those that express all of the 6 series shown in the plot and no others.
HTML(plt_obj.get_html())
The data shows that 58 chromosomes satisfy the plot's criteria. All but 1 of them are East Asian.
The next 6 plots show the series expressed by the 6 chromosome samples that express the first 29 193_843 SNPs, series 4_815, but not series 5_684. Every one of those chromosome samples expresses a different association of series.
plt_obj = dm.superset_allele_mask(am_29, [dm.di_4_815, dm.di_32_1361, dm.di_81_857, dm.di_5_47], [dm.di_5_684], min_match=0.001)
plt = plt_obj.do_plot()
show(plt)
HTML(plt_obj.get_html())
plt_obj = dm.superset_allele_mask(am_29, [dm.di_4_815, dm.di_32_1361, dm.di_81_857], [dm.di_5_47, dm.di_5_684], min_match=0.001)
plt = plt_obj.do_plot()
show(plt)
HTML(plt_obj.get_html())
plt_obj = dm.superset_allele_mask(am_29, [dm.di_4_815, dm.di_32_1361], [dm.di_81_857, dm.di_5_684], min_match=0.001)
plt = plt_obj.do_plot()
show(plt)
HTML(plt_obj.get_html())
plt_obj = dm.superset_allele_mask(am_29, [dm.di_4_815, dm.di_5_1460, dm.di_117_1685],
[dm.di_32_1361, dm.di_5_684], min_match=0.001)
plt = plt_obj.do_plot()
show(plt)
HTML(plt_obj.get_html())
plt_obj = dm.superset_allele_mask(am_29, [dm.di_4_815, dm.di_5_1460],
[dm.di_117_1685, dm.di_32_1361, dm.di_5_684], min_match=0.001)
plt = plt_obj.do_plot()
show(plt)
HTML(plt_obj.get_html())
plt_obj = dm.superset_allele_mask(am_29, [dm.di_4_815], [dm.di_5_1460, dm.di_32_1361, dm.di_5_684], min_match=0.001)
plt = plt_obj.do_plot()
show(plt)
HTML(plt_obj.get_html())
The following plot shows the series expressed by the 1 chromosome sample that expresses 29 193_843 SNPs but does not express the series 4_815
plt_obj = dm.superset_allele_mask(am_29, [dm.di_117_1685], [dm.di_4_815], min_match=0.001)
plt = plt_obj.do_plot()
am_not_4_815 = plt_obj.plot_context.yes_allele_mask
show(plt)
HTML(plt_obj.get_html())
an.sa_193_843.alleles_per_snp(am_not_4_815)
This data confirms that this chromosome sample is the one that expresses the divergent series of 29 193_843 SNPs.
aps_170, am_170 = an.sa_193_843.snps_from_aps_value(170, nam_193_843)
aps_170
This data shows the expressed SNPs and for the chromosome samples that express 170 193_843 SNPs. In this case, all 41 of the chromosome samples express the first 170 SNPs.
plt_obj = dm.superset_allele_mask(am_170, min_match=0.5)
plt = plt_obj.do_plot()
show(plt)
This plot indicates that a selection process associated with the emergence of the series 9_39 has resulted in overexpression of the first 170 193_843 SNPs. The series 9_39 appears to be rooted in an hierachy formed by a recombination of 8_267 from the South Asian tree and the series including 32_1361 and 81_857 at the root of the upper region of the East Asian tree
HTML(plt_obj.get_html())
The data associated with the plot shows that 38 of the chromosome samples that express the series 9_38 express the first 170 SNPs of series 193_843.
aps_169, am_169 = an.sa_193_843.snps_from_aps_value(169, nam_193_843)
aps_169
This data is for the 1 chromosome sample that expresses 169 series 193_843 SNPs. It expresses all but 1 of the first 170 SNPs.
plt_obj = dm.superset_allele_mask(am_169, min_match=0.5)
plt = plt_obj.do_plot()
show(plt)
This plot shows that the 1 chromosome sample that expresses 169 series 193_843 SNPs is the remaining chromosome sample that expresses the series 9_39
HTML(plt_obj.get_html())
This data confirms that the 1 chromosome sample has matched all of the plotted series.
plt_obj = dm.superset_allele_mask(am_170, [dm.di_8_267, dm.di_81_857], [dm.di_9_39], min_match=0.001)
plt = plt_obj.do_plot()
show(plt)
This plot shows the series association for chromosome samples that do express 170 series 193_843 SNPs, the series 8_267, the series 81_857, but not the series 9_39
HTML(plt_obj.get_html())
This data shows that this plot accounts for 2 of the chromsome samples that express 170 series 193_843 SNPs but do not express series 9_39.
plt_obj = dm.superset_allele_mask(am_170, [dm.di_8_267], [dm.di_9_39, dm.di_81_857], min_match=0.001)
plt = plt_obj.do_plot()
show(plt)
This plot shows the series expressed by the 1 remaining chromosome sample that expresses 170 series 193_843 SNPs but does not express the series 9_39.
HTML(plt_obj.get_html())
The data confirms that the plotted association of series is expressed by 1 chromosome sample.
SNP series are documented with a table that shows the properties of the series and the regional population samples that express it. The individual series SNPs are documented in a separate table that shows the properties of each individual SNP.
The column in the SNP data tables labeled "niv" means not_expressed_is_variant. The SNP allele with the lowest frequency in the thousand genome data is always considered to be the variant even when the reference genome labeled the more common allele as the variant. An niv value of 1 means the allele considered as the variant for this analysis is the one that the reference genome considers to be the standard. Note that the word allele is used to indicate a single chromosome sample that expresses an individual SNP or an SNP series.
HTML(dm.series_html(dm.di_11_765))