3. DMP and DMR
[1]:
from pylluminator.samples import Samples
from pylluminator.visualizations import manhattan_plot_dmr, plot_dmp_heatmap, visualize_gene
from pylluminator.dm import get_dmp, get_dmr, get_top_dmrs
from pylluminator.utils import save_object, load_object
from pylluminator.utils import set_logger
set_logger('WARNING') # set the verbosity level, can be DEBUG, INFO, WARNING, ERROR
3.1. Load pylluminator Samples
We assume that you have already processed the .idat files according to your preferences and saved them. If not, please refer to notebook 1 - Read data and get beta values before going any further.
[2]:
my_samples = Samples.load('preprocessed_samples')
Here, we want to filter out the probes on the X or Y chromosomes.
[3]:
my_samples.mask_xy_probes()
To speed up the demo, we will only calculate DMP and DMR on 10% of the probes
[4]:
ten_pct_probes = int(0.1 * my_samples.nb_probes)
probe_ids = my_samples.probe_ids[:ten_pct_probes]
print(f'Selected {ten_pct_probes:,} first probes')
Selected 93,768 first probes
3.2. Differentially Methylated Probes
The second parameter of get_dmp() is a R-like formula used in the design matrix to describe the statistical model, e.g. ‘~age + sex’. The names must be the column names of the sample sheet provided as third parameter
More info on design matrices and formulas:
[5]:
my_samples.sample_sheet
[5]:
| sample_id | sample_name | sample_type | |
|---|---|---|---|
| 0 | GSM7698438 | LNCAP_500_1 | LNCAP |
| 1 | GSM7698446 | LNCAP_500_2 | LNCAP |
| 2 | GSM7698462 | LNCAP_500_3 | LNCAP |
| 3 | GSM7698435 | PREC_500_1 | PREC |
| 4 | GSM7698443 | PREC_500_2 | PREC |
| 5 | GSM7698459 | PREC_500_3 | PREC |
[6]:
dmps, contrasts = get_dmp(my_samples, '~ sample_type', probe_ids=probe_ids)
contrasts
[6]:
['sample_type[T.PREC]']
You can now plot the results, for the 25 most variable probes:
[7]:
plot_dmp_heatmap(dmps, my_samples, contrasts[0], nb_probes=25)
3.3. Differentially Methylated Regions
[8]:
dmrs = get_dmr(my_samples, dmps, contrasts, probe_ids=probe_ids)
[9]:
manhattan_plot_dmr(dmrs, contrasts[0], annotation=my_samples.annotation)
[10]:
get_top_dmrs(dmrs, my_samples.annotation, contrasts[0])
[10]:
| genes | |||
|---|---|---|---|
| segment_sample_type[T.PREC]_p_value | chromosome | segment_id | |
| 7.581671e-63 | 6 | 33585 | ENSG00000237669;HLA-J;ZNRD1ASP |
| 8.426901e-61 | 3 | 29009 | FOXL2NB;FOXL2;LINC01391 |
| 7.147789e-59 | 20 | 25761 | ISM1 |
| 6.430889e-58 | 5 | 31695 | PDE4D;ENSG00000247345 |
| 1.662641e-56 | 6 | 33849 | HSPA1A;HSPA1L |
| 5.611748e-56 | 11 | 6788 | ADM;ADM-DT;SBF2 |
| 8.651960e-55 | 6 | 33533 | OR2I1P |
| 1.668132e-51 | 11 | 8198 | PIWIL4-AS1;AMOTL1 |
| 2.642552e-51 | 7 | 36088 | ENSG00000276496;ENSG00000278020;ENSG0000027796... |
| 8.382595e-50 | 12 | 10408 | DNAH10 |
3.4. Gene visualization
Let’s have a look at a particular gene identified as differentially methylated
[11]:
visualize_gene(my_samples, 'FOXL2')
