<jupyter_start><jupyter_text>Using pandasGEXpress with .gct/x files -------------------------------------------------------------------------------- A sample tutorial using data from [GEO GSE70138](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE70138)Author: <NAME>Last updated March 13th, 2018For questions and issue reports, please file an issue at github.com/cmap/cmapPy. --------------------------------------------------------------------------------  Table of ContentsWe have developed this tutorial to loosely follow a workflow we recommend for .gct/x files downloaded from GEO GSE70138. [Background info](background)[Obtaining data from GEO](getdata)[Activate your conda environment](conda)[Parse in metadata only](metadataparse)[Parsing in a gct/x file](parsegctx)[Add/delete annotations](annotations)[Subsetting gct/x files and GCToo instances](subsetting)[Merging gct/x files](merging)[Construct a GCToo instance](newGCToo)[Using GCToo instances with pandas tools](pandas)[Plotting with GCToo objects](plotting)[Write a GCToo instance to file](writing)[Converting .gct  .gctx](converting)--------------------------------------------------------------------------------  Background: What's .gct, .gctx, and GCToo? For more details on each of these topics, please visit clue.io/code.  The .gct format ![Sample GCT file](GCT_mockup.png)The GCT file format, a tab-delimited text-based format pairing matrix expression values with row and column metadata, allowing comparison of both transcriptional and contextual differences across samples. A schematic of a sample GCT file is pictured above.  The .gctx format![Sample GCTX file](GCTX_mockup.png)The GCTX file format was developed to improve serialization performance of GCT, particularly for large files. GCTX is an HDF5-based data format that can contain the same content as one or more GCT files in a hierarchically-structured, compressed, and easily accessible form (full file specification available at clue.io/code). Amongst other features, the HDF5 format allows users the ability to parse a subset of a large data set without reading the entire data set into memory first (using hyperslab selection). A schematic of a sample GCTX file is pictured above. The GCToo object `GCToo` objects are instances of the GCToo class we developed to represent .gct and .gctx files in a form compatible with Python's [PANDAS](http://pandas.pydata.org/) library. In essence, GCToo is a container of 3 pandas DataFrames (optionally, can also include a multi-index DataFrame of these 3 DataFrames combined) as well as two attributes representing the file source and version respectively. For more information, please see [our documentation](https://clue.io/cmapPy/pandasGEXpress.htmlgctoo-class).   Obtaining Data from Gene Expression Omnibus (GEO) [Gene Expression Omnibus](https://www.ncbi.nlm.nih.gov/geo/), or GEO, is a NCBI-maintained public repository of array and sequence genomics data. As of April 2017, we have provided 1.7 million profiles to GEO at the following accessions:* [GSE70138](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE70138)* [GSE92742](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE92742)* [GSE92743](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE92743)Please note that this does not represent the only data available in .gct and .gctx format; we simply mention it as a public resource for freely available data in these formats, and will use it as an example in this tutorial. Data can be downloaded from GEO using either using FTP or HTTP; see accession pages above for details. Details on what content is within each file published is clarified in the [GEO User Guide](https://docs.google.com/document/d/1q2gciWRhVCAAnlvF2iRLuJ7whrGP6QjpsCMq1yWz7dU/edit?usp=sharing).   Activate your cmapPy conda environment If you haven't done so already, [create a conda environment to use for your analyses](https://clue.io/cmapPy/build.htmlinstall).For this example, we've named our conda environment ``my_cmapPy_env``. **Be sure to activate your environment before proceeding!**``source activate my_cmapPy_env`` What version of cmapPy does this tutorial use?<jupyter_code>%load_ext autoreload
%autoreload 2

import pkg_resources

# Print version of cmapPy being used in current conda environment 
pkg_resources.get_distribution("cmapPy").version
# downloads example data, ~5Gb
import urllib.request
import os

urls = ['https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE70138&format=file&file=GSE70138%5FBroad%5FLINCS%5Fsig%5Finfo%5F2017%2D03%2D06%2Etxt%2Egz',
        'https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE70138&format=file&file=GSE70138%5FBroad%5FLINCS%5Fgene%5Finfo%5F2017%2D03%2D06%2Etxt%2Egz',
        'https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE70138&format=file&file=GSE70138%5FBroad%5FLINCS%5FLevel5%5FCOMPZ%5Fn118050x12328%5F2017%2D03%2D06%2Egctx%2Egz']
files = ['GSE70138_Broad_LINCS_sig_info_2017-03-06.txt.gz',
         'GSE70138_Broad_LINCS_gene_info_2017-03-06.txt.gz',
         'GSE70138_Broad_LINCS_Level5_COMPZ_n118050x12328_2017-03-06.gctx.gz']


for url,filename in zip(urls,files):
    if not os.path.isfile(filename):
        urllib.request.urlretrieve(url, filename)

        
import gzip
import shutil
if not os.path.isfile('GSE70138_Broad_LINCS_Level5_COMPZ_n118050x12328_2017-03-06.gctx'):
    with gzip.open('GSE70138_Broad_LINCS_Level5_COMPZ_n118050x12328_2017-03-06.gctx.gz', 'rb') as f_in:
        with open('GSE70138_Broad_LINCS_Level5_COMPZ_n118050x12328_2017-03-06.gctx', 'wb') as f_out:
            shutil.copyfileobj(f_in, f_out)<jupyter_output><empty_output><jupyter_text>Parse in metadata only Although we are in the process of developing tools to provide out of core functionality for larger that memory .gctx files, for now some .gctx files may be larger than the RAM you have available on your computer (note, we do not recommend the use of .gct for files above ~1000 samples). If this is the case, you have two options:    1. Obtain access to a server/instance that has more RAM.     2. Read the metadata in first, find the sample (column) id values ("cids") that you're actually interested in, and only parse those cids in.  How to parse in metadata only from GEOThis would be applicable if you were interested in only certain samples (or certain compound names, cell types, etc). Here, we parse in the metadata from the level 5 (replicate-collapsed signature) data on GSE70138. Sample (column) annotations On GEO, because of submission requirements annotations file are stored as separate tab-delimited text files. For level 5 data, we're interested in the "sig" related data tables. For this example, we previously downloaded and unzipped "GSE70138_Broad_LINCS_sig_info_2017-03-06.txt.gz" (corresponding to identifying information about the samples; if you wanted metrics-related information for filtering, explore the equivalently labelled "sig metrics" file). Then, we read this in to a pandas DataFrame:<jupyter_code>import pandas as pd

sig_info = pd.read_csv("GSE70138_Broad_LINCS_sig_info_2017-03-06.txt.gz", sep="\t")
# look at the metadata fields provided 
sig_info.columns<jupyter_output><empty_output><jupyter_text>As an illustrative example, say we're only interested in looking at expression data corresponding to cells treated as vorinostat. So, we first want to find the sig_id values that correspond to vorinostat-treated samples:<jupyter_code># get all samples (across all cell types, doses, and other treatment conditions) that were treated with vorinostat
vorinostat_ids = sig_info["sig_id"][sig_info["pert_iname"] == "vorinostat"]
# how many samples are there in this data set?
print("number of samples treated with vorinostat:")
len(vorinostat_ids)<jupyter_output>number of samples treated with vorinostat:<jupyter_text>Gene (row) annotations Analogously, you might be interested in only parsing in a specific row space; for example, you might be interested at only only examining patterns in expression of [landmark genes](https://clue.io/connectopedia/what_are_landmark_genes). The metadata column identifying landmark membership is `pr_is_lm` (for more details on other metadata header meanings, please see the [GEO User Guide](https://docs.google.com/document/d/1q2gciWRhVCAAnlvF2iRLuJ7whrGP6QjpsCMq1yWz7dU/editheading=h.l6bq0r1aih50)). So, we would parse in the gene annotations and obtain a list of only the row identifiers corresponding to landmark genes. **Note that for Entrez ids (as below), you need to specify that the data type annotations are strings (dtype=str argument below); otherwise, the parser currently behaves as if you're subsetting by index.**<jupyter_code>gene_info = pd.read_csv("GSE70138_Broad_LINCS_gene_info_2017-03-06.txt.gz", sep="\t", dtype=str)
gene_info.columns
# again, note that because gene_info dtypes are str, you need to find ids equal to "1" not 1 
landmark_gene_row_ids = gene_info["pr_gene_id"][gene_info["pr_is_lm"] == "1"]<jupyter_output><empty_output><jupyter_text>Once you have the row and/or column identifiers you're interested in, parse only these ids from the GCTX file (see [Parse in a gctx](parsegctx) below).   Side note: How to parse in metadata only from a metadata-annotated gctxNote. *This is not specifically relevant to the data available on GEO*, but if you have a .gctx file with row/column annotations and want to only look at the metadata annotations, you can also read only these fields in from the gctx. You would do this as follows:``from cmapPy.pandasGEXpress.parse import parse`` ``gctx_metadata_only = parse("path/to/gctx.gctx", metadata_only=True)``<jupyter_code>from cmapPy.pandasGEXpress.parse import parse

my_col_metadata = parse("GSE70138_Broad_LINCS_Level5_COMPZ_n118050x12328_2017-03-06.gctx", col_meta_only=True)
my_col_metadata[0:5]<jupyter_output><empty_output><jupyter_text>Parse in a gctxNow we're ready to read in expression data! For parsing, here are your options:1. Read in an entire gctx file2. Read in a subset of ids Generally, we recommend choice 2.  Example 1: Parse in a subset of samples/columns We'll demonstrate how to do so here using the ``vorinostat_ids`` we obtained [in the previous section](metadataparse).<jupyter_code>from cmapPy.pandasGEXpress.parse import parse
vorinostat_only_gctoo = parse("GSE70138_Broad_LINCS_Level5_COMPZ_n118050x12328_2017-03-06.gctx", cid=vorinostat_ids)
# Check to make sure we have all 210 colums of expression data: 
vorinostat_only_gctoo.meth_df.shape #yup!<jupyter_output><empty_output><jupyter_text>Example 2: Parse in a subset of genes/rowsWe'll demonstrate how to do so here using the ``landmark_gene_row_ids`` we obtained [in the previous section](metadataparse).<jupyter_code>landmark_only_gctoo = parse("GSE70138_Broad_LINCS_Level5_COMPZ_n118050x12328_2017-03-06.gctx", rid = landmark_gene_row_ids)
# Check to make sure we only have 978 rows of expression data: 
print(landmark_only_gctoo.cov_df.shape, landmark_only_gctoo.meth_df.shape) #yup!<jupyter_output>(978, 118050) (978, 118050)<jupyter_text>Notes:* If you did want to read in the entire file (option 1 above), you would just parse in the gct/x without the ``cid`` argument.* Subsetting as demonstrated above is also available for row ids (rids). Here, this would correspond to only looking at certain genes of interest.   Add/delete annotations to an existing GCToo instance  Adding annotationsNow we have our sig_id annotations, and also a GCToo instance of expression data. This would probably be more useful to us as one unified object, so let's annotate our `vorinostat_only_gctoo` with the id information relevant to it:<jupyter_code># first, we need to subset all the metadata information from our larger metadata DataFrame
vorinostat_sig_id_info =  sig_info[sig_info["pert_iname"] == "vorinostat"]
# let's look at a couple rows! 
vorinostat_sig_id_info[0:5]
# to add this data as our col_metadata_df, 
# the sig_id column needs to be the index of the metadata DataFrame
# For details on why this is, please refer to GCTX requirements: https://clue.io/connectopedia/gctx_format
vorinostat_sig_id_info.set_index("sig_id", inplace=True)
# now the data frame is indexed by sig_ids consistent with those of the data_df:
vorinostat_sig_id_info.index
# set the relevant annotations as a col_metadata_df: 
vorinostat_only_gctoo.col_metadata_df = vorinostat_sig_id_info<jupyter_output><empty_output><jupyter_text>Deleting Annotations Hm, but now that we look at it we probably don't need the "pert_type" field; since we're only looking at vorinostat data, it's "trt_cp" (compound) for all samples. We could keep it anyway, but just for fun let's delete it:<jupyter_code>del vorinostat_only_gctoo.col_metadata_df["pert_type"]
vorinostat_only_gctoo.col_metadata_df[0:5]<jupyter_output><empty_output><jupyter_text>Subsetting GCToo instances Let's say that we're actually only interested in looking at vorinostat's signatures in A375 cells. First, are there vorinostat samples treated in this cell line in this GCToo instance? Let's check by first listing all unique cell lines present:<jupyter_code>vorinostat_only_gctoo.col_metadata_df["cell_id"].unique()<jupyter_output><empty_output><jupyter_text>Let's subset our GCToo instance to only include samples/columns corresponding to vorinostat-treated A375 cells.<jupyter_code>import cmapPy.pandasGEXpress.subset_gctoo as sg

# first, get the cid values we want to keep 
vorinostat_a375_ids = vorinostat_only_gctoo.col_metadata_df.index[vorinostat_only_gctoo.col_metadata_df["cell_id"] == "A375"]

# subset GCToo instance to only include these values 
vorinostat_a375_gctoo = sg.subset_gctoo(vorinostat_only_gctoo, cid=list(vorinostat_a375_ids.values))

# Now, our samples only correspond to one cell type: A375
vorinostat_a375_gctoo.col_metadata_df["cell_id"].unique()<jupyter_output><empty_output><jupyter_text>Note that you can also subset by excluding rids/cids. For instance, if we wanted to include all cell lines **except** for A375, we could do so as follows:<jupyter_code># subset vorinostat_only_gctoo to exclude A375 samples only
vorinostat_except_a375_gctoo = sg.subset_gctoo(vorinostat_only_gctoo, exclude_cid=vorinostat_a375_ids)

# Now we have all cell lines except for A375: 
vorinostat_except_a375_gctoo.col_metadata_df["cell_id"].unique()<jupyter_output><empty_output><jupyter_text>You could subset genes/rows of the GCToo analogously.   Merging GCToo objects/.gct/.gctx files If you'd like to concatenate two or more objects/files, you have two options: 1. Concatenate files from the command line 2. Concatenate GCToo objects within a Python session. To continue with our GEO-based example, here we'll demonstrate option 2; for details on how to concatenate on the command line, please see [our documentation](https://clue.io/cmapPy/pandasGEXpress.html). Concatenating GCToo objects in a Python session We have our vorinostat expression data in A375; let's say we want to also get belinostat expression data in A375 (separately, for purposes of demonstration) and then combine these two GCToos. First, some prep: let's get the belinostat expression data in A375 that we need to do this operation:<jupyter_code># get sig info corresponding to belinostat-treated A375 samples
belinostat_a375_sig_info = sig_info[(sig_info["pert_iname"] == "belinostat") & (sig_info["cell_id"] == "A375")]
# set sig_id as index 
belinostat_a375_sig_info.set_index("sig_id", inplace=True)
# parse in only expression data corresponding to belinostat-treated A375 samples
belinostat_a375_gctoo = parse("GSE70138_Broad_LINCS_Level5_COMPZ_n118050x12328_2017-03-06.gctx", cid=belinostat_a375_sig_info.index)
# annotate belinostat_a375_gctoo with its corresponding sample annotations
belinostat_a375_gctoo.col_metadata_df = belinostat_a375_sig_info.copy()<jupyter_output><empty_output><jupyter_text>And now we can concatenate our belinostat & vorinostat GCToos as follows:<jupyter_code>import cmapPy.pandasGEXpress.concat as cg 
# horizontally stack the two GCToo instances (note, you can also vertically stack)
HDAC_a375_gctoo = cg.hstack([vorinostat_a375_gctoo, belinostat_a375_gctoo])

# Now we have a gctoo corresponding to two compounds, in one cell line: 
print("unique compound data present: ") 
print(HDAC_a375_gctoo.col_metadata_df["pert_iname"].unique())
print("unique cell line data present: ")
print(HDAC_a375_gctoo.col_metadata_df["cell_id"].unique())

# Also note we still have the same genes as before, just with additional samples of data: 
print("combined gctoo expression data shape:")
print(HDAC_a375_gctoo.meth_df.shape,HDAC_a375_gctoo.cov_df.shape)<jupyter_output>unique compound data present: 
['vorinostat' 'belinostat']
unique cell line data present: 
['A375']
combined gctoo expression data shape:
(12328, 33) (12328, 33)<jupyter_text>Concatenating two .gct/x files on the command lineWe provide the command line tool `concat` to concatenate either .gcts or .gctx files on the command line. Sample usage: ` concat--concat_direction horiz --input_filepaths first_file.gct second_file.gct` will output a file named concated.gctx; you can specify a different outname if you prefer. On the command line, type `concat -h` to see help options.   Constructing a new GCToo instance programatically If you have your own DataFrames of data, you can also construct your own GCToo instance (note that to do so, you need to fulfill the [requirements of a valid gctx](clue.io/code)). Just for illustrative purposes, let's make one using copies of our HDAC GCToo DataFrames, so that we can also associate a single [multi-index DataFrame](http://pandas.pydata.org/pandas-docs/stable/advanced.html) of our combined data_df, row_metadata_df, and column_data_df with this GCToo instance. **Note**: By default, the GCToo constructor does not make a multi_index_df; details on this and other optional arguments [here](http://cmappy.readthedocs.io/en/latest/pandasGEXpress.htmlgctoo-class).<jupyter_code>import cmapPy.pandasGEXpress.GCToo as GCToo

HDAC_a375_gctoo_with_multiindex = GCToo.GCToo(meth_df=HDAC_a375_gctoo.meth_df.copy(), cov_df=HDAC_a375_gctoo.cov_df.copy(),
                                              row_metadata_df=HDAC_a375_gctoo.row_metadata_df.copy(), 
                                              col_metadata_df=HDAC_a375_gctoo.col_metadata_df.copy(), 
                                              make_multiindex=True)

# Let's take a look at a subset of this multiindex DataFrame; this can be useful for indexing etc.
HDAC_a375_gctoo_with_multiindex.multi_index_meth_df[0:5]

# Let's take a look at a subset of this multiindex DataFrame; this can be useful for indexing etc.
HDAC_a375_gctoo_with_multiindex.multi_index_cov_df[0:5]<jupyter_output><empty_output><jupyter_text>For more information on multi-index DataFrames, please refer to the PANDAS project documentation.   Using pandas tools with your GCToo instance A big reason why we wrote `pandasGEXpress` was to enable users to maximize their use of `pandas`, `scikit-learn`, and other powerful and robustly-developed Python data science packages. Interacting with such tools simply requires referencing the relevant DataFrame attributes of your GCToo instance; here, we'll demonstrate how to do some basic arithmetic and grouping tasks on components of a GCToo instance. It should go without saying that the examples below are common but also fairly trivial examples of the wide functionality that libraries like pandas and sci-kit learn provide. Another `pandas` functionality that is likely to be of use is its `groupby` tools; as long as your data/annotations aren't too big, it can work nicely with your GCToo instance's `multi_index_df`. Column/row-wise meansLet's go back to our vorinostat-treated A375 cells. What is the mean expression across all samples?<jupyter_code>import pandas as pd 

vorinostat_a375_gctoo.meth_df.mean(axis = 0) #easy-peasy<jupyter_output><empty_output><jupyter_text>If we wanted to determine the mean expression across all samples for individual genes, we would set axis=1.  Pairwise correlation of columns How correlated are our samples?<jupyter_code>vorinostat_a375_cors = vorinostat_a375_gctoo.meth_df.corr()
vorinostat_a375_cors[0:5] #only showing first 5 rows<jupyter_output><empty_output><jupyter_text>Plotting with GCToo instances Because the attributes of a GCToo instance are pandas DataFrames, you can use any plotting libraries compatible with pandas DataFrames to generate plots; we frequently use [seaborn](https://seaborn.pydata.org/) and [matplotlib](https://matplotlib.org/). For simple plots, pandas also provides built-in plotting functionality for DataFrame objects.  As an illustrative example, we calculate the correlations across samples again, but this time also plot it as a heatmap--all in one!<jupyter_code>import seaborn as sns
import matplotlib.pyplot as plt
sns.set()


sns.heatmap(vorinostat_a375_gctoo.meth_df.corr(), cmap='plasma')
plt.show()<jupyter_output><empty_output><jupyter_text>Writing GCToo instances to file Here's how you would write your GCToo instance as a .gctx; note that if you'd like to write to .gct, the syntax is identical but you need to import ``write_gct`` instead.<jupyter_code>import cmapPy.pandasGEXpress.write_gctx as wg

# write GCToo instance to .gctx
wg.write(belinostat_a375_gctoo, "belinostat_a375")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08494280892555785, 'fraction_numerical': 0.031361335083442715, 'mean_word_length': 5.381393957523183, 'pattern_counts': {'":': 0, '<': 75, '<?xml version=': 0, '>': 75, 'https://': 17, 'lorem ipsum': 0, 'www.': 8, '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': '19', 'dir': 'jupyter-structured-clean-dedup', 'id': '35412', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 7180, 'n_tokens_neox': 6493, 'n_words': 2294, 'path': 'tutorials/cmapPy_pandasGEXpress_tutorial.ipynb', 'repo_name': 'Cellular-Longevity/cmapPy', 'size': '21337'}