Внедрение ДНК кодон таблицы в C


Я решил переделать первую программу я написал (после Здравствуйте, дело мира). Так, что она делает это, чтобы преобразовать последовательность ДНК строку на массив строк, которые представляют собой аминокислоты, которые вносятся в клетке на рибосомах при том, что последовательности ДНК привыкает (ДНК кодон таблица).

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

char* codon_table[4][4][4] =
{
  {
    {"Phe", "Phe", "Len", "Len"},
    {"Ser", "Ser", "Ser", "Ser"},
    {"Tyr", "Tyr", "Stop", "Stop"},
    {"Cys", "Cys", "Stop", "Trp"}
  },
  {
    {"Len", "Len", "Len", "Len"},
    {"Pro", "Pro", "Pro", "Pro"},
    {"His", "His", "Gln", "Gln"},
    {"Arg", "Arg", "Arg", "Arg"}
  },
  {
    {"Ile", "Ile", "Ile", "Met"},
    {"Thr", "Thr", "Thr", "Thr"},
    {"Asn", "Asn", "Lys", "Lys"},
    {"Ser", "Ser", "Arg", "Arg"}
  },
  {
    {"Val", "Val", "Val", "Val"},
    {"Ala", "Ala", "Ala", "Ala"},
    {"Asp", "Asp", "Glu", "Glu"},
    {"Gly", "Gly", "Gly", "Gly"}
  }
};

int check_nucleobase_string(char* nucleobase_string)
{
  if(nucleobase_string==NULL)
  {
    return 1;
  }
  for(size_t i = 0; nucleobase_string[i] != '\0'; i+=1)
  {
    if(
      nucleobase_string[i] != 't' &&
      nucleobase_string[i] != 'c' &&
      nucleobase_string[i] != 'a' &&
      nucleobase_string[i] != 'g'
    )
    {
      return 0;
    }
  }
  return 1;
}

size_t nucleobase_to_aminoacid(char* nucleobase_string, char*** aminoacid_string)
{
  if(!check_nucleobase_string(nucleobase_string))
  {
    printf("Erroneous input: %s\n", nucleobase_string);
    return 0;
  }

  size_t nucleobase_count ;
  for(nucleobase_count = 0; nucleobase_string[nucleobase_count] != '\0'; nucleobase_count+=1);
  size_t aminoacid_count = nucleobase_count/3;

  (*aminoacid_string) = malloc(aminoacid_count);

  for(size_t i = 0, slow_i = 0; slow_i < aminoacid_count; i+=3, slow_i+=1)
  {
    (*aminoacid_string)[slow_i] =
      codon_table[
        nucleobase_string[i] == 't' ? 0 :
        nucleobase_string[i] == 'c' ? 1 :
        nucleobase_string[i] == 'a' ? 2 :
        nucleobase_string[i] == 'g' ? 3 : -1
      ][
        nucleobase_string[i+1] == 't' ? 0 :
        nucleobase_string[i+1] == 'c' ? 1 :
        nucleobase_string[i+1] == 'a' ? 2 :
        nucleobase_string[i+1] == 'g' ? 3 : -1
      ][
        nucleobase_string[i+2] == 't' ? 0 :
        nucleobase_string[i+2] == 'c' ? 1 :
        nucleobase_string[i+2] == 'a' ? 2 :
        nucleobase_string[i+2] == 'g' ? 3 : -1
      ];
  }

  return aminoacid_count;
}

int main()
{
  printf("Enter DNA sequence in 5' → 3' direction.\n");

  char* input_string = NULL;
  size_t n = 0;
  ssize_t line_length = getline(&input_string, &n, stdin);
  if(line_length == -1)
  {
    printf("Failed to read a line.\n");
    return 0;
  }

  char* nucleobase_string = malloc(line_length);
  memcpy(nucleobase_string, input_string, line_length-1);
  nucleobase_string[line_length-1] = '\0';

  char** aminoacid_string = NULL;
  size_t aminoacid_count = nucleobase_to_aminoacid(nucleobase_string, &aminoacid_string);

  for(size_t i = 0; i < aminoacid_count; i+=1)
  {
    if(i!=0)
    {
      printf(" + ");
    }
    printf("%s", aminoacid_string[i]);
  }
  printf("\n");

  return 0;
};


1514
9
задан 11 апреля 2018 в 01:04 Источник Поделиться
Комментарии
6 ответов

Массив char* Не многие структуры данных для хранения AA или кодон последовательности данных. Указатель занимает 4 или 8 байт на каждый кодон, но есть только 64 возможных кодонов, чтобы вы могли легко хранить один кодон в байте uint8_t. И лишь около 22 аминокислот, который близок к тому, чтобы соответствовать 2 ААС на байт, но он не совсем работает, и это наверняка проще все-таки использовать uint8_t.

Затем, когда вы хотите напечатать, используйте это число в качестве индекса в таблице строк.

Или если все вы хотите сделать, это распечатать строку, тогда не построить массив строк вовсе, а просто printf, как вы преобразовать. Или memcpy 4-байтовые строки, как "Phe " в буфер, чтобы создать строку самостоятельно, вместо вызова printf отдельно для каждого АА имя.

(Компиляторы хороши, в 4-байтовое копии, потому что это размер целочисленный регистр. Особенно от постоянных данных, потому что они могут встраивать информацию в качестве непосредственного, если это оптимально. Так что делайте таблицу кодонов static const)


slow_i странное имя переменной. Использовать aa_idx или еще что-то значимое, а может j. Или использовать один i переменной и использовать i*3 + 0, i*3 + 1и i*3 + 2 чтобы дать компилятору понять, что он хочет сделать.


Это очень длинная строка: for(nucleobase_count = 0; nucleobase_string[nucleobase_count] != '\0'; nucleobase_count+=1); должно быть nucleobase_count = strlen(nucleobase_string);. Использовать библиотечные функции вместо того, чтобы повторно изобретать колесо трудно читать длинные линии. Строковые функции библиотеки libc, как правило, вручную оптимизированные реализации ПКР, которые будут быстро, особенно для длинных строк.

Или еще лучше, иметь свой нуклеотидных проверки функция возвращает число , если не найти проблему! Ты уже обходишь всю строку; вы не должны делать это снова, чтобы подсчитать длину.


Или еще лучше, у вашего собеседника передать вам длины , так что вы можете выделить память сразу и совместить проверки с преобразованием. Вы могли бы абонент мимо вас уже выделил буфер вывода, а также длину, поэтому он имеет возможность использовать С99 переменной длины массива, повторно использовать буфер, он уже выделил или что-то, вместо выпечки malloc на ваше обращение рутины.


 codon_table[  // your code
nucleobase_string[i] == 't' ? 0 :
nucleobase_string[i] == 'c' ? 1 :
nucleobase_string[i] == 'a' ? 2 :
nucleobase_string[i] == 'g' ? 3 : -1
][...][...]

Если плохой нуклеотидных сделал скольжения через трещины и сделать его на этот код, вы бы неопределенное поведение от чтения codon_table[-1][0][0]вне границ допуска. Вы могли бы снять == g часть и просто 3 последний шаг в тернарный оператор, с комментарием, чтобы указать, что это g если бы не один из других 3. Кроме того, что блок тернарные операторы должны идти в крошечной вспомогательной функции, как int base_to_idx(char base); вместо того, чтобы быть повторен 3 раза.

Затем вы можете повторно реализовать его с помощью таблиц или идеальной хэш-функции. например

// probably not faster than a table lookup, unless you make a SIMD version
/*
'a' 0b1100001 -> 0b00
'c' 0b1100011 -> 0b10
'g' 0b1100111 -> 0b01
't' 0b1110100 -> 0b11
*/
unsigned nuc_to_idx(unsigned base) {
if (base == 'a' || base == 'c' || base == 'g' || base == 't')
return (base*3 & 0x3<<2)>>2;
return -128U;
}

или таблицу подстановки версия, которая делает все кодона и возвращает целочисленный код. Я с места-инициализатора синтаксис, плюс расширение GNU для диапазонов , чтобы сделать таблицу static const. Вы могли бы также инициализации во время выполнения с memset(table, -1, UCHAR_MAX); а затем назначить 4 допустимые значения.

#include <stdint.h>
#include <limits.h>

// int8_t sign extends to the width of int / unsigned int
// so all the high bits get set before treating as unsigned (assuming 2's complement), for BAD_NUC inputs.
// But we don't depend on that: we still get an unsigned value > 64
// on sign/magnitude or one's complement machines if any of the 3 chars is bogus.

#define BAD_NUC -128
static const int8_t nuc_table[UCHAR_MAX+1] = {
[0 ... 255] = BAD_NUC, // ranges are a GNU extension
// last init takes precedence https://gcc.gnu.org/onlinedocs/gcc/Designated-Inits.html
['a'] = 0,
['c'] = 1,
['g'] = 2,
['t'] = 3,
};
// Fun fact: doesn't depend on ASCII encoding because of the designated-initializers.

// unsigned char* so high-ASCII -> 128..255, not negative,
// and works as an index into a 256 entry table
unsigned codon_to_idx_LUT(unsigned char *p) {

unsigned idx = nuc_table[p[0]];
idx = idx*4 + nuc_table[p[1]];
idx = idx*4 + nuc_table[p[2]];
return idx;
// 0..63 for in-range values, otherwise has higher bits set
// so you only have to range-check once per codon, speeding up the common case of no error.

// to map to AA codes, maybe use a lookup table on this 0..63 value
}

В таблице версия компилирует очень эффективно на x86 с НКУ -О3 (на Godbolt компилятор проводника), или, Надеюсь, с любой вменяемый компилятор. (Условие if в идеальной хэш-версия компиляции к немедленному точечный рисунок, который аккуратно, но это еще предстоит проверить/филиал несколько раз для каждого персонажа, а не просто делать 2 нагрузок.)

codon_to_idx_LUT:                # gcc7.3 -O3 for x86-64 System V calling convention: pointer arg in RDI
movzx eax, BYTE PTR [rdi]
movsx edx, BYTE PTR nuc_table[rax]
movzx eax, BYTE PTR [rdi+1]
movsx eax, BYTE PTR nuc_table[rax]
lea edx, [rax+rdx*4] # this is why I wrote it as idx*4 + nuc_table[...] instead of some other way of multiplying and adding the results.
movzx eax, BYTE PTR [rdi+2]
movsx eax, BYTE PTR nuc_table[rax]
lea eax, [rax+rdx*4]
ret

Это будет работать даже быстрее, если компилятор сделал один широкий нагрузки и распаковали его с инструкциями АЛУ вместо загрузки каждой строки байт отдельно. Работает это в цикле будет узким местом на 2 нагрузок в часы на основные процессоров Intel/AMD, так что о 1 нуклеотид в такт, что тоже не плохо. Накладных расходов на проверку индекса в диапазон (для выявления фиктивных нуклеотидных букв) должна быть незначительной по сравнению с строку нагрузок + поиск по таблице. Особенно на процессорах, которые можете сделать только одно нагрузка на часы.

(Связанные с этим вопросы и ответы Все о настройки вашего источника для рук компилятор сделать лучше АСМ. Вы могли бы сделать это здесь, но только ценой читабельности, я думаю. Вы, вероятно, нужно сделать хитроумный указатель-литой или безопасный memcpy и читать 4 байта как uint32_t и распаковывать.)

Вы можете использовать это 0..63 кодон как индекс в таблице строк для печати. (Просто расплющит текущей 3D-массив, хотя я не делал никаких усилий, чтобы быть в соответствии с которой номера идут с какой базы ДНК в хэш-таблице или против вашей версии.)

11
ответ дан 11 апреля 2018 в 10:04 Источник Поделиться


  • Линии

        (*aminoacid_string) = malloc(aminoacid_count);

    выделяет aminoacid_count из байтов. Код должен, что многие указатели. Правильный способ-это

        *aminoacid_string = malloc(aminoacid_count * sizeof(**aminoacid_string));

  • Проверка на правильный ввод (т. е. вызов check_nucleobase_string) относится к основным.

  • Я рекомендую перевести азотистые строку в числовой массив (сопоставление t для 0 и т. д.), и избежать некрасивой каскад тернарных операторов.

  • Цикл

        for(nucleobase_count = 0; nucleobase_string[nucleobase_count] != '\0'; nucleobase_count+=1);

    очень долго (и неэффективный) способ сказать

        nucleobase_count = strlen(nucleobase_string);

  • Я не вижу причин

        memcpy(nucleobase_string, input_string, line_length-1);

    Вы можете работать напрямую с input_string.


7
ответ дан 11 апреля 2018 в 01:04 Источник Поделиться

Вы должны попытаться выполнить рефакторинг кода, чтобы менее жестко заданные константы. Е. Г. nucleobase_to_aminoacid имеет как tcag и codon_table жестко.
Это в общем-то, что препятствует повторному использованию.


  • Вы должны предпочтительно иметь отдельную функцию для печати codon_table - в стандартном формате. Это бы обнаружил, что вы использовали "лень" вместо "Ноу" как сокращение от "лейцин".

  • Если расшифровывать tcag был отдельной функции или ввода - вы могли бы повторно использовать ту же функцию для мРНК.

  • Если codon_table был вход в nucleobase_to_aminoacid вы можете повторно использовать тот же код для нечетных случаях с разными codon_table.

7
ответ дан 11 апреля 2018 в 09:04 Источник Поделиться

В дополнение к текущих (и будущих) ответы:


  1. Вы используете malloc()но от того, что я вижу, вы не free() он позже. В моем случае с Valgrind жалуется, что 125 байт в 2 кварталах "однозначно потеряли", - значит у вас есть утечки памяти.

  2. Это распространенный способ использования ++i вместо i += 1особенно в для петли.

  3. Используйте const классификатор, где это применимо. В то время как компилятор может делать такие оптимизации для вас (хотя не так часто, как можно подумать), это полезно, чтобы указать, что что-то не изменится.

    Например, char *codon_table может быть char const *const codon_table. Если вы используете -Wconversion флаг при компиляции с помощью GCC, вы можете поймать всевозможные ошибки связано с типом преобразований.


6
ответ дан 11 апреля 2018 в 01:04 Источник Поделиться

Не много осталось пересмотреть с учетом других хороших ответов. Некоторые дополнительные идеи:

Срубая в '\n'

Конечно, ОП хотел отрубить обычный трейлинг '\n'. Ниже отрежут последний символ из строки ввода, независимо от того, если это '\n' или нет.

nucleobase_string[line_length-1] = '\0';  // May lop off something else.

Рассмотрим конца трубного ввода пользователя могут не закончиться '\n'так, хотя трейлинг '\n' очень часто, это не cerrtainty. Рекомендую проверить '\n'.

if (nucleobase_string[line_length-1] == '\n') {
nucleobase_string[line_length-1] = '\0';
}
// or the one liner (slower)
nucleobase_string[strcspn(nucleobase_string, "\n")] = '\0';

Дома держать

Хорошая кодирование практика бесплатно input_string в конце концов. Конечно, такое распределение неявно освобождается, когда программа заканчивается. Еще с явным free()кто-то повторно, используя этот фрагмент кода будет выгодно.

free(input_string);
return 0;

3
ответ дан 12 апреля 2018 в 02:04 Источник Поделиться

Использование stderr

В нескольких местах вашей программы вы используете printf(...) к:


  • спросить пользователя, чтобы обеспечить ввод

  • вывод сообщения об ошибке

  • конечный результат

printf функция пишет в stdout. Вы могли бы рассмотреть, используя fprintf(stderr, ...) спрашиваю для ввода данных и журнала ошибок.

Статус выхода

Также он является общим, чтобы возвратить ненулевое информация о программной ошибке вызова exit(EXIT_FAILURE); или return EXIT_FAILURE; от main().

Таким образом можно использовать вашу программу, как:

#!/bin/sh
result="$(echo 'whatever' | ./your_program 2>/dev/null)"
if [ $? -ne 0 ]; then
# failed
fi

3
ответ дан 12 апреля 2018 в 03:04 Источник Поделиться