Open In Colab

ISB-CGC Community Notebooks

Check out more notebooks at our Community Notebooks Repository!

Title:   How to Perform Nearest Centroid Classification with_BigQuery
Author:  Lauren Hagen
Created: 2019-12-17
URL:     https://github.com/isb-cgc/Community-Notebooks/blob/master/Notebooks/How_to_preform_Nearest_Centroid_Classification_with_BigQuery.ipynb
Purpose: Demonnstrate using BigQuery to categorize patients based on gene expression using the Nearest Centroid Classification.
Notes:

Introduction

Overview

This notebook is to demonstrate how to use BigQuery to categorize patients based on gene expression using the Nearest Centroid Classifier. We will be using the Kidney renal papillary cell carcinoma (KIRP) study from The Cancer Genome Atlas (TCGA) as an example dataset and will use the RNA Sequence Gene Expression Data.

What is the Nearest Centroid Classifier?

Nearest Centroid Classifier assigns a label based on the nearest mean (centroid) of the training samples that the observation is closest to1. This classifier is an example of supervised learning and does not create a model for use later2.

Before we get started, we will need to load the BigQuery module, authenticate ourselves, create a client variable, and load necessary libraries.

In [0]:
# Load the BigQuery Module
from google.cloud import bigquery
In [2]:
# Authenticate ourselves
!gcloud auth application-default login
Go to the following link in your browser:

    https://accounts.google.com/o/oauth2/auth?code_challenge=vczbWqjPWlXINDmstsMXSC_dgiQ7OUkYLzS4LUn6f5o&prompt=select_account&code_challenge_method=S256&access_type=offline&redirect_uri=urn%3Aietf%3Awg%3Aoauth%3A2.0%3Aoob&response_type=code&client_id=764086051850-6qr4p6gpi6hn506pt8ejuq83di341hur.apps.googleusercontent.com&scope=https%3A%2F%2Fwww.googleapis.com%2Fauth%2Fuserinfo.email+https%3A%2F%2Fwww.googleapis.com%2Fauth%2Fcloud-platform+https%3A%2F%2Fwww.googleapis.com%2Fauth%2Faccounts.reauth


Enter verification code: 4/vwEmwq_ebAQ_L2OTt_nQjSnrLyC1fYuTMR1PH5CClYdxCP9qLhKUIdA
WARNING: Cannot find a project to insert into application default credentials (ADC) as a quota project.
Run $gcloud auth application-default set-quota-project to insert a quota project to ADC.

Credentials saved to file: [/content/.config/application_default_credentials.json]

These credentials will be used by any library that requests Application Default Credentials (ADC).


To take a quick anonymous survey, run:
  $ gcloud survey

In [0]:
# Create a variable for which client to use with BigQuery
project_num = 'your_project_number'# Replace with your project ID
if project_num == 'your_project_number':
    print('Please update the project number with your Google Cloud Project')
else:
    client = bigquery.Client(project_num)
In [0]:
# Required Libraries
import pandas
import pandas_confusion

Select Genes

We will be using the top three genes with the highest mean expression as the genes we will use to categorize the clinical stage for each case.

In [0]:
# Find top 20 genes by mean with training set
top_20_highest_mean_gene_exp = """SELECT gene_name, AVG(HTSeq__FPKM_UQ) as m
FROM `isb-cgc.TCGA_hg38_data_v0.RNAseq_Gene_Expression`
WHERE project_short_name = 'TCGA-KIRP'
GROUP BY gene_name
ORDER BY m DESC
LIMIT 20"""
In [0]:
# Run the query
top_20 = client.query(top_20_highest_mean_gene_exp).to_dataframe()
print(top_20)
   gene_name             m
0     MT-CO3  5.275677e+08
1     MT-CO1  4.638104e+08
2     MT-CO2  4.374221e+08
3     MT-ND4  4.188905e+08
4    MT-RNR2  3.383655e+08
5    MT-ATP6  2.715383e+08
6     MT-ND3  2.477042e+08
7     MT-CYB  2.397069e+08
8        FTL  1.631893e+08
9     MT-ND2  1.626764e+08
10    MT-ND1  1.602064e+08
11   MT-ND4L  1.299254e+08
12   MT-ATP8  1.064019e+08
13    MT-ND6  1.018257e+08
14    MT-ND5  8.410137e+07
15   MT-RNR1  8.201815e+07
16    TMSB10  7.268314e+07
17     MT-TP  6.806825e+07
18      SPP1  6.308363e+07
19  MTATP6P1  4.023022e+07

Nearest Centroid in BigQuery

What are we going to use Nearest Centroids for?

We will be attempting to predict the clinical stage for a case based on a specific subset of genes and their expression levels.

First, we will filter the TCGA clinical table for our target disease (KIRP) and not include the cases where the clinical stage is missing. We need there to be no missing clinical stage features due to the nature of k nearest algorithms.

In [0]:
filtered_clin = """
WITH
  clin AS (
  SELECT
    case_barcode,
    clinical_stage    
  FROM
    `isb-cgc.TCGA_bioclin_v0.Clinical`
  WHERE
    disease_code = 'KIRP'
    AND clinical_stage <> 'NULL'
  ), """

Next, we will want to get the gene expressions for the genes that we are going to use to attempt to identify the clinical stage. We will filter for the 3 genes that we want, then we will label each case with whether it will be in the training or testing group. For more information on randomization in BigQuery, please see the How to Create a Random Sample in BigQuery Notebook in the ISB-CGC Community Notebook Repository.

In [0]:
filtered_expr = """
expr AS (
  SELECT
    case_barcode,
    IF ( MOD(ABS(FARM_FINGERPRINT(case_barcode)),10) > 1, 'TRAIN', 'TEST') as class,
    gene_name,
    HTSeq__FPKM_UQ
  FROM
    `isb-cgc.TCGA_hg38_data_v0.RNAseq_Gene_Expression`
  WHERE
    project_short_name = 'TCGA-KIRP'
    AND gene_name IN ('MT-CO3','MT-CO1','MT-CO2')
    AND (
    case_barcode IN (select case_barcode from clin)
    )
  GROUP BY case_barcode, gene_name, HTSeq__FPKM_UQ
  ),"""

We will then find the centroids or means for each clinical stage within the training data.

In [0]:
centroids = """
 mean AS (
  SELECT
    expr.class,
    clin.clinical_stage,
    AVG(CASE
        WHEN gene_name = 'MT-CO3' THEN HTSeq__FPKM_UQ
    END
      ) AS gene1,
    AVG(CASE
        WHEN gene_name = 'MT-CO1' THEN HTSeq__FPKM_UQ
    END
      ) AS gene2,
    AVG(CASE
        WHEN gene_name = 'MT-CO2' THEN HTSeq__FPKM_UQ
    END
      ) AS gene3
  FROM
    expr
  JOIN
    clin
  ON
    expr.case_barcode = clin.case_barcode
  WHERE
    expr.class = 'TRAIN'
  GROUP BY
    expr.class,
    clin.clinical_stage),"""

We will also need the testing data to have the genes that we are going to use in the classifier in separate columns for the test cases.

In [0]:
test = """
test AS (
  SELECT
    case_barcode,
    class,
    MAX(CASE WHEN gene_name = 'MT-CO3' THEN HTSeq__FPKM_UQ END) AS gene1,
    MAX(CASE WHEN gene_name = 'MT-CO1' THEN HTSeq__FPKM_UQ END) AS gene2,
    MAX(CASE WHEN gene_name = 'MT-CO2' THEN HTSeq__FPKM_UQ END) AS gene3
  FROM
    expr
  WHERE
    class = 'TEST'
  GROUP BY case_barcode, class),"""

This next section of the query will find the euclidean distance for each case in the testing data set3. The euclidean distance can be found by the following equation4:

$d(q,p) = \sqrt{(q_1-p_1)^2+(q_2-p_2)^2+...+(q_n-p_n)^2}$

In [0]:
dist = """
distance AS (SELECT
  case_barcode,
  gene1,
  gene2,
  SQRT(POWER((SELECT gene1 FROM mean WHERE clinical_stage = 'Stage I')-gene1, 2) + POWER((SELECT gene2 FROM mean WHERE clinical_stage = 'Stage I')-gene2, 2) + POWER((SELECT gene3 FROM mean WHERE clinical_stage = 'Stage I')-gene3, 2)) as stage1,
  SQRT(POWER((SELECT gene1 FROM mean WHERE clinical_stage = 'Stage II')- gene1, 2) + POWER((SELECT gene2 FROM mean WHERE clinical_stage = 'Stage II')-gene2, 2) + POWER((SELECT gene3 FROM mean WHERE clinical_stage = 'Stage II')-gene3, 2)) as stage2,
  SQRT(POWER((SELECT gene1 FROM mean WHERE clinical_stage = 'Stage III')-gene1, 2) + POWER((SELECT gene2 FROM mean WHERE clinical_stage = 'Stage III')-gene2, 2) + POWER((SELECT gene3 FROM mean WHERE clinical_stage = 'Stage III')-gene3, 2)) as stage3,
  SQRT(POWER((SELECT gene1 FROM mean WHERE clinical_stage = 'Stage IV')-gene1, 2) + POWER((SELECT gene2 FROM mean WHERE clinical_stage = 'Stage IV')-gene2, 2) + POWER((SELECT gene3 FROM mean WHERE clinical_stage = 'Stage IV')-gene3, 2)) as stage4,
FROM test),
"""

Finally, we will calculate which category or stage has the smallest distance from the centroid and assign that stage to the case. We will also include the case barcode and actual clinical stage as columns in our final result to assist with calculating the accuracy of the classification.

In [0]:
pred = """
pred AS (SELECT case_barcode,
  (CASE WHEN stage1<stage2 AND stage1<stage3 AND stage1<stage4 THEN 'Stage I'
  WHEN stage2<stage1 AND stage2<stage3 AND stage2<stage4 THEN 'Stage II'
  WHEN stage3<stage1 AND stage3<stage2 AND stage3<stage4 THEN 'Stage III'
  ELSE 'Stage IV' END) AS prediction
FROM distance)

SELECT a.case_barcode, a.prediction, b.clinical_stage
FROM pred AS a
JOIN clin AS b
ON a.case_barcode = b.case_barcode
"""

Now that we have our full query laid out, we can combine all of the sub-queries into the full string and query the data in BigQuery.

In [15]:
# Build the query with all of the subqueries
final_query = filtered_clin + filtered_expr + centroids + test + dist + pred

# Run the query
cn = client.query(final_query).to_dataframe()
cn.head()
Out[15]:
case_barcode prediction clinical_stage
0 TCGA-BQ-7055 Stage I Stage I
1 TCGA-2Z-A9JK Stage IV Stage III
2 TCGA-EV-5903 Stage IV Stage I
3 TCGA-2Z-A9JE Stage I Stage I
4 TCGA-G7-6797 Stage IV Stage III

Find Accuracy of the Model

To find the accuracy of the model, we weil use a confusion matrix and then will calculate the accuracy with the following formula3:

$Accuracy = \frac{M_{ii}}{\sum_{j}M_{ji}}$

In [17]:
# Create a confusion matrix for the clinical predictions
confusion = pandas.crosstab(cn['prediction'], cn['clinical_stage'])
print(confusion)
clinical_stage  Stage I  Stage II  Stage III  Stage IV
prediction                                            
Stage I              12         0          1         2
Stage II              3         2          1         0
Stage III             1         0          0         0
Stage IV             13         2          3         0
In [39]:
# Calculate the accuracy of the model
print("Overall Accuracy is", (confusion['Stage I'][0] + 
                              confusion['Stage II'][1] + 
                              confusion['Stage III'][2] + 
                              confusion['Stage IV'][3]) / cn.count()[0])
print("Stage 1 Accuracy is", round(confusion['Stage I'][0]/sum(confusion['Stage I']), 2))
print("Stage 2 Accuracy is", round(confusion['Stage II'][1]/sum(confusion['Stage II']), 2))
print("Stage 3 Accuracy is", round(confusion['Stage III'][2]/sum(confusion['Stage III']), 2))
print("Stage 4 Accuracy is", round(confusion['Stage IV'][3]/sum(confusion['Stage IV']), 2))
Overall Accuracy is 0.35
Stage 1 Accuracy is 0.41
Stage 2 Accuracy is 0.5
Stage 3 Accuracy is 0.0
Stage 4 Accuracy is 0.0

It seems that this model has an overall accuracy of 35% and was able to get close to 41% accuracy for Stage 1 but it was not able to accurately predict any of the other classes. We can tell that this is not a good model for this data set though it could easily be updated to different genes or data sets.

References

1Statistical classification - Wikipedia. n.d. https://en.wikipedia.org/wiki/Statistical_classification.

2Euclidean distance - Wikipedia. n.d. https://en.wikipedia.org/wiki/Euclidean_distance.

3Finding Nearest Neighbors in SQL | Sisense. n.d. https://www.sisense.com/blog/finding-nearest-neighbors-in-sql/.

4Machine Learning with Python: Confusion Matrix in Machine Learning with Python. n.d. https://www.python-course.eu/confusion_matrix.php.