<jupyter_start><jupyter_text>The purpose of this notebook is to compare our median-VAF-based purity estimates with ones from the other well-established tools. The comparison uses the Cancer Genome Atlas (TCGA) where each method, depending on what data type it uses for the estimation, was run on the same set of patients. Purity estimation tools of interest here are: * `ESTIMATE`: gene expression data * `ABSOLUTE`: copy-number variation data * `LUMP`: methylation data * `IHC`: pathology images * `CPE`: consensus metric combining others (median value after normalization)All these estimates are coming from the following study: [Systematic pan-cancer analysis of tumour purity](http://www.nature.com/ncomms/2015/151204/ncomms9971/full/ncomms9971.html). We will be comparing these to our own simplistic metric, called `MPurity`, that we calculate by normalizing the median VAF for a sample by `0.5` (see [first post in the series](0 - What is a variant allele frequency?.ipynb) for details). For this, we will be using the [normalized MAF files from the TCGA project](https://www.biostars.org/p/161298/161777) that contain VAF information for each of the variants:<jupyter_code>!mkdir -p data/downloads \
    && cd data/downloads \
    && wget -c "https://dl.dropboxusercontent.com/u/8680991/mafs/tcga_pancancer_dcc_mafs_082115.tar.gz" \
    && tar zxf "tcga_pancancer_dcc_mafs_082115.tar.gz" \
    && cd ../..<jupyter_output>--2016-06-01 00:24:43--  https://dl.dropboxusercontent.com/u/8680991/mafs/tcga_pancancer_dcc_mafs_082115.tar.gz
Resolving dl.dropboxusercontent.com... 192.168.3.11
Connecting to dl.dropboxusercontent.com|192.168.3.11|:443... connected.
HTTP request sent, awaiting response... 206 PARTIAL CONTENT
Length: 354970179 (339M), 0 remaining [text/plain]
Saving to: 'tcga_pancancer_dcc_mafs_082115.tar.gz'

tcga_pancancer_dcc_ 100%[++++++++++++++++++++++] 338.53M  --.-KB/s   in 0s     

2016-06-01 00:24:44 (0.00 B/s) - 'tcga_pancancer_dcc_mafs_082115.tar.gz' saved [354970179/354970179]<jupyter_text>Now that we have all the normalized MAF files, let's extract the VAF information from studies that have this information and aggregate these values within a single data-frame:<jupyter_code>import glob
import pandas as pd

folder = "data/downloads/mafs" ## Extracted archive folder that contains all the MAFs
df = pd.DataFrame()
for maf_file in glob.glob("{}/*.maf".format(folder)):
    mdf = pd.read_csv(maf_file, delimiter="\t", comment="#", low_memory=False)
    if "t_alt_count" in mdf.columns:  ## Not all MAFs have this. If not, skip the study MAF
        # example file name: `tcga_laml_from_dcc.maf`
        file_name = maf_file.split("/")[-1].split("_")[1].upper()
        # Also corresponds to the study abbreviation
        mdf['Study'] = file_name
        # Shorten the ID to make it compatible with the other data set (see below)
        mdf['Tumor_Sample_Barcode'] = ["-".join(s.split("-")[0:4]) for s in mdf['Tumor_Sample_Barcode']]
        mdf['read_count'] = mdf.t_alt_count + mdf.t_ref_count
        mdf['VAF'] = mdf.t_alt_count / mdf.read_count
        # Only save columns that might be of interest to us
        mdf = mdf[['Study', 'Tumor_Sample_Barcode', 'VAF']]
        df = df.append(mdf)
df.head()<jupyter_output><empty_output><jupyter_text>Our `df` object contains VAF information for all the variants called from the patient data. But we need to summarize this data by taking the median value for each patient and normalizing it against the expected value: `0.5`:<jupyter_code>purities = df.groupby(['Study', 'Tumor_Sample_Barcode']).median()
purities['MPurity'] = [min(v, 1) for v in (purities.VAF / 0.5)]  ## min -> safety check for right-skewed data
purities['Study'] = [i[0] for i in purities.index]
purities['Sample'] = [i[1] for i in purities.index]
purities.index = purities.Sample
purities.head()<jupyter_output><empty_output><jupyter_text>and if we plot these purity estimated for each study, we get this:<jupyter_code>import seaborn as sns

%matplotlib inline
import matplotlib.pyplot as plt
sns.set(color_codes=True)

sns.factorplot(x="Study", y="MPurity", hue="Study", data=purities, kind="box", aspect=4.5)
plt.show()<jupyter_output><empty_output><jupyter_text>For comparison, let's load the TCGA estimates from various methods as well:<jupyter_code>tcga = pd.read_csv("data/2-ncomms9971_s2.csv")  ## From the supplementary material
tcga.index = tcga["Sample ID"]
tcga.head()
# Plot the consensus purity estimate
sns.factorplot(x="Cancer type", y="CPE", hue="Cancer type", data=tcga, kind="box", aspect=4.5)
plt.show()<jupyter_output><empty_output><jupyter_text>The plot above definitely looks somewhat different than ours, but before jumping to any conclusions, let's prepare the data for a fair pairwise-method comparison:<jupyter_code>all_ests = pd.concat([tcga.drop(["Sample ID", "Cancer type"], axis=1), purities["MPurity"]], axis=1, join="inner")
# We want to work with complete data, so let's drop the patients that lack at least one of the estimates
all_ests = all_ests.dropna()
all_ests.head()<jupyter_output><empty_output><jupyter_text>we can use this data-frame to compare all methods in a pairwise manner:<jupyter_code>sns.pairplot(all_ests, dropna=True, kind='reg')<jupyter_output><empty_output><jupyter_text>Surprising our simplistic `MPurity` estimates are not that bad (see all those correlated metrics), especially considering that we are doing this using a completely different data set than others: mutation data. To get a sense of which methods behave similarly and make estimates that are close to one another, we can calculate the Median Squared Error (MSE) value between two groups of estimates and use this a distance metric to a hierarchical clustering method:<jupyter_code>distance = lambda column1, column2: ((column1 - column2) ** 2).mean(axis=None)  ## MSE
mse_matrix = all_ests.apply(lambda col1: all_ests.apply(lambda col2: distance(col1, col2)))
sns.clustermap(mse_matrix)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11461461461461461, 'fraction_numerical': 0.03403403403403404, 'mean_word_length': 5.129856850715746, 'pattern_counts': {'":': 0, '<': 32, '<?xml version=': 0, '>': 33, 'https://': 3, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 2, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '139149', 'license': 'Apache-2.0', 'n_tokens_mistral': 2123, 'n_tokens_neox': 1941, 'n_words': 617, 'path': '2 - Estimating tumor purity from a VAF distribution.ipynb', 'repo_name': 'armish/vaf-experiments', 'size': '5992'}