<jupyter_start><jupyter_text>How to build an RNA-seq logistic regression classifier with BigQuery MLCheck out other notebooks at our [Community Notebooks Repository](https://github.com/isb-cgc/Community-Notebooks)!- **Title:** How to build an RNA-seq logistic regression classifier with BigQuery ML- **Author:** <NAME>- **Created:** 2021-07-19- **Purpose:** Demonstrate use of BigQuery ML to predict a cancer endpoint using gene expression data.- **URL:** https://github.com/isb-cgc/Community-Notebooks/blob/master/MachineLearning/How_to_build_an_RNAseq_logistic_regression_classifier_with_BigQuery_ML.ipynb- **Note:** This example is based on the work published by [Bosquet et al.](https://molecular-cancer.biomedcentral.com/articles/10.1186/s12943-016-0548-9)This notebook builds upon the [scikit-learn notebook](https://github.com/isb-cgc/Community-Notebooks/blob/master/MachineLearning/How_to_build_an_RNAseq_logistic_regression_classifier.ipynb) and demonstrates how to build a machine learning model using BigQuery ML to predict ovarian cancer treatment outcome. BigQuery is used to create a temporary data table that contains both training and testing data. These datasets are then used to fit and evaluate a Logistic Regression classifier.   Import Dependencies<jupyter_code># GCP libraries
from google.cloud import bigquery
from google.colab import auth<jupyter_output><empty_output><jupyter_text>AuthenticateBefore using BigQuery, we need to get authorization for access to BigQuery and the Google Cloud. For more information see ['Quick Start Guide to ISB-CGC'](https://isb-cancer-genomics-cloud.readthedocs.io/en/latest/sections/HowToGetStartedonISB-CGC.html). Alternative authentication methods can be found [here](https://googleapis.dev/python/google-api-core/latest/auth.html)<jupyter_code># if you're using Google Colab, authenticate to gcloud with the following
auth.authenticate_user()

# alternatively, use the gcloud SDK
#!gcloud auth application-default login<jupyter_output><empty_output><jupyter_text>ParametersCustomize the following parameters based on your notebook, execution environment, or project. BigQuery ML must create and store classification models, so be sure that you have write access to the locations stored in the "bq_dataset" and "bq_project" variables.<jupyter_code># set the google project that will be billed for this notebook's computations
google_project = 'google-project' ## CHANGE ME

# bq project for storing ML model
bq_project = 'bq-project' ## CHANGE ME

# bq dataset for storing ML model
bq_dataset = 'scratch' ## CHANGE ME

# name of temporary table for data
bq_tmp_table = 'tmp_data'

# name of ML model
bq_ml_model = 'tcga_ov_therapy_ml_lr_model'

# in this example, we'll be using the Ovarian cancer TCGA dataset
cancer_type = 'TCGA-OV'

# genes used for prediction model, taken from Bosquet et al.
genes = "'RHOT1','MYO7A','ZBTB10','MATK','ST18','RPS23','GCNT1','DROSHA','NUAK1','CCPG1',\
'PDGFD','KLRAP1','MTAP','RNF13','THBS1','MLX','FAP','TIMP3','PRSS1','SLC7A11',\
'OLFML3','RPS20','MCM5','POLE','STEAP4','LRRC8D','WBP1L','ENTPD5','SYNE1','DPT',\
'COPZ2','TRIO','PDPR'"

# clinical data table
clinical_table = 'isb-cgc-bq.TCGA_versioned.clinical_gdc_2019_06'

# RNA seq data table
rnaseq_table = 'isb-cgc-bq.TCGA.RNAseq_hg38_gdc_current'<jupyter_output><empty_output><jupyter_text>BigQuery ClientCreate the BigQuery client.<jupyter_code># Create a client to access the data within BigQuery
client = bigquery.Client(google_project)<jupyter_output><empty_output><jupyter_text>Create a Table with a Subset of the Gene Expression DataPull RNA-seq gene expression data from the TCGA RNA-seq BigQuery table, join it with clinical labels, and pivot the table so that it can be used with BigQuery ML. In this example, we will label the samples based on therapy outcome. "Complete Remission/Response" will be labeled as "1" while all other therapy outcomes will be labeled as "0". This prepares the data for binary classification. Prediction modeling with RNA-seq data typically requires a feature selection step to reduce the dimensionality of the data before training a classifier. However, to simplify this example, we will use a pre-identified set of 33 genes (Bosquet et al. identified 34 genes, but PRSS2 and its aliases are not available in the hg38 RNA-seq data). Creation of a BQ table with only the data of interest reduces the size of the data passed to BQ ML and can significantly reduce the cost of running BQ ML queries. This query also randomly splits the dataset into "training" and "testing" sets using the "FARM_FINGERPRINT" hash function in BigQuery. "FARM_FINGERPRINT" generates an integer from the input string. More information can be found [here](https://cloud.google.com/bigquery/docs/reference/standard-sql/hash_functions).<jupyter_code>tmp_table_query = client.query(("""
  BEGIN
  CREATE OR REPLACE TABLE `{bq_project}.{bq_dataset}.{bq_tmp_table}` AS
  SELECT * FROM (
    SELECT
      labels.case_barcode as sample,
      labels.data_partition as data_partition,
      labels.response_label AS label,
      ge.gene_name AS gene_name,
      -- Multiple samples may exist per case, take the max value
      MAX(LOG(ge.HTSeq__FPKM_UQ+1)) AS gene_expression
    FROM `{rnaseq_table}` AS ge
    INNER JOIN (
      SELECT
        *
      FROM (
        SELECT
          case_barcode,
          primary_therapy_outcome_success,
          CASE
            -- Complete Reponse    --> label as 1
            -- All other responses --> label as 0
            WHEN primary_therapy_outcome_success = 'Complete Remission/Response' THEN 1
            WHEN (primary_therapy_outcome_success IN (
              'Partial Remission/Response','Progressive Disease','Stable Disease'
            )) THEN 0
          END AS response_label,
          CASE 
            WHEN MOD(ABS(FARM_FINGERPRINT(case_barcode)), 10) < 5 THEN 'training'
            WHEN MOD(ABS(FARM_FINGERPRINT(case_barcode)), 10) >= 5 THEN 'testing'
          END AS data_partition
          FROM `{clinical_table}`
          WHERE
            project_short_name = '{cancer_type}'
            AND primary_therapy_outcome_success IS NOT NULL
      )
    ) labels
    ON labels.case_barcode = ge.case_barcode
    WHERE gene_name IN ({genes})
    GROUP BY sample, label, data_partition, gene_name
  )
  PIVOT (
    MAX(gene_expression) FOR gene_name IN ({genes})
  );
  END;
""").format(
  bq_project=bq_project,
  bq_dataset=bq_dataset,
  bq_tmp_table=bq_tmp_table,
  rnaseq_table=rnaseq_table,
  clinical_table=clinical_table,
  cancer_type=cancer_type,
  genes=genes
)).result()

print(tmp_table_query)<jupyter_output><google.cloud.bigquery.table._EmptyRowIterator object at 0x7f3894001250><jupyter_text>Let's take a look at this subset table. The data has been pivoted such that each of the 33 genes is available as a column that can be "SELECTED" in a query. In addition, the "label" and "data_partition" columns simplify data handling for classifier training and evaluation.<jupyter_code>tmp_table_data = client.query(("""
  SELECT
    * --usually not recommended to use *, but in this case, we want to see all of the 33 genes
  FROM `{bq_project}.{bq_dataset}.{bq_tmp_table}`
""").format(
    bq_project=bq_project,
    bq_dataset=bq_dataset,
    bq_tmp_table=bq_tmp_table
)).result().to_dataframe()

print(tmp_table_data.info())
tmp_table_data<jupyter_output><class 'pandas.core.frame.DataFrame'>
RangeIndex: 264 entries, 0 to 263
Data columns (total 36 columns):
 #   Column          Non-Null Count  Dtype  
---  ------          --------------  -----  
 0   sample          264 non-null    object 
 1   data_partition  264 non-null    object 
 2   label           264 non-null    int64  
 3   RHOT1           264 non-null    float64
 4   MYO7A           264 non-null    float64
 5   ZBTB10          264 non-null    float64
 6   MATK            264 non-null    float64
 7   ST18            264 non-null    float64
 8   RPS23           264 non-null    float64
 9   GCNT1           264 non-null    float64
 10  DROSHA          264 non-null    float64
 11  NUAK1           264 non-null    float64
 12  CCPG1           264 non-null    float64
 13  PDGFD           264 non-null    float64
 14  KLRAP1          264 non-null    float64
 15  MTAP            264 non-null    float64
 16  RNF13           264 non-null    float64
 17  THBS1           264 non-null    flo[...]<jupyter_text>Train the Machine Learning ModelNow we can train a classifier using BigQuery ML with the data stored in the subset table. This model will be stored in the location specified by the "bq_ml_model" variable, and can be reused to predict samples in the future.We pass three options to the BQ ML model: model_type, auto_class_weights, and input_label_cols. Model_type specifies the classifier model type. In this case, we use "LOGISTIC_REG" to train a logistic regression classifier. Other classifier options are documented [here](https://cloud.google.com/bigquery-ml/docs/reference/standard-sql/bigqueryml-syntax-create). Auto_class_weights indicates whether samples should be weighted to balance the classes. For example, if the dataset happens to have more samples labeled as "Complete Response", those samples would be less weighted to ensure that the model is not biased towards predicting those samples. Input_label_cols tells BigQuery that the "label" column should be used to determine each sample's label. **Warning**: BigQuery ML models can be very time-consuming and expensive to train. Please check your data size before running BigQuery ML commands. Information about BigQuery ML costs can be found [here](https://cloud.google.com/bigquery-ml/pricing).<jupyter_code># create ML model using BigQuery
ml_model_query = client.query(("""
  CREATE OR REPLACE MODEL `{bq_project}.{bq_dataset}.{bq_ml_model}`
  OPTIONS
    (
      model_type='LOGISTIC_REG',
      auto_class_weights=TRUE,
      input_label_cols=['label']
    ) AS
  SELECT * EXCEPT(sample, data_partition)  -- when training, we only the labels and feature columns
  FROM `{bq_project}.{bq_dataset}.{bq_tmp_table}`
  WHERE data_partition = 'training' -- using training data only
""").format(
  bq_project=bq_project,
  bq_dataset=bq_dataset,
  bq_ml_model=bq_ml_model,
  bq_tmp_table=bq_tmp_table
)).result()
print(ml_model_query)

# now get the model metadata
ml_model = client.get_model('{}.{}.{}'.format(bq_project, bq_dataset, bq_ml_model))
print(ml_model)<jupyter_output><google.cloud.bigquery.table._EmptyRowIterator object at 0x7f3893663810>
Model(reference=ModelReference(project='isb-project-zero', dataset_id='jhp_scratch', project_id='tcga_ov_therapy_ml_lr_model'))<jupyter_text>Evaluate the Machine Learning ModelOnce the model has been trained and stored, we can evaluate the model's performance using the "testing" dataset from our subset table. Evaluating a BQ ML model is generally less expensive than training. Use the following query to evaluate the BQ ML model. Note that we're using the "data_partition = 'testing'" clause to ensure that we're only evaluating the model with test samples from the subset table.  BigQuery's ML.EVALUATE function returns several performance metrics: precision, recall, accuracy, f1_score, log_loss, and roc_auc. More details about these performance metrics are available from [Google's ML Crash Course](https://developers.google.com/machine-learning/crash-course/classification/video-lecture). Specific topics can be found at the following URLs: [precision and recall](https://developers.google.com/machine-learning/crash-course/classification/precision-and-recall), [accuracy](https://developers.google.com/machine-learning/crash-course/classification/accuracy), [ROC and AUC](https://developers.google.com/machine-learning/crash-course/classification/roc-and-auc).<jupyter_code>ml_eval = client.query(("""
SELECT * FROM ML.EVALUATE (MODEL `{bq_project}.{bq_dataset}.{bq_ml_model}`, 
  (
    SELECT * EXCEPT(sample, data_partition)
    FROM `{bq_project}.{bq_dataset}.{bq_tmp_table}`
    WHERE data_partition = 'testing'
  )
)
""").format(
  bq_project=bq_project,
  bq_dataset=bq_dataset,
  bq_ml_model=bq_ml_model,
  bq_tmp_table=bq_tmp_table
)).result().to_dataframe()
# Display the table of evaluation results
ml_eval<jupyter_output><empty_output><jupyter_text>Predict Outcome for One or More SamplesML.EVALUATE evaluates a model's performance, but does not produce actual predictions for each sample. In order to do that, we need to use the ML.PREDICT function. The syntax is similar to that of the ML.EVALUATE function and returns "label", "predicted_label", "predicted_label_probs", and all feature columns. Since the feature columns are unchanged from the input dataset, we select only the original label, predicted label, and probabilities for each sample. Note that the input dataset can include one or more samples, and must include the same set of features as the training dataset.<jupyter_code>ml_predict = client.query(("""
SELECT
  label,
  predicted_label,
  predicted_label_probs
FROM ML.PREDICT (MODEL `{bq_project}.{bq_dataset}.{bq_ml_model}`, 
  (
    SELECT * EXCEPT(sample, data_partition)
    FROM `{bq_project}.{bq_dataset}.{bq_tmp_table}`
    WHERE data_partition = 'testing' -- Use the testing dataset
  )
)
""").format(
  bq_project=bq_project,
  bq_dataset=bq_dataset,
  bq_ml_model=bq_ml_model,
  bq_tmp_table=bq_tmp_table
)).result().to_dataframe()
# Display the table of prediction results
ml_predict
# Calculate the accuracy of prediction, which should match the result of ML.EVALUATE
accuracy = 1-sum(abs(ml_predict['label']-ml_predict['predicted_label']))/len(ml_predict)
print('Accuracy: ', accuracy)<jupyter_output>Accuracy:  0.6230769230769231
{'fraction_non_alphanumeric': 0.08723559445660102, 'fraction_numerical': 0.019401896425966448, 'mean_word_length': 4.324660194174757, 'pattern_counts': {'":': 0, '<': 38, '<?xml version=': 0, '>': 40, 'https://': 13, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '176', 'license': 'Apache-2.0', 'n_tokens_mistral': 4444, 'n_tokens_neox': 4102, 'n_words': 1350, 'path': 'MachineLearning/How_to_build_an_RNAseq_logistic_regression_classifier_with_BigQuery_ML.ipynb', 'repo_name': 'rpatil524/Community-Notebooks', 'size': '13713'}