<jupyter_start><jupyter_text>Copyright 2019 Google LLCLicensed under the Apache License, Version 2.0 (the "License");you may not use this file except in compliance with the License.You may obtain a copy of the License at    https://www.apache.org/licenses/LICENSE-2.0Unless required by applicable law or agreed to in writing, softwaredistributed under the License is distributed on an "AS IS" BASIS,WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.See the License for the specific language governing permissions andlimitations under the License.<jupyter_code>import os
import pandas as pd
import numpy as np
import scipy.stats
import sklearn.linear_model
# Replace this with the location of the downloaded data and models
datadir = '/path/to/datafiles'<jupyter_output><empty_output><jupyter_text>**Common functions**  **Auxiliary function for design matrices.**<jupyter_code>def construct_difference_operators(times, only_x0=False):
  """Constuct matrix operators for differencing.

  The times arg is the concatention of strictly increasing sequences, each
  starting from zero.  The output operators are all len(times) X len(times)
  matrices that operate on vectors of gene expression values that are len(times)
  long.

  Args:
    times: np.array of time values
    only_x0: if true, return the x0_operator, None, None.

  Returns:
    x0: operator that returns the time zero point for the timecourse
    diff: operator that returns the finite-difference for the timescourse
    diff_inv: operator that returns the timecourse from the finite differences

  Raises:
    ValueError:  First time point isn't zero or timeseries is garbled.
  """
  num_times = len(times)
  # Tag the zeros
  zero_f = times == 0.0
  if not zero_f[0]:
    # First time point is not zero, something is wrong.
    raise ValueError('First time point is non-zero, exiting.')
  else:
    old_t = -1
    for t in times:
      if t != 0 and t <= old_t:
        raise ValueError('Timeseries not sorted properly.')
      old_t = t

  # The x0 operator acting on an expression vector returns a vector of the same
  # length containing the time zero point replicated for the full timecourse,
  # for each timecourse concatenated.
  x0 = np.zeros((num_times, num_times))
  for i, z in enumerate(zero_f):
    if z:
      col = i
    x0[i, col] = 1.0
  if only_x0:
    return x0, None, None

  # Shift times backward, add a zero at the end
  shifted_times = np.append(times[1:], 0.0)
  # Time deltas
  ddt = 0.5 / (shifted_times - times)
  # Tag the endpoints
  last_f = ddt < 0.0
  # The negative values correspond to the "boundaries" between experiments.
  # Set these to zero.
  ddt[last_f] = 0.0
  # Prepend with a zero so that the first experiment can see "backward" to
  # a previous dummy experiment.
  ddt = np.append(0.0, ddt)
  # This is the difference operator
  diff = np.zeros((num_times, num_times))
  for i in range(num_times):
    # Typical case along the diagonal.
    # If the ddt values are equal, this is 0.
    # This corresponds to the usual symmetric difference.
    diff[i, i] = ddt[i] - ddt[i+1]
    # Backward derivative.
    if i - 1 >= 0:
      diff[i, i-1] = - ddt[i]
    # Forward derivative.
    if i + 1 < num_times:
      diff[i, i+1] = ddt[i+1]
  # This is the inverse of the difference operator
  diff_inv = np.zeros((num_times, num_times))
  # It is not invertible.  However, removing the first column
  # and last row from each timecourse makes it invertible.
  # We do that because we don't care about the zero point of the sequence,
  # thus filtering columns for non-zero time points, and we don't care about
  # the final derivative, thus filtering rows for "not last".
  d_inv = np.linalg.inv(diff[np.ix_(~last_f, ~zero_f)])
  # Fill in the values.
  # Note that in constructing the data from the derivatives,
  # we do not try to reconstruct x0, thus the row filter on non-zero times
  # and we do not use the final derivative of the timecourse due to linear
  # dependence, thus the column filter on "not last".
  # diff_inv . diff . X == X - X0
  diff_inv[np.ix_(~zero_f, ~last_f)] = d_inv

  # Return the three operators.
  return x0, diff, diff_inv
def add_gene_block(df_base, df_add, block_name, genes):
  # Column-wise bind df_base (dataframe) and df_add (may be np.array).
  # Make a multindex of (blockname, gene) where gene is taken from genes.
  df_add.columns = pd.MultiIndex.from_tuples(
      [(block_name, x) for x in genes])
  return pd.concat([df_base, df_add], axis=1)
TIMECOURSE_IDX = ['TF', 'strain', 'date', 'restriction', 'mechanism']
CHIP_IDX = TIMECOURSE_IDX + ['time']
TIMESERIES_IDX = ['GeneName', 'SystematicName'] + CHIP_IDX
DESIGN_IDX = CHIP_IDX + ['seq', 'rseq']<jupyter_output><empty_output><jupyter_text>**Auxiliary functions for noise model**<jupyter_code>def normalize_to_t0(df_dm, x0_op):
  """Constructs variables relative to t0.

  Args:
    df_dm: dataframe of the design_matrix
    x0_op: operator to produce the t=0 point for the timeseries

  Returns:
    design matrix dataframe with ratio block.
  """
  genes = df_dm['r_g_median'].columns

  # The basic stat used is the median.
  # Here we compare with the zero time
  df_dm = add_gene_block(
      df_dm,
      np.log(df_dm['r_g_median'] /
             np.matmul(x0_op, df_dm['r_g_median']).values),
      'log_ratio',
      genes)

  return df_dm
def repair_crosstalk_outliers(df, x0_op):
  """Corrects for timeseries with severe red - green crosstalk.

  Args:
    df: dataframe of the design_matrix
    x0_op: operator to produce the t=0 point for the timeseries

  Returns:
    design matrix dataframe with corrected ratio.
  """
  # Pull out the zero points of the green channel, convert back to DataFrame.
  green0 = pd.DataFrame(
      np.matmul(x0_op, df['green_median']).values,
      index=df.index,
      columns=df['green_median'].columns)

  # Removing the t0 point, these are the 30% quantiles of each
  # timecourse + gene pair.  This covers 2/7 of the remaining points.
  # The stack() turns the gene columns into rows.
  gratio = ((df['green_median'] / green0).query('time > 0').
            groupby(level=['TF', 'strain']).quantile(0.3).stack())
  lratio = (df['log_ratio'].query('time > 0').
            groupby(level=['TF', 'strain']).quantile(0.3).stack())

  # Flag the genes that both have green going way up (8x) and also the
  # uncorrected ratio going up 2x beyond that as "crosstalk" genes.
  # When iterated on, this index is (TF, strain, gene) tuples.
  crosstalk_genes = gratio[(gratio > 8.0) & (lratio > np.log(2.0))].index
  for tf, strain, gene in crosstalk_genes:
    print(
        'Repairing %s, %s timecourse for gene %s.' % (tf, strain, gene))
    tc_f = ((df.index.get_level_values('TF') == tf) &
            (df.index.get_level_values('strain') == strain))
    # This adds log(green / green0), which yields
    # log_ratio = log(red) - log(green0)
    # (the green0 is propagated across the timecourse).
    df.loc[tc_f, ('log_ratio', gene)] += np.log(
        df.loc[tc_f, ('green_median', gene)] / green0.loc[tc_f, gene])

  return df
# The typical timecourse was 0, 5, 10, 15, 20, 30, 45, 90
TIME_POINTS = (0, 6, 12, 17.5, 25, 35, 85, 1000)
TIME_BINS = zip(TIME_POINTS[:-1], TIME_POINTS[1:])
def clean_time_trends(df_dm):
  """Corrects for median timeseries behavior.

  Args:
    df_dm: dataframe of the design_matrix

  Returns:
    design matrix dataframe with new blocks.
  """
  genes = df_dm['log_ratio'].columns

  dfr = df_dm['log_ratio']
  dfrc = dfr.copy()
  times = dfr.index.get_level_values('time')
  gev = dfr.index.get_level_values('mechanism') == 'GEV'
  for mn, mx in TIME_BINS:
    for g in (gev, ~gev):
      filt = (times > mn) & (times <= mx) & g
      dfrc[filt] -= dfr[filt].median()

  df_dm = add_gene_block(
      df_dm,
      dfrc,
      'log_cleaned_ratio',
      genes)

  return df_dm
def build_cleaned_table(df_dm):
  """Augment design table and write to disk.

  Args:
    df_dm: design table dataframe

  Returns:
    design table dataframe with normalized, cleaned, and differentiated columns
  """
  x0_op, _, _ = construct_difference_operators(
      np.array(df_dm.index.get_level_values('time')), only_x0=True)

  df_dm = normalize_to_t0(df_dm, x0_op)

  df_dm = repair_crosstalk_outliers(df_dm, x0_op)

  df_dm = clean_time_trends(df_dm)

  df_dm = construct_noise_model(df_dm)

  return df_dm
# For the normal distribution, the central quintile is +/- 0.253, approximately.
# We scale to infer the standard deviation.
CENTRAL_QUINTILE_SCALE = 1.0 / (scipy.stats.norm.ppf(0.6) -
                                scipy.stats.norm.ppf(0.4))

def central_quintile_error(dfc, axis=0):
  # Use central quintile to define the noise level.
  return (dfc.quantile(0.6, axis=axis) -
          dfc.quantile(0.4, axis=axis)) * CENTRAL_QUINTILE_SCALE


def log_cleaned_gene_error(df):
  return central_quintile_error(df.query('time > 0')['log_cleaned_ratio'])


def chip_residual_error(df):
  return central_quintile_error(df[df != 0], axis=1)
def apply_soft_thresholding(df, noise):
  return np.median(
      np.array([np.array(df + noise),
                np.array(df - noise),
                np.zeros(df.shape)]),
      axis=0)


def apply_hard_thresholding(df, noise):
  hth = df.copy()
  hth[np.abs(df) < noise] = 0.0
  return hth
def count_timecourse_detections(df):
  # Select active genes for each timecourse.
  # A gene is active in a timecourse if it ever passes hard thresholding.
  # Drop the unused names from the initial index.
  drop_names = [x for x in df.index.names if x not in TIMECOURSE_IDX]

  return (
      # Significant data points (bool)
      ((df != 0).
       # Expanded index
       reset_index().
       # Group by timecourses
       groupby(TIMECOURSE_IDX).
       # Sum just counts significant points
       sum()).
      # Drop unwanted terms in expanded index
      drop(drop_names, axis=1, errors='ignore'))


def flag_non_zero_timecourses(df):
  return count_timecourse_detections(df) != 0


def count_consecutive_timecourse_detections(df):
  # Check for consecutive points passing thresholding.
  # Shift forward
  shifted = np.zeros_like(df)
  shifted[:-1] = (df[1:] != 0)
  # consec == 0 for first and last points, automatically.
  consec = (df != 0) * shifted
  return count_timecourse_detections(consec)


def strip_timecourses(df, sig_tc, significant=True):
  """Strip timecourses according to a set of significant ones.

  Args:
    df: design matrix dataframe to be stripped.
    sig_tc: boolean design with timecourse index and gene columns.
    significant: if true, strip the significant timecourses,
                 otherwise, strip the insignificant ones.

  Returns:
    dataframe of stripped design matrix
  """
  dfn = df.copy()
  for idx, sig in sig_tc.iterrows():
    # This builds the query of index_name==index_value, joined with ' and '.
    query = ' and '.join(['%s=="%s"' % t for t in zip(TIMECOURSE_IDX, idx)])
    filt = dfn.query(query)
    dfn.loc[filt.index] *= np.array(~sig if significant else sig, dtype=float)
  return dfn
MINIMUM_THRESHOLDING = 0.1
CHIP_WEIGHT = 0.2
def construct_noise_model(df_dm):
  """Builds a gene X chip noise model.

  Args:
    df_dm: dataframe of the design_matrix

  Returns:
    design matrix dataframe with new blocks.
  """
  genes = df_dm['log_cleaned_ratio'].columns
  gene_error = log_cleaned_gene_error(df_dm)

  # Thresholding factor.
  threshold_factor = np.sqrt(2.0 * np.log(
      len(df_dm.query('time > 0'))))
  threshold_error = np.maximum(
      gene_error * threshold_factor,
      MINIMUM_THRESHOLDING)

  hth = apply_hard_thresholding(
      df_dm['log_cleaned_ratio'], threshold_error)
  # For each timecourse, zero out the significant genes.
  sig_tc = flag_non_zero_timecourses(hth)
  dfn = strip_timecourses(
      df_dm['log_cleaned_ratio'], sig_tc)
  # Compute chip-level error via central quintile.
  chip_error = chip_residual_error(dfn)
  filt = chip_error.index.get_level_values('time') == 0.0
  # Fill in 0.0 for the t=0 points. t=0 will thus have gene-level noise
  chip_error[filt] = 0.0

  noise_df = pd.DataFrame(
      np.sqrt(np.ones(dfn.shape) * gene_error.values.reshape(1, -1)**2 +
              CHIP_WEIGHT *
              np.ones(dfn.shape) * chip_error.values.reshape(-1, 1)**2),
      index=chip_error.index,
      columns=gene_error.index)

  # Normalize the non-zero points properly.
  nonzero = noise_df.query('seq != 0')
  noise_df.loc[nonzero.index] /= np.sqrt(1.0 + CHIP_WEIGHT)

  # Output the noise model.
  df_dm = add_gene_block(
      df_dm,
      noise_df,
      'log_noise_model',
      genes)

  return df_dm
LOG_SUFFIXES = ['noise_model',
                'ratio',
                'cleaned_ratio',
                'cleaned_ratio_zth',
                'cleaned_ratio_hth',
                'cleaned_ratio_th',
                'cleaned_ratio_zth2d',
                'cleaned_ratio_hth2d',
                'cleaned_ratio_th2d',
                'cleaned_ratio_zth2dfilt',
                'cleaned_ratio_hth2dfilt',
                'cleaned_ratio_th2dfilt']<jupyter_output><empty_output><jupyter_text>**Process "raw" data**  **Read in the raw data and verify that the red to green ratio column is correctly computed.**<jupyter_code># Read in the raw data.

with open(os.path.join(datadir, 'yeast_raw_data_table_20180826.tsv')) as f:
  df_raw = pd.read_csv(f, sep='\t')
# The counts are red(mutant) and green(control).  r_g_ratio is the ratio of these, clipped at values of 2.0.

d = np.maximum(df_raw['rProcessedSignal'], 2.0) / np.maximum(df_raw['gProcessedSignal'], 2.0) - df_raw['r_g_ratio']
min(d), max(d)<jupyter_output><empty_output><jupyter_text>**Aggregate across spots on the microarrays.**<jupyter_code># Aggregate across the multiple spots (usually 2) on the microarrays.

df_ts = df_raw.pivot_table(index=TIMESERIES_IDX,
                           values=['r_g_ratio',
                                   'rProcessedSignal',
                                   'gProcessedSignal'],
                           aggfunc=[len, np.median, np.min, np.max, np.std])
# Save memory.
del df_raw<jupyter_output><empty_output><jupyter_text>**Trim the most egregious outliers.**<jupyter_code>MAX_SPOT_RATIO = 4.0
JUMP_THRESHOLD = 4.0
def trim_timeseries_outliers(df_ts):
  """Replace outliers in timeseries with bracketing averages.

  Args:
    df_ts: dataframe with gene X experiment X time per row

  Returns:
    dataframe where timeseries outliers have been replaced by
    bracketing averages.
  """
  # Shift times down by one to get "next" time.
  next_t = pd.Series(
      np.append(np.array(df_ts.iloc[1:].index.get_level_values('time')), 0),
      index=df_ts.index)
  # Values shifted forward and backward by one in time.
  next_val = pd.Series(
      np.append(np.array(df_ts.iloc[1:][('median', 'r_g_ratio')]), 0),
      index=df_ts.index)
  prev_val = pd.Series(
      np.append(0, np.array(df_ts.iloc[:-1][('median', 'r_g_ratio')])),
      index=df_ts.index)
  orig_val = df_ts[('median', 'r_g_ratio')]
  geometric_mean = np.sqrt(prev_val * next_val)

  # filters for first and last in a timecourse.
  first_t = df_ts.index.get_level_values('time') == 0
  last_t = df_ts.index.get_level_values('time') > next_t

  # Filter for genes on a chip where the spots disagree.
  filt = (df_ts[('amax', 'r_g_ratio')] /
          np.maximum(df_ts[('amin', 'r_g_ratio')], 0.001) >
          MAX_SPOT_RATIO)

  # Replace the cross-spot medians for the inconsistent genes.
  m_filt = filt & ~first_t & ~last_t
  df_ts.loc[m_filt, ('median', 'r_g_ratio')] = geometric_mean[m_filt]

  f_filt = filt & first_t
  df_ts.loc[f_filt, ('median', 'r_g_ratio')] = next_val[f_filt]

  l_filt = filt & last_t
  df_ts.loc[l_filt, ('median', 'r_g_ratio')] = prev_val[l_filt]

  print('%d / %d spot-discrepant outliers trimmed. ' % (
      sum(filt), len(filt)))

  # Check for outliers from the time-trend in r_g_ratio
  m_filt = ~first_t & ~last_t
  ratio1 = orig_val / prev_val
  ratio2 = next_val / orig_val
  # outlier is a non-endpoint location where there was movement
  # up then down or down then up, each over threshold.
  outlier_filt = m_filt & (
      ((ratio1 > JUMP_THRESHOLD) & (ratio2 < 1. / JUMP_THRESHOLD) |
       (ratio1 < 1. / JUMP_THRESHOLD) & (ratio2 > JUMP_THRESHOLD)))
  print('%d / %d time-series trend outliers trimmed.' % (
               np.sum(outlier_filt), len(outlier_filt)))
  # Replace them with the geometric mean of their neighbors.
  df_ts.loc[outlier_filt, ('median', 'r_g_ratio')] = geometric_mean[
      outlier_filt]

  return df_ts

df_ts = trim_timeseries_outliers(df_ts)<jupyter_output><empty_output><jupyter_text>**Build a design table with rows as experiment time points and columns as transcripts.**<jupyter_code>def design_table(df_tsm):
  """Generates the design matrix for linear modeling.

  Args:
    df_tsm: dataframe of the full dataset

  Returns:
    dataframe containing design matrix
  """
  df_tsmf = df_tsm[[('median', 'r_g_ratio'),
                    ('median', 'rProcessedSignal'),
                    ('median', 'gProcessedSignal')]]
  df_tsmf.columns = ['r_g_median', 'red_median', 'green_median']

  df_dm = pd.pivot_table(
      df_tsmf.reset_index(),
      index=CHIP_IDX,
      columns=['GeneName'])

  return df_dm
def timeseries_sequence(times):
  """Compute the timecourse sequence number of the time points, and reverse."""
  diffs = np.append(times[1:], 0) - times
  i = 0
  seq = []
  for x in diffs:
    seq.append(i)
    if x <= 0:
      i = 0
    else:
      i += 1

  # Now in reverse.
  rtimes = times[::-1]
  rdiffs = np.append(0, rtimes[:-1]) - rtimes
  i = 0
  rseq = []
  for x in rdiffs:
    if x <= 0:
      i = 0
    else:
      i += 1
    rseq.append(i)

  # Don't forget to re-reverse rseq
  return seq, rseq[::-1]

df_dm = design_table(df_ts)

seq, rseq = timeseries_sequence(
      np.array(df_dm.index.get_level_values('time')))

df_dm['seq'] = seq
df_dm['rseq'] = rseq

df_dm = df_dm.reset_index().set_index(DESIGN_IDX).sort_index()
df_dm.shape
# There a few mismatches between the transcripts on GEV vs. ZEV microarrays.
# Reduce to the common set.

df_dm.dropna(axis=1, inplace=True)
df_dm.shape
# Save memory.
del df_ts<jupyter_output><empty_output><jupyter_text>**Build the noise model and "cleaned" dataset.**<jupyter_code># Data cleaning and noise model depend on the full dataset.
# To reproduce the first round of results, we need to filter out
# the rounds of validation experiments.

validation_f = pd.to_datetime(df_dm.index.get_level_values('date')) > '2017-07-01'
# This will recapitulate the original training set,
# before validation experiments were performed.
# The outputs of cleaning and noise model are slightly different.
#
# df_dm = build_cleaned_table(df_dm[~validation_f])

# This recapitulates the full dataset.
df_dm = build_cleaned_table(df_dm)
df_dm.shape<jupyter_output><empty_output><jupyter_text>**Apply noise model to generate "thresholded" dataset.**<jupyter_code>def apply_all_thresholding(df_dm, min_thresholding=0.1):
  """Adds thresholded data blocks to design matrix.

  Args:
    df_dm: dataframe of the design_matrix
    min_thresholding: floor for amount of thresholding

  Returns:
    design matrix dataframe with new blocks.
  """
  genes = df_dm['log_ratio'].columns
  # Soft thresholding.
  # Consider the max of N draws from a normal distribution.
  n_eff = len(df_dm.query('time > 0'))
  soft_threshold_factor = np.sqrt(2.0 * np.log(n_eff))

  print(
      'Thresholding at %.3f sigma for %d samples. ' % (
      soft_threshold_factor,
      n_eff))

  gene_error = df_dm['log_noise_model'].query('time==0').median()

  threshold_error = np.maximum(
      gene_error * soft_threshold_factor,
      min_thresholding)

  df_dm = apply_hard_and_soft_thresholding(
      df_dm, threshold_error, genes, 'th')

  threshold_error_2d = np.maximum(
      df_dm['log_noise_model'] * soft_threshold_factor,
      min_thresholding)

  df_dm = apply_hard_and_soft_thresholding(
      df_dm, threshold_error_2d, genes, 'th2d')

  df_dm = apply_hard_and_soft_thresholding(
      df_dm, threshold_error_2d, genes, 'th2dfilt', aggressive=True)

  return df_dm
def apply_hard_and_soft_thresholding(df_dm, noise, genes, suffix,
                                     aggressive=False):
  """Adds thresholded blocks to the design matrix.

  Args:
    df_dm: dataframe of the design_matrix
    noise: compatible dataframe (perhaps broadcast) with noise model
    genes: list of gene names
    suffix: suffix for the names of the new blocks
    aggressive: if true, apply simple timecourse-level filter as well

  Returns:
    design matrix dataframe with new blocks.
  """
  # Here we hard threshold the values
  ratio_hth = apply_hard_thresholding(df_dm['log_cleaned_ratio'], noise)

  tc_counts = count_timecourse_detections(ratio_hth)
  tc_consec_counts = count_consecutive_timecourse_detections(ratio_hth)
  last_counts = count_timecourse_detections(ratio_hth.query('rseq==0'))

  print('%d timecourses with signal (%s)' % (
               (tc_counts > 0).sum().sum(),
               suffix))
  print('%d timecourses with one detection (%s)' % (
               (tc_counts == 1).sum().sum(),
               suffix))
  print('%d timecourses with only final detection (%s)' % (
               ((tc_counts == 1) & last_counts).sum().sum(),
               suffix))
  print('%d timecourses with no detection gaps (%s)' % (
               (tc_consec_counts == (tc_counts - 1)).sum().sum(),
               suffix))

  if aggressive:
    # Non-final single detection is bad.
    notlast = (tc_counts == 1) & (last_counts == 0)
    # Non-consecutive detection is bad.
    nonconsec = tc_consec_counts < (tc_counts - 1)

    bad_tc = notlast | nonconsec

    good_tc = (tc_counts > 0) & ~bad_tc
    print('%d timecourses with no detection gaps or singleton not last (%s)' % (
              good_tc.sum().sum(),
              suffix))

    stripped_tc = strip_timecourses(df_dm['log_cleaned_ratio'],
                                    good_tc,
                                    significant=False)
    # We need the hard thresholded data without the bad timecourses.
    ratio_hth = apply_hard_thresholding(stripped_tc, noise)

  else:
    stripped_tc = strip_timecourses(df_dm['log_cleaned_ratio'],
                                    tc_counts > 0,
                                    significant=False)

  df_dm = add_gene_block(
      df_dm,
      stripped_tc.copy(),
      'log_cleaned_ratio_z' + suffix,
      genes)

  df_dm = add_gene_block(
      df_dm,
      ratio_hth,
      'log_cleaned_ratio_h' + suffix,
      genes)

  # Here we soft threshold the values
  ratio_th = apply_soft_thresholding(stripped_tc, noise)

  df_dm = add_gene_block(
      df_dm,
      pd.DataFrame(ratio_th, index=df_dm.index, columns=genes),
      'log_cleaned_ratio_' + suffix,
      genes)

  return df_dm
df_dm = apply_all_thresholding(
        df_dm,
        min_thresholding=MINIMUM_THRESHOLDING)
df_dm['log_cleaned_ratio_th'].shape<jupyter_output><empty_output><jupyter_text>**Generate the data table in log2 units.**<jupyter_code>def biologist_table(df):
  """Output a table in biologist format - stacked and log base 2."""
  base_columns = (DESIGN_IDX +
                  ['GeneName',
                   'green_median',
                   'red_median',
                   'r_g_median'])

  dfs = df.stack().reset_index()
  log_columns = ['log_' + x for x in LOG_SUFFIXES]
  dfs[log_columns] /= np.log(2.0)

  log2_columns = ['log2_' + x for x in LOG_SUFFIXES]
  dfb = dfs[base_columns + log_columns]
  dfb.columns = [base_columns + log2_columns]
  return dfb
df_bio = biologist_table(df_dm)
df_bio.shape<jupyter_output><empty_output><jupyter_text>**Read processed data**  **Read in the processed table.**<jupyter_code># Read in the final data.

with open(os.path.join(datadir, 'yeast_data_table_20180826.tsv')) as f:
  df_bio_read = pd.read_csv(f, sep='\t')
df_bio_read.shape<jupyter_output><empty_output><jupyter_text>**Verify the processing procedure (only if *Process "raw" data* was run)**<jupyter_code>v = np.array(df_bio['log2_cleaned_ratio']).flatten() - df_bio_read['log2_cleaned_ratio']
min(v), max(v)
v = np.array(df_bio['log2_cleaned_ratio_th2d']).flatten() - df_bio_read['log2_cleaned_ratio_th2d']
min(v), max(v)
del df_bio<jupyter_output><empty_output><jupyter_text>**Convert processed table back to design matrix format.**<jupyter_code>def design_from_biologist(df):
  df.set_index(DESIGN_IDX+['GeneName'], inplace=True)
  log2_columns = ['log2_' + x for x in LOG_SUFFIXES]
  df[log2_columns] *= np.log(2.0)
  cols = [x.replace('log2', 'log') for x in df.columns]
  df.columns = cols
  return df.unstack()
df_dm_read = design_from_biologist(df_bio_read)
df_dm = df_dm_read
del df_bio_read<jupyter_output><empty_output><jupyter_text>**Add expression level variables (not logarithmic).**<jupyter_code>def add_linear_column_blocks(df):
  genes = df['r_g_median'].columns
  for block in df.columns.levels[0]:
    if block.startswith('log_'):
      df = add_gene_block(df, np.exp(df[block]), block[4:], genes)
  return df
df_dm = add_linear_column_blocks(df_dm)<jupyter_output><empty_output><jupyter_text>**Linear modeling**  **Construct the design for modeling.**<jupyter_code>def construct_coefficient_labels(df, model_quadratic=True,
                                 corrected_intercept=False):
  """Compute the column names for the design matrix."""
  # These are the names of the genes.
  gene_labels = df['r_g_median'].columns

  # These are the names of the coefficients.
  # The quadratic terms need to be added.
  coef_labels = []
  coef_labels.extend(gene_labels)
  if model_quadratic:
    coef_labels.extend(['quad_' + x for x in gene_labels])
  if corrected_intercept:
    coef_labels.append('corrected_intercept')
  return pd.Series(coef_labels)
def construct_full_design(df,
                          in_sample_f,
                          out_of_sample_f,
                          gene,
                          delta_model=False,
                          thresholded='',
                          model_quadratic_terms=True,
                          log_weights=False,
                          inverse_design=False,
                          corrected_intercept=False,
                          x0_op=None,
                          diff_op=None,
                          inv_diff_op=None):
  """Construct the design matrix and dependent variable.

  For a given gene, the self promoted experiments are filtered out.  Then,
  the design matrix for the quadratic model is constructed.

  Args:
    df: data frame with time points as rows, genes as columns.
    in_sample_f: Build model on this filter
    out_of_sample_f: Compute out of sample deviance on this filter
    gene: gene to build model for.
    delta_model: whether to subtract assumed steady state solution.
    thresholded: suffix tag for thresholded data.
    model_quadratic_terms: whether to include quadratic terms
    log_weights: whether to apply weights to make log model
    inverse_design: (boolean) whether the design matrix is inverse differenced.
    corrected_intercept: (boolean) Fit intercept with "integrate" and "log",
    x0_op: matrix to produce vector of time zero points for each timecourse
    diff_op: matrix for differencing operator
    inv_diff_op: matrix for inverse differencing operator
  Returns:
    design matrix, dependent variable, time=0 values,
    weights, row_index, col_index, in_sample_f, out_of_sample_f
  """
  notlast_f = df.index.get_level_values('rseq') != 0
  notzero_f = df.index.get_level_values('time') != 0
  row_index = df.index

  col_index = pd.Series(
      construct_coefficient_labels(
          df, model_quadratic=model_quadratic_terms,
          corrected_intercept=corrected_intercept))
  print('Promoter experiment filtered, shape = %s' % (df.shape,))

  # Select column sets according to input booleans.
  val = 'cleaned_ratio' + thresholded
  yval = val
  # Construct the design matrix d_m.
  #
  # df has multi-indexed columns, where the final piece is "gene",
  # e.g. df['ratio'] has one column per gene.
  #
  # Matrices, with rows (strain X time) and columns for genes:
  #
  # u_m is the matrix of r/g values, divided by the t0 ratio
  # uu is a matrix of the same shape with the quadratic terms.
  #
  # u is the column of r/g values for the gene in question.
  u_m = np.array(df[val])
  if model_quadratic_terms:
    u = np.array(df[val][gene]).reshape(-1, 1)
    uu = u * u_m
    d_m = np.hstack([u_m, uu])
  else:
    d_m = u_m

  # We now subtract 1.0 from the design matrix, giving it
  # column blocks like (u-1)(uu-1).
  # The intercept is now the deviation from steady state at t=0.
  # For a delta model we assume steady state, with intercept=0.
  d_m -= 1.0

  # We need to add a column for a non-standard intercept.
  # This will be transformed by at least one of the log_weights and
  # inverse_design clauses immediately below.
  if corrected_intercept:
    d_m = np.append(d_m, np.ones((d_m.shape[0], 1)), axis=1)

  # For log models, we now divide by the gene in question.
  if log_weights:
    u = np.array(df[val][gene]).reshape(-1, 1)
    d_m /= u
    yval = 'log_cleaned_ratio' + thresholded

  # For inverse models, we now "integrate" the design matrix.
  if inverse_design:
    d_m = np.dot(inv_diff_op, d_m)

  if inverse_design:
    yfull = df[yval][gene]
    # The inverse difference operator returns 0 for all timeseries 0 points.
    # We thus subtract out the zero values, and filter the 0 rows.
    yzero = yfull - np.dot(x0_op, np.array(yfull))
    y = yzero[notzero_f]
    d_mat = d_m[notzero_f]
    row_index = df.index[notzero_f]
    weights = None
    in_sample_f = in_sample_f[notzero_f]
    out_of_sample_f = out_of_sample_f[notzero_f]
  else:
    # The final derivative is linearly dependent on the previous ones,
    # because there are only N-1 differences.  Thus, we filter on rows that are
    # not the final point in their timecourse.
    y = np.dot(diff_op, df[yval][gene])[notlast_f]
    d_mat = d_m[notlast_f]
    row_index = df.index[notlast_f]
    weights = 1.0 / (diff_op ** 2).sum(axis=1)[notlast_f]
    in_sample_f = in_sample_f[notlast_f]
    out_of_sample_f = out_of_sample_f[notlast_f]

  y_base = np.dot(x0_op, np.array(df['r_g_median'][gene]))

  return (d_mat, y, y_base, weights, row_index, col_index,
          in_sample_f, out_of_sample_f)
# Picking a gene to model - choose one that is a transcption factor
gene = 'FKH1'
promoted = df_dm.index.get_level_values('TF') == gene

# Only hold out the self-experiment
in_sample_f = np.full(len(df_dm.index), True, dtype=bool)
out_of_sample_f = np.full(len(df_dm.index), False, dtype=bool)

in_sample_f = in_sample_f[~promoted]

out_of_sample_original_size_f = pd.Series(
      out_of_sample_f & ~promoted, index=df_dm.index)
out_of_sample_f = out_of_sample_f[~promoted]

df = df_dm[~promoted]

# Build the differencing operators
times = np.array(df.index.get_level_values('time'))
x0_op, diff_op, inv_diff_op = construct_difference_operators(times)
(design_m, y_dep, y_base, weights, row_index,
 col_index, in_sample_f, out_of_sample_f) = construct_full_design(
       df,
       in_sample_f,
       out_of_sample_f,
       gene,
       delta_model=True,
       thresholded='_zth2dfilt',
       model_quadratic_terms=True,
       log_weights=True,
       inverse_design=False,
       corrected_intercept=False,
       x0_op=x0_op,
       diff_op=diff_op,
       inv_diff_op=inv_diff_op)
design_m.shape
y_dep.shape<jupyter_output><empty_output><jupyter_text>**Build a simple linear model (sklearn, not glmnet).**<jupyter_code># We are attempting to match our weighting scheme with GLMnet, imperfectly.
larsmodel = sklearn.linear_model.Lars(fit_intercept=False, normalize=False)
larsmodel.fit(design_m, y_dep)
lambdacol = 40
f = larsmodel.coef_path_[:, lambdacol] != 0
coefs = (larsmodel.coef_path_)[f, lambdacol]
for gene, coef in zip(col_index[f], coefs):
  print('%15s % .4e' % (gene, coef))<jupyter_output><empty_output><jupyter_text>**Prediction model**  **Read in the prediction model.**<jupyter_code># Read in the prediction model.

# Before validation experiments.
coeffs = pd.read_csv(os.path.join(datadir, 'insample_coefs_20170601.csv'))

# After validation experiments.
#coeffs = pd.read_csv(os.path.join(datadir, 'insample_coefs_20180826.csv'))<jupyter_output><empty_output><jupyter_text>**Construct predictions from model.**<jupyter_code>design = df_dm['cleaned_ratio_zth2dfilt']
design.shape
x0op, diff, diff_inv = construct_difference_operators(np.array(design.index.get_level_values('time')))
# Predict for for each gene in turn.
gene = 'YHB1'

# Don't predict the TF in its own experiment!
tf_filt = design.index.get_level_values('TF') != gene
# Data for gene in question
actualg = design[gene]
actualg_diff = np.matmul(diff, np.log(actualg))
# Consider coefficients affecting the gene.
f = coeffs['effect'] == gene
n_c = sum(f)
# 2 extra columns for data, full prediction
marginal_matrix = np.zeros((len(actualg), n_c + 2))
# Fill in the data
marginal_matrix[:, 0] = actualg_diff
# Initialize the full prediction.
prediction = 0. * actualg
# Iterate over the coefficients that affect this gene.
for i, x in enumerate(coeffs[f].iterrows()):
  cause = x[1]['cause']
  coef = x[1]['coef'] * tf_filt
  # The selected model enforces intercept == 0, and is "log".
  # Compute the marginal contribution to (d/dt) ln(gene)
  if cause.startswith('quad'):
    marginal = (design[cause[5:]] - 1.0 / design[gene]) * coef
  else:
    marginal = (design[cause] - 1.0) / design[gene] * coef
  marginal_matrix[:, i + 2] = marginal
  prediction += marginal
# "Integrate" with the diff_inv operator to get ln(gene), then exponentiate
data_pred = np.exp(np.matmul(diff_inv, prediction))
# Fill in the predicted results
marginal_matrix[:, 1] = prediction
# Write out marginals
columns = ['data', 'prediction'] + list(coeffs[f]['cause'])
marginal_diff_df = pd.DataFrame(marginal_matrix, index=design.index, columns=columns)
marginal_df = pd.DataFrame(np.exp(np.matmul(diff_inv, marginal_matrix)), index=design.index, columns=columns)<jupyter_output><empty_output><jupyter_text>**Investigate individual timecourse.**<jupyter_code>aft1 = marginal_df.query('TF=="AFT1"')
hits = np.sum(np.abs(np.log(aft1)) > 0.2) > 0
aft1[aft1.columns[hits]]<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10790740848175437, 'fraction_numerical': 0.010790740848175436, 'mean_word_length': 3.4825120270445975, 'pattern_counts': {'":': 0, '<': 87, '<?xml version=': 0, '>': 96, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '19', 'dir': 'jupyter-structured-clean-dedup', 'id': '92081', 'license': 'Apache-2.0', 'n_tokens_mistral': 12264, 'n_tokens_neox': 11735, 'n_words': 3192, 'path': 'yeast_transcription_network/YeastNetwork.ipynb', 'repo_name': 'ojInc/google-research', 'size': '34474'}