{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# The NIH LINCS data is very large\nThe dataset contains 12328 rows (genes) by 118050 columns (samples) for a total of 1,455,320,400 entries\nIf we could use this dataset as a training dataset for our models, we could train for much longer, and then use the Kaggle data just for finetuning. \n\n*By the the way, I am looking for a team. If you think I'd be a good fit, feel free to reach out*\n\nYou can download all relevant files on the [GEO Website](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE70138).\nFor this notebook, you'll need:\n* Four metadata files:\n    * GSE70138_Broad_LINCS_pert_info_2017-03-06.txt\n    * GSE70138_Broad_LINCS_sig_info_2017-03-06.txt\n    * GSE70138_Broad_LINCS_cell_info_2017-04-28.txt\n    * GSE70138_Broad_LINCS_gene_info_2017-03-06.txt\n* The L5 LINCS dataset: \n    * GSE70138_Broad_LINCS_Level5_COMPZ_n118050x12328_2017-03-06.gctx","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can't load in the entire L5 LINCS dataset on most machines, because it wouldn't fit in RAM:","metadata":{}},{"cell_type":"code","source":"# Assuming an 8 byte float\nbase_rows = 12328\nbase_columns = 118050\nf\"Thats {base_rows * base_columns * 8 / 1e9} gigabytes\"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**This is too large for practical use**, and might be too big to fit in working memory (RAM) for many computers.\n\nBut we aren't interested in all the samples. Many samples here measure data for experiments we are not interested in. Luckily, the NIH provides a number of metadata files we can use to help use decide which experiments we are interested in.\n\n### Determining data of interest\n\nThe reagent (chemical or genetic) being studied in a given experiment is called the perturbagen, and there a variety of types of perturbagens. \n\nLet's load in the metadata file that contains info on the perturbagens.","metadata":{}},{"cell_type":"code","source":"pert_info = pd.read_csv(\"data/GSE70138_Broad_LINCS_pert_info_2017-03-06.txt\", sep=\"\\t\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's see what kinds of perturbagens there are in the dataset:","metadata":{}},{"cell_type":"code","source":"pert_info[\"pert_type\"].drop_duplicates()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Looking at the Connectopedia [entry on perturbagens](https://clue.io/connectopedia/perturbagen_types_and_controls), we can see that the `pert_type` for drugs (compounds) is `trt_cp`. \n\nLet's see some examples:","metadata":{}},{"cell_type":"code","source":"# Also want to include controls\ncompound_perturbagens = pert_info[pert_info[\"pert_type\"].isin([\"trt_cp\",\"ctl_vehicle\"])]\nprint(f\"Found {compound_perturbagens.shape[0]} different compounds.\")\ncompound_perturbagens[:5]","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Note that `pert_iname` in this dataset corresponds with `sm_name` in the Kaggle dataset (`de_train.parquet`). The same holds true `canonical_smiles` and `SMILES`, respectively.\n\nBy the way, the negative control (DMSO) exists in this dataset too, with a special `pert_type` called `ctl_vehicle`.","metadata":{}},{"cell_type":"code","source":"control_perturbagen = pert_info[pert_info[\"pert_type\"]==\"ctl_vehicle\"]\ncontrol_perturbagen","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"And so are the positive controls: dabrafenib and belinostat.","metadata":{}},{"cell_type":"code","source":"positive_perturbagens = pert_info[pert_info[\"pert_iname\"].isin([\"dabrafenib\",\"belinostat\"])]\npositive_perturbagens","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Building an index\nThe work we just did tells us which `pert_id`'s we are interested in, but we don't quite have an index into the dataset yet. \n\nLet's load in the sample metadata. Note that there is 1 entry for every sample in the dataset.","metadata":{}},{"cell_type":"code","source":"sig_info = pd.read_csv(\"data/GSE70138_Broad_LINCS_sig_info_2017-03-06.txt\", sep=\"\\t\")\nsig_info.shape","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This sample metadata contains the same `pert_type` column as above perturbagen metadata, but it doesn't have any info on what the perturbagen is:","metadata":{}},{"cell_type":"code","source":"sig_info.columns","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Luckily we can combine all the work we've done so far. \nWe want to annotate every sample which is uses a compound perturbagen with it's SMILES and International Chemical Identifier (INCHI).\n\n*Hold on tight!*","metadata":{}},{"cell_type":"code","source":"compound_perturbagens = compound_perturbagens[['pert_id', 'canonical_smiles']]\nkey = \"pert_id\"\ncomp_sig_info = sig_info.join(compound_perturbagens.set_index(key),on=key,how=\"right\")\ncomp_sig_info","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Lookss great, but I want info about the cells. Let's load in the metadata.","metadata":{}},{"cell_type":"code","source":"gene_info = pd.read_csv(\"data/GSE70138_Broad_LINCS_cell_info_2017-04-28.txt\", sep=\"\\t\")\ngene_info","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Out of curiosity, do we have any Henrietta Lacks cells?","metadata":{}},{"cell_type":"code","source":"gene_info[gene_info[\"base_cell_id\"]==\"HELA\"]","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Fascinating. Regardless, lets kinds of cells are available:","metadata":{}},{"cell_type":"code","source":"gene_info[\"cell_type\"].drop_duplicates()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's dig into the primary and differentiated cells:","metadata":{}},{"cell_type":"code","source":"gene_info[gene_info[\"cell_type\"].isin([\"differentiated\",\"primary\"])]","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"No blood cells. Let's look a little further:","metadata":{}},{"cell_type":"code","source":"gene_info[gene_info[\"primary_site\"]==\"blood\"]","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Not helpful. Let's just keep the info we gathered about the compounds and move on.\nDid we save any space?","metadata":{}},{"cell_type":"code","source":"f\"We still have {12328 * comp_sig_info.shape[0] * 8 / 1e9} gigabytes\"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"But we can include only the genes we care about.\nSee this [notebook for more info](https://www.kaggle.com/code/laurasisson/exploring-the-lincs-gene-metadata).","metadata":{}},{"cell_type":"code","source":"gene_info = pd.read_csv(\"data/GSE70138_Broad_LINCS_gene_info_2017-03-06.txt\", sep=\"\\t\")\ngene_info","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_df = pd.read_parquet(\"data/de_train.parquet\")\ntrain_df","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"shared_gene_info = gene_info[gene_info[\"pr_gene_symbol\"].isin(train_df)]\nshared_gene_info","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's see how much space now:","metadata":{}},{"cell_type":"code","source":"f\"We still have {shared_gene_info.shape[0] * comp_sig_info.shape[0] * 8 / 1e9} gigabytes\"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It doesn't work (below). Let's do landmarks.","metadata":{}},{"cell_type":"code","source":"# from cmapPy.pandasGEXpress.parse import parse\n\n# gene_ids = shared_gene_info[\"pr_gene_id\"].astype(str)\n# sig_ids = comp_sig_info[\"sig_id\"]\n\n# l5_data = parse(\"data/GSE70138_Broad_LINCS_Level5_COMPZ_n118050x12328_2017-03-06.gctx\", cid = sig_ids, rid = gene_ids)\n# l5_data.data_df.shape","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"landmark_info = gene_info[gene_info[\"pr_is_lm\"]==1]\nlandmark_gene_row_ids = gene_info[\"pr_gene_id\"][gene_info[\"pr_is_lm\"] == 1]\nlandmark_info","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Did that save space?","metadata":{}},{"cell_type":"code","source":"f\"We have {landmark_info.shape[0] * comp_sig_info.shape[0] * 8 / 1e9} gigabytes\"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Great! Let's load in the dataset:","metadata":{}},{"cell_type":"code","source":"from cmapPy.pandasGEXpress.parse import parse\ngene_ids = landmark_info[\"pr_gene_id\"].astype(str)\nsig_ids = comp_sig_info[\"sig_id\"]\nl5_data = parse(\"data/GSE70138_Broad_LINCS_Level5_COMPZ_n118050x12328_2017-03-06.gctx\", cid = sig_ids, rid = gene_ids)\nl5_data.data_df.shape","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's prepare metadata for annotation.","metadata":{}},{"cell_type":"code","source":"for sk in [\"pert_id\",\"pert_type\",\"distil_id\"]:\n    del comp_sig_info[sk]\ncomp_sig_info.set_index(\"sig_id\", inplace=True)\ncomp_sig_info","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for gk in [\"pr_gene_title\",\"pr_is_bing\",\"pr_is_lm\"]:\n    del landmark_info[gk]\nlandmark_info.set_index(\"pr_gene_id\",inplace=True)\nlandmark_info.index = landmark_info.index.map(str)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Time to annotate!","metadata":{}},{"cell_type":"code","source":"l5_data.col_metadata_df = comp_sig_info\nl5_data.row_metadata_df = landmark_info\nl5_data.data_df","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Before we finish the labels, let's add a control column","metadata":{}},{"cell_type":"code","source":"# Add a control column\ncomp_sig_info[\"control\"] = comp_sig_info[\"pert_iname\"] == \"DMSO\"\ncomp_sig_info[comp_sig_info[\"control\"]]","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"These are some heavy duty operations, so let's add tqdm for progress bars","metadata":{}},{"cell_type":"code","source":"from tqdm.auto import tqdm\ntqdm.pandas()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Join on the index (experiment ID):","metadata":{}},{"cell_type":"code","source":"final_data = comp_sig_info.join(l5_data.data_df.T)\nfinal_data = final_data.rename(columns=landmark_info.to_dict()[\"pr_gene_symbol\"])\nfinal_data","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Build individual datasets","metadata":{}},{"cell_type":"markdown","source":"Calculate common genes","metadata":{}},{"cell_type":"code","source":"LINCS_TSM_IDX = 6\nlincs_genes = final_data.columns[LINCS_TSM_IDX:]\nlincs_genes","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"TRAIN_TSM_IDX = 5\nshared_genes = train_df.columns[train_df.columns.isin(lincs_genes)]\nshared_genes","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import tqdm\nimport torch","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"lincs_cmpd_df = final_data[~final_data[\"control\"]]\nlincs_control_pert = final_data[final_data[\"control\"]][shared_genes.union([\"cell_id\"])].groupby(by=\"cell_id\").mean()\nlincs_control_pert","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Join each value with the average of the controls, by cell_id\nlincs_joined_df = lincs_cmpd_df.join(lincs_control_pert,on=\"cell_id\",lsuffix=\" post_treatment\",rsuffix=\" pre_treatment\")\n# Split labels to use join suffix as top level for multindex. For non, overlapping columns, use 'label'\nlincs_joined_df.columns = pd.MultiIndex.from_tuples([(y if not pd.isnull(y) else \"label\",x) for (x,y) in lincs_joined_df.columns.str.split(expand=True)])\n# Sort each multindex\nlincs_joined_df = lincs_joined_df.sort_index(axis=1,level=[0,1])\n# The non-shared genes are kept in the df under \"label\". Drop them.\nlincs_joined_df.drop(lincs_joined_df[\"label\"].select_dtypes('number').columns, axis = 1, level=1,inplace = True)\nlincs_joined_df","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Time to build the datasets. There could be a way to do this vectorized but its ok.","metadata":{}},{"cell_type":"code","source":"kaggle_cmpd_df = train_df[~train_df[\"control\"]]\nkaggle_control_pert = train_df[train_df[\"control\"]][shared_genes.union([\"cell_type\"])].groupby(by=\"cell_type\").mean()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Join and label as above\nkaggle_joined_df = kaggle_cmpd_df.join(kaggle_control_pert,on=\"cell_type\",lsuffix=\" post_treatment\",rsuffix=\" pre_treatment\")\n# Split labels to use join suffix as top level for multindex. For non, overlapping columns, use 'label'\nkaggle_joined_df.columns = pd.MultiIndex.from_tuples([(y if not pd.isnull(y) else \"label\",x) for (x,y) in kaggle_joined_df.columns.str.split(expand=True)])\n# Sort each multindex\nkaggle_joined_df = kaggle_joined_df.sort_index(axis=1,level=[0,1])\n# The non-shared genes are kept in the df under \"label\". Drop them.\nkaggle_joined_df.drop(kaggle_joined_df[\"label\"].select_dtypes('number').columns, axis = 1, level=1,inplace = True)\nkaggle_joined_df","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_df = pd.read_csv(\"data/id_map.csv\")\ntest_joined_df = test_df.join(kaggle_control_pert,on=\"cell_type\")\n# Because there are no overlapping cells, we must create multiindex manually\ngene_cols = shared_genes + \" pre_treatment\"\nlabel_cols = test_joined_df.columns[~test_joined_df.columns.isin(shared_genes)] + \" label\"\n\ntest_joined_df.columns = pd.MultiIndex.from_tuples([(y,x) for (x,y) in label_cols.union(gene_cols,sort=False).str.split(expand=True)])\ntest_joined_df","metadata":{"editable":true,"slideshow":{"slide_type":""},"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's do some quick sanity checks","metadata":{}},{"cell_type":"code","source":"# Check column ordering is the same for pre_treatment\nassert (kaggle_joined_df[\"pre_treatment\"].columns ==  lincs_joined_df[\"pre_treatment\"].columns).all()\nassert (kaggle_joined_df[\"pre_treatment\"].columns ==  test_joined_df[\"pre_treatment\"].columns).all()\n# For post_treatment\nassert (kaggle_joined_df[\"post_treatment\"].columns ==  lincs_joined_df[\"post_treatment\"].columns).all()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Change some labels in the lincs data so we can compare more easily.\nlincs_joined_df = lincs_joined_df.rename(columns={\"canonical_smiles\":\"SMILES\",\"pert_iname\":\"sm_name\",\"cell_id\":\"cell_type\"})","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"lincs_joined_df.to_parquet(\"data/lincs_pretreatment.parquet\")\nkaggle_joined_df.to_parquet(\"data/kaggle_pretreatment.parquet\")\ntest_joined_df.to_parquet(\"data/test_pretreatment.parquet\")","metadata":{},"execution_count":null,"outputs":[]}]}