<jupyter_start><jupyter_text>Mixing local and remote processingIn this tutorial you will learn:1. How to connect to a remote GMQL instance and run remote queries2. How to interleave remote and local GMQL processing3. How to manipulate the GDataframe data structure to plot an heatmap4. How to perform a Principal Component Analysis using the GDataframe data structure5. How to transform a pure Pandas Dataframe back to a GMQLDataset to perform a new query  Connecting to a remote GMQL instancePyGMQL can interact with a remote server hosting a GMQL instance and repository.In order to connect to the remote instance the user has to:1. Set the IP address of the service2. Login to the serverIt must be noticed that GMQL offers both authenticated and un-authenticated login: by default PyGMQL executes an un-authenticated login. This means that the remote instance will provide PyGMQL with limited resources and not permanent result storage. In order to login to the GMQL server as an authenticated user, one has to:1. Register on the remote service by going to its IP address from the web browser2. Use the `pygmql_login` tool, which is installed together with the libraryIn the following, we are going to use the publicly accessible GMQL service offered by the Genomic Computing group at Politecnico di Milano (accessible at http://gmql.eu/gmql-rest/). Since we did not registered to the service, when we the `login` call is encountered, the library will login as GUEST user with limited capabilities (in terms of computation and disk space).<jupyter_code>import gmql as gl

gl.set_remote_address("http://gmql.eu/gmql-rest/")
gl.login()<jupyter_output>[PyGMQL] Logging using stored authentication token<jupyter_text>From this point on, all the remote calls (explained later) will be forwarded to the GMQL service hosted at the specified address as a guest user. Remote executionWith PyGMQL the user can decide if the next query (sequence of statements on GMQLDatasets ended with a materialization) will be executed locally or on the remote service.This can be done using the `set_mode` function. If the user sets the mode to `remote` the following will happen:1. The local datasets that are used are uploaded to the remote service. Nothing is done to the remote datasets used in the query (if present) since they are already on the server.2. A compressed representation of the query is sent to the remote service, decoded and executed3. Once the execution is complete, the results are downloaded, stored and loaded into a `GDataframe` . The following image describes the steps of remote processing with PyGMQL:  The exampleWe will showcase these features through a simple example. We want to understand the enrichment of histone modifications on promotorial regions of REFSEQ genes. In particular:* The gene RefSeq annotation dataset holds the information about gene symbols and coordinates and it is stored in the remote server* The Chip-Seq data about histone modifications will come from ENCODE and it is as well stored in the remote system.The result of the analysis will be an heatmap $genes \times histone\_mark$ displaying the average signal value of each HM for each of the selected genes.  Finding the promotorial regionsWe will use as definition of promoter the region situated at $\left[gene_{start} - 2000; gene_{start} + 2000\right]$.We will perform this operation remotely.<jupyter_code>gl.set_mode("remote")<jupyter_output><empty_output><jupyter_text>We "load" the remote RefSeq dataset using the `load_from_remote` function. Since this dataset is available to all the users of the system, we will specify `public` as the ownership of the datasets. In the case in which we are loading our private datasets stored on the service, we can omit the `owner` parameter.<jupyter_code>genes = gl.load_from_remote("HG19_ANNOTATION_REFSEQ", owner="public")<jupyter_output><empty_output><jupyter_text>Now we have to:1. Select only the samples from the `HG19_ANNOTATION_REFSEQ` dataset specific of genes2. Select only the regions related to the genes (TP53, MYC, MAX, MAZ, PTEN)This operation can be done using the `select` function on the `gene` GMQLDataset variable.<jupyter_code>genes = genes.select(meta_predicate=genes["annotation_type"] == 'gene')<jupyter_output><empty_output><jupyter_text>To extract the promotorial regions we will use the `reg_project` (region projection) operation. In particular the following call will:1. Return a dataset with only the `chr, start, stop, gene_symbol` region fields2. For each gene, return the genomic position at $\left[gene_{start} - 2000; gene_{start} + 2000\right]$ taking also into account the strand of the gene.<jupyter_code>promoters = genes.reg_project(field_list=['gene_symbol'], 
                              new_field_dict={'start': genes.start - 2000, 
                                              'stop': genes.start + 2000})<jupyter_output><empty_output><jupyter_text>Since we are interested in keeping the result of this part of the query for future usage, we will materialize the result locally.<jupyter_code>promoters = promoters.materialize("./promoters", all_load=False)<jupyter_output>RUNNING<jupyter_text>Selecting the histone marksNext we want to select the H3K27me3, H3K36me3, H3K4me1, H3K79me2, H3K9ac and H3K9me3 from the `HG19_ENCODE_BROAD_NOV_2017` dataset, also in this case stored remotely.As before, we will load the dataset using `load_from_remote`. The following operation is a *meta selection*, which means that the selection will be done on the metadata attributes of the dataset. The resulting dataset will have all the samples for which the metadata satisfy the specified predicate.In PyGMQL the metadata selection follows the same syntax as the Pandas Python data analysis library ```selected_dataset = dataset[dataset['metadata_attribute'] == 'value']```Since we want to test for multiple histone marks we can use the `isin()` operator.<jupyter_code>hms = gl.load_from_remote("HG19_ENCODE_BROAD_NOV_2017", owner="public")
hms = hms[hms['experiment_target'].isin(['H3K27me3-human', 'H3K36me3-human', 
                                         "H3K4me1-human", "H3K79me2-human", 
                                         "H3K9ac-human", "H3K9me3-human"])]<jupyter_output><empty_output><jupyter_text>Grouping replicatesThe dataset resulting from the previous selection will have several duplicates for each region, due to the several biological and technical replicates offered by ENCODE. Therefore we need a way to group together the overlapping regions and to aggregate the relative signal values.This can be done using the `normal_cover` operation. In particular this function enables you to specify:1. The number of minimum overlappings for a region to be kept2. The maximum number of overlappings for a region to be kept3. If you want to groupby the resulting samples by some metadata attribute4. The aggregation function for the input region fields to be applied to the region groupingsIn our case we will keep all the regions (minimum overallpping = 1, maximum overlapping = "ANY") and group by the `experiment_target` metadata attribute. The signals of overlapping regions will be averaged.<jupyter_code>hms = hms.normal_cover(1, "ANY", groupBy=['experiment_target'], 
                       new_reg_fields={'avg_signal': gl.AVG("signal")})<jupyter_output><empty_output><jupyter_text>Also in this case we want to keep the resulting dataset for later processing.<jupyter_code>hms = hms.materialize("./hms", all_load=False)<jupyter_output>RUNNING<jupyter_text>Mapping the histone marks signal to the promotorial regionsNow that we have our set of gene promoters and our set of histone marks, we want to map the signal of the latter to the regions of the former. We can do that using the `map` operation. In particular, since on a promotorial region can be mapped multiple HMs, we will aggregate the signal by averaging it.We decide to perform the following operations on the local machine.<jupyter_code>gl.set_mode("local")
mapping = promoters.map(hms, refName='prom', 
                        expName="hm", 
                        new_reg_fields={'avg_signal': gl.AVG("avg_signal")})<jupyter_output><empty_output><jupyter_text>We finally materialize the result and load it in a GDataframe.<jupyter_code>mapping = mapping.materialize("./mapping")
mapping.regs.head()
mapping.meta.head()<jupyter_output><empty_output><jupyter_text>Visualizing the HM signal heatmap on the gene promotersFinally we visualize the heatmap (in the form of clusters) with the assistance of the Seaborn visualization package.In order to visualize the data we have to transform the GDataframe data structure in the `mapping` variable in a matricial form where we have $genes$ on the rows and $histone\_marks$ on the columns. This can be done with the `to_matrix` function offered by the GDataframe.The `to_matrix` method is particularly versatile since it enables the user to project region and metadata at the same time. For example, in our case the information about the histone modification name is stored as a metadata attribute `hm.experiment_target` while we want to use the `avg_signal` region field as matrix value. The `to_matrix` methods therefore enable us to aggregate the information of both regions and metadata in a single structure. Underneeth this method leverages on the `pivot_table` Pandas function and can accept all its parameter together with its specific ones.<jupyter_code>!pip install seaborn
!pip install scikit-learn
import matplotlib.pyplot as plt
import seaborn as sns
plt.figure()
sns.clustermap(mapping.to_matrix(index_regs=['gene_symbol'], 
                                 columns_meta=['hm.experiment_target'], 
                                 values_regs=['avg_signal'], fill_value=0, dropna=False).sample(100))
plt.show()<jupyter_output><empty_output><jupyter_text>Performing Principal Component Analysis<jupyter_code>from sklearn.decomposition import PCA
import pandas as pd
pca = PCA(n_components=2)
promoter_region_VS_hm = mapping.to_matrix(index_regs=['chr', 'start', 'stop', 'gene_symbol'], 
                                                 columns_meta=['hm.experiment_target'], 
                                                 values_regs=['avg_signal'], fill_value=0)
components = pca.fit_transform(promoter_region_VS_hm.values)
PCA_regions = pd.DataFrame(index=promoter_region_VS_hm.index, columns=['1st_component', '2nd_component'], data=components).reset_index()
plt.figure(figsize=(10, 10))
sns.regplot(data=PCA_regions, x='1st_component', y='2nd_component', fit_reg=False)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.06461625282167043, 'fraction_numerical': 0.01109857035364936, 'mean_word_length': 4.756903086085544, 'pattern_counts': {'":': 0, '<': 50, '<?xml version=': 0, '>': 50, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '201446', 'license': 'Apache-2.0', 'n_tokens_mistral': 2903, 'n_tokens_neox': 2745, 'n_words': 1267, 'path': 'examples/notebooks/02b_Mixing_Local_Remote_Processing_COMPLETE.ipynb', 'repo_name': 'DEIB-GECO/PyGMQL', 'size': '10632'}