méthode pour encoder une séquence ADN

Version imprimable

Voir 40 message(s) de cette discussion en une page

15/03/2016, 18h47
bacelar

Citation:

Le "__" est utilisé pour les conventions d'appel sous Zindows (MSVC) ex. __cdecl, __fastcall ...

Ils se permettent cela justement parce que la norme leur en donne le droit en C++.
VC++ est un compilateur C++, pas un compilateur C.
Mais je pense qu'il y a les même restrictions en C (tout OS, tout compilateur, même s'ils ne le vérifient pas).
15/03/2016, 18h51
Bousk

Quitte à disgresser jusqu'au bout:
http://www.doc.ic.ac.uk/lab/cplus/c++.rules/chap5.html

Citation:

The use of two underscores (`__') in identifiers is reserved for the compiler's internal use according to the ANSI-C standard.

Underscores (`_') are often used in names of library functions (such as "_main" and "_exit"). In order to avoid collisions, do not begin an identifier with an underscore.

http://www.gnu.org/software/libc/man...ved-Names.html

Citation:

In addition to the names documented in this manual, reserved names include all external identifiers (global functions and variables) that begin with an underscore (‘_’) and all identifiers regardless of use that begin with either two underscores or an underscore followed by a capital letter are reserved names. This is so that the library and header files can define functions, variables, and macros for internal purposes without risk of conflict with names in user programs.

En gros c'est réservé à l'implémentation et aux mot-clés du compilo.
16/03/2016, 07h14
FryHandiz

Woaw*! D’habitude quand je poste une question sur un forum j’ai soit jamais de réponses, soit on me répond que mon message n’est pas «*conforme*»*:D
Bon ok, j’avais aussi perdu le chemin pour venir ici …*! (non je ne suis pas blonde du tout*! Ahem.)

Bon alors, j’ai demandé un peu plus d’éclaircissement, et donc il faut pouvoir «*encoder*» de sorte à ce qu’on utilise 8 bits pour 4 lettres dans chacune des cases de notre tableau où vont se retrouver nos A,T,C,G.

@Koala01 > Dans mon cas, on souhaite surtout travailler sur des séquences même si, plus tard, on risque de vouloir rechercher si, entre deux séquences, on ne retrouve pas de similarité, etc. Mais pour moi, ça ce sera peut être plus tard ou jamais vu ma très faible connaissance en C++ xD (C’est la première année où je fais de l’informatique et on nous demande de coder une API de mapping de séquence… AHEM. Vive la reconversion…*! Heureusement que je suis un peu maso*:p).

@RPGamer > Si je ne dis pas de bétises, pour le complément/manipuler les séquences, il est préférable de le faire une fois qu’on a encodé. Car certaines fois, les séquences de nucléotides peuvent aller jusqu'à… 120 000 nucléotides, ce qui apparemment peut saturer la mémoire si on ne passe pas avant à l’encodage.

Me frappez pas mais pour la suite je n’ai pas trop compris le coup du A = 0x01, etc.*:s
Pourquoi devoir faire 0x01…*? Et pourquoi pas 1*? Et aussi, pourquoi de 2 en 2…*?
Je suis vraiment désolée si mes questions sont stupides mais je suis malheureusement loin d’être une informaticienne >.<’

Encore merci pour votre aide*!
16/03/2016, 09h07
RPGamer

Quand tu parles d'encoder les nucleotides sur un octet c'est dans cette direction que tu dois aller. On utilise la notation hexadécimal (0x..) pour coder un octet pour des questions pratiques (2⁴ = 16 donc 2 digits hexa permettent de coder 1 octet).

Donc pour coder une amine on te propose d'utiliser le premier bit :

0x01 = 00000001

pour coder une cétone le 2e bit :

0x02 = 00000010

pour coder une purine le 3e bit :

0x04 = 00000100

et pour coder une pyrimidine le 4e bit :

0x08 = 00001000

tous les bits du premier digit à 1 pour le masque de complément :

0x0f = 00001111

Avec ces valeurs tu vas pouvoir coder A, C, G et T par une simple opération logique "OU" (|), c-à-d pour chaque position si au moins 1 bit sur l'un des octets vaut 1, le bit du résultat à la même position vaut 1. Table de vérité :

1 OU 1 = 1
1 OU 0 = 1
0 OU 1 = 1
0 OU 0 = 0

Exemple :

A = amine + purine = 0x01 | 0x04 =

00000001
00000100 |
________
00000101 = 0x05

on a donc A = 0x05.

Pour le complément il suffit d'appliquer un "OU exclusif" (^) avec le masque de complément, c-à-d une opération OU dont un bit du résultat à une position donnée ne peut valoir 1 que si un seul des deux opérande a un bit 1 à cette même position. Table de vérité :

1 OUX 1 = 0
1 OUX 0 = 1
0 OUX 1 = 1
0 OUX 0 = 0

Le complément de A peut se déterminer comme suit : A OUX complément = 0x05 ^ 0x0f =

00000101
00001111 ^
________
00001010 = 0x0a

Le complément de A est donc T (cétone + pyrimidine).

Etant donné qu'il suffit de 4 bits pour coder tous les nucleotides, si tu veux absolument économiser de la mémoire il est toujours possible de coder 2 nucleotides par octets comme ceci :

A = 0x05 = 00000101
C = 0x09 = 00001001

AC = (A << 4) | C = 0x59 = 01011001

et pour extraire C par exemple, il faut appliquer un opération "ET" (&), c'est à dire que le bit du résultat vaut 1 seulement si les bits des 2 opérandes de cette position valent 1. Table de vérité :

1 ET 1 = 1
1 ET 0 = 0
0 ET 1 = 0
0 ET 0 = 0

C = AC ET complément = 0x59 & 0x0f =

01011001
00001111 &
________
00001001 = 0x09

Pour extraire A, étant donné qu'il est positionné au niveau des bits de poids fort, un simple décalage de 4 bits vers la droite suffit, les ordinateurs sont très performants pour le faire :

A = AC >> 4 = 0x59 >> 4 = 00000101 = 0x05

Après si ton max de mémoire occupée c'est 120 Ko, pas de quoi faire tousser nos machines ;)
16/03/2016, 09h15
ternel

D'une manière générale, on préfère utiliser la notation hexadécimale pour les valeurs subissant des opérations bit-à-bit, c'est à dire &, |, ^ et ~ (inverse).
Cela souligne que ce sont les bits du nombre qui nous intéresse, et qu'on va faire des choses non algébriques dessus.
16/03/2016, 11h50
fenkys

Bonjour,

120 000 nucléotide ce n'est rien avec un ordinateur moderne.

Pour répondre à RPGamer, pourquoi utiliser un bit différent pour une amine et une cétone. Le nucléotide est soit l'un, soit l'autre. Idem pour les purines / pyrimidines. Tu peux donc réduire à 2 bits par bases et 4 bases par octets.

Ceci étant, rechercher un séquence avec ce genre de stockage peut être un poil compliqué. Mais si la recherche n'est pas sa priorité, c'est jouable.

Citation:

on utilise 8 bits pour 4 lettres dans chacune des cases de notre tableau

Et si j'ai bien traduit cette phrase, c'est ce qui est demandé.

OK bah peu importe, la méthode reste la même, il suffit juste de diviser le nombre d'énumérés par 2. Si je reste sur une implémentation avec 4 codes :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
#include <iostream>
#include <vector>
#include <map>
#include <algorithm>
 
enum Nucleotide : char
{
    AMINE      = 0x01,
    KETONE     = 0x02,
    PURINE     = 0x04,
    PYRIMIDINE = 0x08,
 
    COMPLEMENT = 0x0f,
 
    A = AMINE | PURINE,
    C = AMINE | PYRIMIDINE,
    G = KETONE | PURINE,
    T = KETONE | PYRIMIDINE
};
 
inline const Nucleotide operator ^=(Nucleotide &x, const Nucleotide y)
{
    return x = static_cast<Nucleotide>(x ^ y);
}
 
inline const Nucleotide operator ^(Nucleotide x, const Nucleotide y)
{
    return x ^= y;
}
 
inline std::ostream &operator <<(std::ostream &os, const Nucleotide n)
{
    static const std::map<Nucleotide, char> MAPPING
    {{
        {A, 'A'},
        {C, 'C'},
        {G, 'G'},
        {T, 'T'}
    }};
 
    return os << MAPPING.at(n);
}
 
class Sequence
{
private:
 
    std::vector<Nucleotide> sequence;
 
public:
 
    Sequence(const std::vector<Nucleotide> &sequence) :
        sequence(sequence)
    {
 
    }
 
    void reverse()
    {
        std::reverse(sequence.begin(), sequence.end());
    }
 
    void complement()
    {
        for (auto &n : sequence)
        {
            n ^= COMPLEMENT;
        }
    }
 
    friend std::ostream &operator <<(std::ostream &os, const Sequence &s)
    {
       for (const auto n : s.sequence)
       {
          os << n;
       }
       return os;
    }
};
 
int main()
{
   Sequence sequence{{A, G, C, T, G, T, A, C}};
   std::cout << "Sequence de base : " << sequence << std::endl;
 
   Sequence sequence2 = sequence;
 
   sequence.reverse();
   std::cout << "Inverse : " << sequence << std::endl;
 
   sequence2.complement();
   std::cout << "Complement : " << sequence2 << std::endl;
 
   Nucleotide complement = A ^ COMPLEMENT;
   std::cout << "Complement de " << A << " : " << complement << std::endl;
 
   return 0;
}

Il faut ensuite l'adapter pour n'utiliser que 2 codes.

On peut aussi mixer macros temporaires et énumérations. Par exemple:

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
//Header file: nucleotide.h
 
#ifndef NUCLEOTIDE_H
#define NUCLEOTIDE_H
 
#include <stdexcept>
 
#ifdef NUCLEOTIDE_F_AMINE
# error "NUCLEOTIDE_F_AMINE is already defined!"
#endif // NUCLEOTIDE_F_AMINE
#ifdef NUCLEOTIDE_F_KETONE
# error "NUCLEOTIDE_F_KETONE is already defined!"
#endif // NUCLEOTIDE_F_KETONE
#ifdef NUCLEOTIDE_F_PURINE
# error "NUCLEOTIDE_F_PURINE is already defined!"
#endif // NUCLEOTIDE_F_PURINE
#ifdef NUCLEOTIDE_F_PYRIMIDINE
# error "NUCLEOTIDE_F_PYRIMIDINE is already defined!"
#endif // NUCLEOTIDE_F_PYRIMIDINE
#ifdef NUCLEOTIDE_M_COMPLEMENT
# error "NUCLEOTIDE_M_COMPLEMENT is already defined!"
#endif // NUCLEOTIDE_M_COMPLEMENT
 
//New temporary macros
#define NUCLEOTIDE_F_AMINE      0x01
#define NUCLEOTIDE_F_KETONE     0x02
#define NUCLEOTIDE_F_PURINE     0x04
#define NUCLEOTIDE_F_PYRIMIDINE 0x08
#define NUCLEOTIDE_M_COMPLEMENT 0x0f
 
//New unscoped enumeration: Nucleotide
enum Nucleotide : char
{
  Adenine  = NUCLEOTIDE_F_AMINE  | NUCLEOTIDE_F_PURINE,
  Cytosine = NUCLEOTIDE_F_AMINE  | NUCLEOTIDE_F_PYRIMIDINE,
  Guanine  = NUCLEOTIDE_F_KETONE | NUCLEOTIDE_F_PURINE,
  Thymine  = NUCLEOTIDE_F_KETONE | NUCLEOTIDE_F_PYRIMIDINE
}
 
//New function: isAmine
inline bool isAmine(Nucleotide const & n) noexcept
{
  return n & NUCLEOTIDE_F_AMINE;
}
 
//New function: isKetone
inline bool isKetone(Nucleotide const & n) noexcept
{
  return n & NUCLEOTIDE_F_KETONE;
}
 
//New function: isPurine
inline bool isPurine(Nucleotide const & n) noexcept
{
  return n & NUCLEOTIDE_F_PURINE;
}
 
//New function: isPyrimidine
inline bool isPyrimidine(Nucleotide const & n) noexcept
{
  return n & NUCLEOTIDE_F_PYRIMIDINE;
}
 
//New function: areComplement
inline bool areComplement(Nucleotide const & n1, Nucleotide const & n2) noexcept
{
  return (n1 ^ n2) == NUCLEOTIDE_M_COMPLEMENT;
}
 
//New function: complementOf
inline Nucleotide complementOf(Nucleotide const & n) noexcept
{
  return static_cast< Nucleotide >(n ^ NUCLEOTIDE_M_COMPLEMENT);
}
 
//New function: charToNucleotide
Nucleotide charToNucleotide(char const & c)
{
  Nucleotide n = Adenine;
  switch(c)
  {
    case 'A':
    case 'a':
      break;
    case 'C':
    case 'c':
      n = Cytosine;
      break;
    case 'G':
    case 'g':
      n = Guanine;
      break;
    case 'T':
    case 't':
      n = Thymine;
      break;
    default:
      throw std::invalid_argument("Not a DNA nucleotide!");
      break;
  }
  return n;
}
 
//New function: nucleotideToChar
inline char nucleotideToChar(Nucleotide const & n) noexcept
{
  return (n & NUCLEOTIDE_F_AMINE) ?
             ((n & NUCLEOTIDE_F_PURINE) ? 'A' : 'C') :
             ((n & NUCLEOTIDE_F_PURINE) ? 'G' : 'T');
}
 
#undef NUCLEOTIDE_F_AMINE
#undef NUCLEOTIDE_F_KETONE
#undef NUCLEOTIDE_F_PURINE
#undef NUCLEOTIDE_F_PYRIMIDINE
#undef NUCLEOTIDE_M_COMPLEMENT
 
#endif // NUCLEOTIDE_H

Je me suis permis de remplacer A, C, G et T par Adenine, Cytosine, Guanine et Thymine du fait que c'est moins risqué en ce qui concerne les collisions de nommage.

17/03/2016, 14h16
bacelar

Citation:

On peut aussi mixer macros temporaires et énumérations.

Moi, je ne vois que des inconvénients à l'utilisation des MACRO ici, perte de typage, namespace aux fraises, etc...
17/03/2016, 14h43
VivienD

Citation:

Envoyé par bacelar

Moi, je ne vois que des inconvénients à l'utilisation des MACRO ici, perte de typage, namespace aux fraises, etc...

C'est vrai que les macros sont insensibles aux espaces de nommage et que le typage avec les macros peut relever de la gageure. Néanmoins, l'énumération Nucleotide n'a, grâce aux macros, que quatre valeurs disponibles, à savoir Adenine, Cytosine, Guanine et Thymine. Ainsi, les chances que l'utilisateur entre une valeur erronée, comme un drapeau seul ou un masque, au lieu d'un nucléotide valide sont réduites à zéro, ou presque (cf. loi de Murphy).

Dans ce cas-là, on préférera utiliser des expressions constantes.

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
#include <iostream>
 
namespace DNA
{
    constexpr char NUCLEOTIDE_F_AMINE      = 0x01;
    constexpr char NUCLEOTIDE_F_KETONE     = 0x02;
    constexpr char NUCLEOTIDE_F_PURINE     = 0x04;
    constexpr char NUCLEOTIDE_F_PYRIMIDINE = 0x08;
    constexpr char NUCLEOTIDE_M_COMPLEMENT = 0x0f;
 
    enum Nucleotide : char
    {
        Adenine  = NUCLEOTIDE_F_AMINE  | NUCLEOTIDE_F_PURINE,
        Cytosine = NUCLEOTIDE_F_AMINE  | NUCLEOTIDE_F_PYRIMIDINE,
        Guanine  = NUCLEOTIDE_F_KETONE | NUCLEOTIDE_F_PURINE,
        Thymine  = NUCLEOTIDE_F_KETONE | NUCLEOTIDE_F_PYRIMIDINE
    };
}
 
int main()
{
    auto n = DNA::Adenine;
    std::cout << "0x0" << std::hex << n << std::endl;
 
    return 0;
}

Citation:
Envoyé par RPGamer
Dans ce cas-là, on préférera utiliser des expressions constantes.
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 #include <iostream> namespace DNA { constexpr char NUCLEOTIDE_F_AMINE = 0x01; constexpr char NUCLEOTIDE_F_KETONE = 0x02; constexpr char NUCLEOTIDE_F_PURINE = 0x04; constexpr char NUCLEOTIDE_F_PYRIMIDINE = 0x08; constexpr char NUCLEOTIDE_M_COMPLEMENT = 0x0f; enum Nucleotide : char { Adenine = NUCLEOTIDE_F_AMINE | NUCLEOTIDE_F_PURINE, Cytosine = NUCLEOTIDE_F_AMINE | NUCLEOTIDE_F_PYRIMIDINE, Guanine = NUCLEOTIDE_F_KETONE | NUCLEOTIDE_F_PURINE, Thymine = NUCLEOTIDE_F_KETONE | NUCLEOTIDE_F_PYRIMIDINE }; } int main() { auto n = DNA::Adenine; std::cout << "0x0" << std::hex << n << std::endl; return 0; }
Certes, ca règle aussi les problèmes de validité mais l'utilisateur final peut encore accéder aux drapeaux seuls et au masque. Si ca ne dérange pas, on peut laisser comme ca (ou même envoyer tout ca dans un sous-espace de nommage du genre DNA::bits puis regrouper les drapeaux et les masques dans les énumérations Nucleotide_Flags et respectivement Nucleotide_Masks); dans le cas contraire, je ne vois que les macros temporaires pour atteindre l'effet voulu.

17/03/2016, 16h01
RPGamer

~~Les macros n'empêchent pas plus l'accès aux valeurs par l'utilisateur, au contraire même.~~ J'avais pas vu que tu les désactive à la fin de ton fichier source.

Il y a une autre solution qui consiste à transformer DNA en classe principale et à placer les expressions constantes comme membres statiques privés de cette classe et les fonctions libres comme fonctions amies.
17/03/2016, 16h12
VivienD

Citation:

Envoyé par RPGamer

Les macros n'empêchent pas plus l'accès aux valeurs par l'utilisateur, au contraire même.

Regarde la fin de mon code aux lignes 112 à 116 et tu verras que je les "détruis" à coup de #undef. Et oui! ca marche: je me suis déjà servi de cette astuce à de maintes et maintes reprises.

EDIT:

Citation:

Envoyé par RPGamer

~~Les macros n'empêchent pas plus l'accès aux valeurs par l'utilisateur, au contraire même.~~ J'avais pas vu que tu les désactive à la fin de ton fichier source.

Il y a une autre solution qui consiste à transformer DNA en classe principale et à placer les expressions constantes comme membres statiques privés de cette classe et les fonctions libres comme fonctions amies.

C'est vrai; je n'y avais pas pensé.:aie:
20/03/2016, 17h51
bowow

Bonjour,

Je ne sais pas trop ce que tu recherches à faire, mais si la finalité est de faire de l'alignement de séquence, on ne peut se permettre de les stocker en mémoire vu la quantité de séquences dans les fichiers FASTQ et FASTA.

Par contre, si c'est un simple exercice pour enregistrer une séquence ADN en mémoire, trouver la séquence inverse, complémentaire, la séquence ARNm et polypeptidique, t'embêtes pas avec ces encodages.

De plus, il existe déjà des algorithmes de texte en bio-informatique et une multitude d'outils déjà préexistants qui feront sûrement ton JOB sans problème.
27/03/2016, 13h02
FryHandiz

Bonjour,
J'ai bien peur que vous alliez un peu trop loin pour moi... En effet, j'ai pas vraiment besoin de parler de kétone ou d'amine car je m'intéresse pour le moment qu'au séquence nucléotidique et pas aux acides aminés. Donc du coup, je suis partie vers quelque chose de plus simple.
Je suis partie du principe que A et T sont complémentaires et G et C le sont aussi.
Du coup, pour coder sur deux bits chaque nucléotide, j'ai dit que A = 00, T=11, G=10 et C=01. Comme ça pour le complémentaire, pas de souci.
Mon projet est à rendre bientôt et je n'aurai jamais le temps de le terminer, mes connaissances sont vraiment très limitées xD Ce devait en effet être pour de l'alignement de séquence, mapper les reads, etc, mais vu que je ne pourrais jamais y arriver je vais me contenter de petites choses x)
En tout cas, merci beaucoup pour votre aide !
27/03/2016, 23h53
VivienD

Citation:

Envoyé par FryHandiz

Bonjour,
J'ai bien peur que vous alliez un peu trop loin pour moi... En effet, j'ai pas vraiment besoin de parler de kétone ou d'amine car je m'intéresse pour le moment qu'au séquence nucléotidique et pas aux acides aminés. Donc du coup, je suis partie vers quelque chose de plus simple.
[...]

C'est vrai que nous nous sommes "un peu" emportés et que nous sommes bien plus loin que demandé en nous adonnant au débat ci-dessus. :aie:

Citation:

Envoyé par FryHandiz

[...]
Je suis partie du principe que A et T sont complémentaires et G et C le sont aussi.
Du coup, pour coder sur deux bits chaque nucléotide, j'ai dit que A = 00, T=11, G=10 et C=01. Comme ça pour le complémentaire, pas de souci.
[...]

Si tu te contentes de la complémentarité des nucléotides seule, cet encodage convient tout à fait. :)

Citation:

Envoyé par FryHandiz

[...]
Mon projet est à rendre bientôt et je n'aurai jamais le temps de le terminer, mes connaissances sont vraiment très limitées xD Ce devait en effet être pour de l'alignement de séquence, mapper les reads, etc, mais vu que je ne pourrais jamais y arriver je vais me contenter de petites choses x)
[...]

C'est vraiment dommage, ça... :? Toutefois, j'espère que ça ne te décourageras pas; dis-toi que ça n'est qu'à charge de revanche. ;)

Citation:

Envoyé par FryHandiz

[...]
En tout cas, merci beaucoup pour votre aide !

Mais, de rien. ;)
Après, je me permets de te conseiller de conserver ce fil de discussion sous le coude (en le mettant dans tes marque-pages, par exemple) car dans notre "délire" nous avons abordé différentes méthodes concernant la représentation de données via des champs de bits; c'est un sujet auquel tout développeur a souvent affaire, du moins plus souvent qu'on ne le croit. :mrgreen:

PS: Personnellement, je vais me garder ce fil dans mes marque-pages. En effet, ça serait regrettable, outrageant voire blasphématoire que de laisser autant de verve et de dévotion sombrer dans les profondeurs abyssales et pleines de microbes de l'oubli. :lun:
Tiens! J'en avais encore à revendre! :weird:

Voir 40 message(s) de cette discussion en une page