Commit 5509c30d authored by Trent Hauck's avatar Trent Hauck
Browse files

Go with first option.

parent 84fe50e8
Loading
Loading
Loading
Loading
+47 −84
Original line number Diff line number Diff line
@@ -8,7 +8,6 @@ from __future__ import unicode_literals
# TODO(rbharath): Use standard joblib once old-data has been regenerated.
import joblib
from sklearn.externals import joblib as old_joblib
from sklearn.preprocessing import LabelEncoder, OneHotEncoder
import gzip
import json
import pickle
@@ -108,7 +107,6 @@ def load_csv_files(filenames, shard_size=None, verbose=True):
        yield df



def seq_one_hot_encode(sequences, letters='ATCGN'):
  """One hot encodes list of genomic sequences.

@@ -131,38 +129,51 @@ def seq_one_hot_encode(sequences, letters='ATCGN'):
  np.ndarray: Shape (N_sequences, 4, sequence_length, 1).
  """

  sequence_length = len(sequences[0])
  letters_length = len(letters)
  # The label encoder is given characters for ACGTN
  letter_encoder = {l: i for i, l in enumerate(letters)}
  alphabet_length = len(letter_encoder)

  # Peak at the first sequence to get the length of the sequence.
  try:
    first_seq = next(sequences)
    tail_seq = sequences
  except TypeError:
    first_seq = sequences[0]
    tail_seq = sequences[1:]

  sequence_length = len(first_seq)

  seqs = []

  # depends on Python version
  integer_type = np.int32
  seqs.append(
      _seq_to_encoded(first_seq, letter_encoder, alphabet_length,
                      sequence_length))

  # The label encoder is given characters for ACGTN
  letters_array = np.array((letters,))
  label_encoder = LabelEncoder().fit(letters).view(integer_type))
  for other_seq in tail_seq:
    if len(other_seq) != sequence_length:
      raise ValueError

    seqs.append(
        _seq_to_encoded(other_seq, letter_encoder, alphabet_length,
                        sequence_length))

  integer_array = []
  return np.expand_dims(np.array(seqs), -1)

  # TODO(rbharath): Unlike the DRAGONN implementation from which this
  # was ported, I couldn't transform the "ACGT..." strings into
  # integers all at once. Had to do one at a time. Might be worth
  # figuring out what's going on under the hood.

  for sequence in sequences:
    if len(sequence) != sequence_length:
      raise ValueError("All sequences must be of same length")
    integer_seq = label_encoder.transform(
        np.array((sequence,)).view(integer_type))
    integer_array.append(integer_seq)
def _seq_to_encoded(seq, letter_encoder, alphabet_length, sequence_length):
  b = np.zeros((alphabet_length, sequence_length))
  seq_ints = [letter_encoder[s] for s in seq]
  b[seq_ints, np.arange(sequence_length)] = 1

  return b

  integer_array = np.concatenate(integer_array)
  integer_array = integer_array.reshape(len(sequences), sequence_length)
  one_hot_encoding = OneHotEncoder(
      sparse=False, n_values=5, dtype=integer_type).fit_transform(integer_array)

  return one_hot_encoding.reshape(len(sequences), sequence_length, letters_length, 1).swapaxes(1, 2)
# This could just be ambiguous_dna_letters, but that would be much higher dim.
class IUPACUnambiguousDNAWithN(Alphabet):
  letters = unambiguous_dna_letters + "N"

def encode_fasta_sequence(fname, letters='ATCGN'):

def encode_fasta_sequence(fname):
  """
  Loads fasta file and returns an array of one-hot sequences.

@@ -171,36 +182,12 @@ def encode_fasta_sequence(fname, letters='ATCGN'):
  fname: str
    Filename of fasta file.
  """
  name, seq_chars = None, []
  sequences = []
  with open(fname) as fp:
    for line in fp:
      line = line.rstrip()
      if line.startswith(">"):
        if name:
          sequences.append(''.join(seq_chars).upper())
        name, seq_chars = line, []
      else:
        seq_chars.append(line)
  if name is not None:
    sequences.append(''.join(seq_chars).upper())

  return seq_one_hot_encode(np.array(sequences), letters)

def encode_sequence(fname, letters='ATCGN', file_type="fasta"):
  dna_alphabet = IUPACUnambiguousDNAWithN()
  return encode_bio_sequence(fname, "fasta", dna_alphabet)

    # TODO: if None, then get from the filename
    sequences = []
    for seq in SeqIO.parse(fname, file_type):
        sequences.append(str(seq.seq).upper())

    return seq_one_hot_encode(np.array(sequences), letters)

# This could just be ambiguous_dna_letters, but that would be much higher dim.
class IUPACUnambiguousDNAWithN(Alphabet):
    letters = unambiguous_dna_letters + "N"

def encode_sequence_with_biopython(fname, file_type="fasta", alphabet=None):
def encode_bio_sequence(fname, file_type="fasta", alphabet=None):
  # np.ndarray: Shape (N_sequences, 4, sequence_length, 1).

  if alphabet is None:
@@ -208,31 +195,7 @@ def encode_sequence_with_biopython(fname, file_type="fasta", alphabet=None):

  # TODO: if None, then get from the filename
  sequences = SeqIO.parse(fname, file_type, alphabet)

    # The label encoder is given characters for ACGTN
    letter_encoder = {l: i for i, l in enumerate(alphabet.letters)}
    alphabet_length = len(letter_encoder)

    # Peak at the first sequence to get the length of the sequence.
    first_seq = next(sequences)
    sequence_length = len(first_seq.seq)

    seqs = []

    seqs.append(_seq_to_encoded(first_seq, letter_encoder, alphabet_length, sequence_length))

    for other_seq in sequences:
      seqs.append(_seq_to_encoded(other_seq, letter_encoder, alphabet_length, sequence_length))

    # return np.expand_dims(np.array(seqs), -1)
    return np.array(seqs)

def _seq_to_encoded(seq, letter_encoder, alphabet_length, sequence_length):
    b = np.zeros((alphabet_length, sequence_length))
    seq_ints = [letter_encoder[s] for s in seq]
    b[seq_ints, np.arange(sequence_length)] = 1

    return b
  return seq_one_hot_encode(sequences, alphabet.letters)


def save_metadata(tasks, metadata_df, data_dir):
+13 −15
Original line number Diff line number Diff line
@@ -20,26 +20,24 @@ class TestSeq(unittest.TestCase):
  def test_one_hot_simple(self):
    sequences = np.array(["ACGT", "GATA", "CGCG"])
    sequences = dc.utils.save.seq_one_hot_encode(sequences)
    assert sequences.shape == (3, 4, 4, 1)
    self.assertEqual(sequences.shape, (3, 5, 4, 1))

  def test_one_hot_mismatch(self):
    # One sequence has length longer than others. This should throw a
    # value error.
    thrown = False
    try:
    # ValueError.

    with self.assertRaises(ValueError):
      sequences = np.array(["ACGTA", "GATA", "CGCG"])
      sequences = dc.utils.save.seq_one_hot_encode(sequences)
    except ValueError:
      thrown = True
    assert thrown

  def test_encode_sequence_with_biopython(self):
  def test_encode_fasta_sequence(self):
    fname = "./data/example.fasta"

      encoded_seqs = dc.utils.save.encode_sequence_with_biopython(fname)
      expected = np.array([
    encoded_seqs = dc.utils.save.encode_fasta_sequence(fname)
    expected = np.expand_dims(
        np.array([
            [[0, 0], [1, 0], [0, 0], [0, 1], [0, 0]],
            [[1, 0], [0, 1], [0, 0], [0, 0], [0, 0]],
      ])
        ]), -1)

    np.testing.assert_array_equal(expected, encoded_seqs)