Almacenamiento de datos de ADN de alta capacidad con Oligonucleótidos de longitud variable utilizando código de acumulación repetida y mapeo híbrido
Un práctico sistema de almacenamiento de datos de ADN con alta capacidad
Comenzamos con la construcción de una arquitectura de almacenamiento de datos y recuperación de datos de un almacenamiento basado en ADN (Fig. 1 A)). Los datos de usuario se segmentaron primero en 11.400 paquetes de usuario binarios con una longitud de paquete de 266 bits por paquete. Para corregir los errores que se producen en cualquier etapa de los procesos de almacenamiento de ADN, incluida la síntesis, la amplificación, el almacenamiento y la preparación de muestras para la secuenciación, aplicamos una codificación RA en paquetes de usuario binarios donde se generaron paquetes redundantes/de paridad del 5%. Con cada uno de los 12.000 paquetes binarios, se agregaron 14 bits para la indexación para ordenar los oligos estocásticos y 20 bits para la Comprobación de Redundancia Cíclica (CRC) para detectar los errores internos en cada paquete. Como resultado, el número total de bits asociados con cada paquete se convirtió en 300 bits (Ver archivo adicional 1: Figura S4). Después, mapeamos todas las secuencias binarias en secuencias de ADN a través del esquema de mapeo híbrido propuesto. Luego las secuencias de ADN fueron enviadas a Twist Bioscience para la síntesis de oligos. Después de recibir el pool de oligos sintetizados, lo amplificamos utilizando la Reacción en Cadena de la Polimerasa (PCR) antes de enviar las muestras a NovogeneAIT para su secuenciación utilizando Illumina HiSeq. En la última etapa, analizamos y decodificamos los datos de secuenciación para convertir los registros de ADN a datos binarios digitales. En primer lugar, muestreamos las lecturas de secuencia de millones del resultado de la secuenciación y realizamos el reverso de la codificación y el mapeo de RA para reconstruir los datos originales del usuario sin errores, validando la viabilidad de nuestro método.
Además de la recuperación completa de los datos utilizando los resultados de secuenciación, también analizamos cuantitativamente el esquema de almacenamiento basado en ADN propuesto y lo comparamos con otros esquemas de última generación, haciendo referencia a una tabla de comparación anterior (Fig. 1 C)). La definición detallada de las métricas de rendimiento de la tabla se describe en el archivo adicional 1: Sección S7. En la tabla, solo comparamos con los esquemas que se diseñaron y probaron con la premisa del formato de almacenamiento oligo pool, donde se sintetizaron los oligos cortos monocatenarios de longitud alrededor de 200 nt. Tenga en cuenta que con la suposición equivalente de almacenar hebras de ADN mucho más largas como , por ejemplo, 1000bp, el esquema de codificación propuesto sigue siendo factible, y la densidad de información neta aumentará con la longitud, logrando una densidad más alta que , por ejemplo, 1.84 bits/base sobre 1.74 bits/base (ver archivo adicional 1: Sección S3).
La alta densidad de información neta de 1.67 bits/nt obtenidos por el esquema de almacenamiento basado en ADN propuesto (Fig. 1 (D)) se debe principalmente a las siguientes dos técnicas que hemos utilizado. En primer lugar, el esquema de mapeo híbrido propuesto exhibe un potencial de mapeo de 1,98 bits/nt con una pequeña brecha del 1% desde el límite superior teórico de 2 bits/nt. En segundo lugar, el código RA optimizado para el control de errores tiene una pequeña redundancia de 1.05. Junto con la indexación de 14 bits y el CRC de 20 bits, el esquema obtiene una densidad de información neta de 1,67 bits/nt, lo que produce el 91% de la capacidad de Shannon (1,83 bits/nt con 0.tasa de abandono escolar del 5%), que es un 6% más que la última tasa más alta reportada (Archivo adicional 1: Sección S3). Teóricamente, en comparación con, el aumento en nuestra densidad de información es el resultado combinado de los oligoelementos de ADN de longitud variable ligeramente más largos (151nt-159nt versus 152nt, excluyendo los sitios de unión de imprimación), la menor redundancia de control de errores (1,05 versus 1,07) y la indexación más corta (14 bits versus 32 bits). La longitud de los oligoelementos de ADN está elaboradamente diseñada para hacer un uso completo de las técnicas de síntesis de ADN ampliamente disponibles (TWIST Bioscience, EE.UU.), que pueden sintetizar eficientemente oligoelementos de 200 nt de largo. El diseño de código RA optimizado ofrece una redundancia de control de errores ligeramente reducida con el supuesto equivalente de abordar una tasa de abandono práctica del 1,3%, mientras que la recuperación completa con cobertura de 10x (10,5 x pulg.) indica que se mantiene la resistencia a los errores. La diferencia más distintiva surge en la indexación, en la que usamos 14 bits únicamente para indicar el orden de 12000 oligos codificados, mientras que usamos 32 bits para representar las semillas necesarias para la transformación de Luby, que establece la base del código fuente, lo que resulta en bits de indexación redundantes.
Para verificar aún más que el rendimiento de alta capacidad del esquema de codificación propuesto se mantiene bien con el aumento del tamaño de los datos (escalabilidad), estimamos la densidad de información neta para codificar el tamaño de los datos con magnitudes más altas in silico, es decir, de 2 MB a 2000 MB. Las densidades estimadas disminuyen ligeramente con los aumentos exponenciales del tamaño de los datos debido al incremento de la longitud de indexación requerida para registrar un tamaño de datos más grande (Archivo adicional 1: Sección S3 y Fig. 1 E)). Se obtiene una densidad de 1,66 bits/nt para almacenar 2 MB de datos de origen, que sigue siendo un 6% mayor que . Además, tanto el código RA como la estrategia de mapeo híbrido que consiste en el esquema de codificación propuesto tienen una baja complejidad que son eficientes de implementar en la práctica. En particular, el uso de código RA evita el posible fallo de decodificación (debido a la pérdida de entradas iniciales para iniciar la decodificación en el proceso de selección) y la redundancia de direcciones que puede surgir en la fuente de ADN, y el mapeo híbrido logra un potencial de mapeo muy alto que es competitivo con la fuente de ADN, al tiempo que evita la alta complejidad que se exhibe en los códigos de bloques restringidos convencionales.
Además, calculamos computacionalmente la densidad física que el esquema propuesto podría exhibir. A través de experimentos de dilución, los autores observaron una tasa de abandono del 4% con una muestra de almacenamiento de ADN de 10 pg, que casi se acercó a su límite de decodificación (que estaba predeterminado por la redundancia de código). El código RA utilizado en nuestro esquema fue diseñado de manera óptima con un nivel de redundancia bajo el mismo supuesto de tasa de abandono considerado . También hemos demostrado que teóricamente nuestro código puede tolerar una tasa de abandono de hasta el 4,75% (archivo adicional 1: Figura S4), que está por encima de la tasa de abandono del 4% observada en la secuenciación de la muestra 10pg. Con un límite de decodificación similar, nuestro esquema propuesto probablemente funcionaría igual que la fuente de ADN en los experimentos de bajo peso molecular (por ejemplo, con muestras de 10 pg) debido al uso de las mismas tuberías, protocolos y estándares de experimento. En otras palabras, el diseño del código en la etapa inicial permite que el sistema propuesto pueda recuperar datos de condiciones propensas a errores en los experimentos de dilución similares a la fuente de ADN. Bajo la suposición de molecules 1300 moléculas por oligo en promedio, una profundidad de secuenciación de 511x y tuberías, protocolos y estándares equivalentes como el experimento de dilución de 10 pg en fuente de ADN, podríamos estimar computacionalmente que nuestro esquema alcanzará una densidad física de 239 PB/g \(\left (\frac {266 * 11400/8\text {byte}}{1300*11400*1.0688*10^{-19}\texto {gramo}}\derecha)\). Sin embargo, se requiere un experimento riguroso para verificar esta densidad física calculada computacionalmente.
Diseño de código RA y esquema de mapeo híbrido para almacenamiento de ADN
Diseñamos un método de codificación que comprende código RA (repetición de nivel oligo) y un esquema de mapeo híbrido eficiente.
Diseño de código RA
En los sistemas de comunicación tradicionales, el código RA se utiliza a nivel de bits, donde se generan bits redundantes para mitigar los errores de sustitución. Sin embargo, el almacenamiento de ADN es propenso no solo a errores de sustitución, sino también a errores de inserción y eliminación. Por lo tanto, en lugar de la codificación RA convencional a nivel de bits, diseñamos una codificación RA a nivel de paquete para el almacenamiento de ADN de modo que un paquete sujeto a errores de inserción, eliminación o sustitución pudiera recuperarse a través del decodificador RA. Como se describió anteriormente, hemos segmentado un archivo digital grande en paquetes más pequeños del mismo tamaño. Estos paquetes se consideraron los paquetes fuente que se utilizaron para generar los paquetes redundantes o de paridad utilizando código RA sistemático Fig. 2 A). Tenga en cuenta que cada paquete se incorporó con CRC para detectar errores en el paquete. Para los paquetes que pasaron la prueba CRC en el decodificador, los consideramos recuperados correctamente, mientras que los demás se consideraron eliminados o borrados. Por lo tanto, el problema general de diseño de código para el almacenamiento de ADN se convirtió en el diseño de código para el canal de borrado. Para garantizar una alta confiabilidad, el diseño del código se realizó considerando una probabilidad de abandono ligeramente mayor que la probabilidad real de abandono. En este trabajo, se consideró que la tasa real de abandono escolar era del 1,3%, que se informó en el periódico fountain . Por lo tanto, diseñamos el código RA de manera que el código resultante exhibiera un umbral asintótico más alto que la probabilidad de abandono de 0,013. Siguiendo el procedimiento de optimización (ver archivo Adicional 1: Sección S2), diseñamos un código RA de tasa 0.95, que da un umbral asintótico de 0.0475. El código resultante muestra solo una brecha de 0.0025 desde el límite de capacidad de Shannon (0.05). El rendimiento de corrección de errores simulado del código RA diseñado se muestra en el archivo adicional 1: Figura S4. Debido a la tarifa 0.Código RA de 95, generamos 600 paquetes redundantes / de paridad basados en 11.400 paquetes de origen, recibiendo un total de 12.000 paquetes binarios después de la codificación.
Esquema de mapeo híbrido
A continuación, consideramos representar los datos digitales en el contexto de ADN que denominamos mapeo de ADN. Una estrategia de mapeo de ADN debe permitir que las secuencias de oligo mapeadas satisfagan las restricciones bioquímicas, brindando así estabilidad al almacenamiento. Hay dos restricciones en los datos de ADN como las siguientes: (i) El contenido de GC (la relación del número total de ‘G’ y ‘C’ contra el número total de nucleótidos en una secuencia) debe ser cercano al 50% (ii) Todas las longitudes de carrera de homopolímeros (la longitud de nucleótidos repetitivamente consecutivos) debe ser inferior a 4 . Tenga en cuenta que el mapeo binario a cuaternario, es decir, el mapeo de dos bits a un nucleótido, que exhibe el potencial de mapeo óptimo (2 bits/nt), no siempre cumple con los requisitos mencionados anteriormente. En su lugar, a menudo no cumple con la restricción de ejecución de homopolímero máximo. Las limitaciones existentes en el almacenamiento de datos de ADN reducen el potencial de mapeo efectivo, afectando negativamente la capacidad de almacenamiento de datos de ADN. Por lo tanto, exploramos el enfoque de diseñar código restringido con alta velocidad de código y desarrollamos una estrategia de mapeo híbrido para garantizar que las secuencias de oligo cumplan con las demandas bioquímicas con un sacrificio mínimo del potencial de mapeo.
Este esquema de asignación consta de dos métodos de asignación diferentes, a saber, la asignación intercalada y la asignación VLC. El primero funciona como mapeo primario debido a su potencial de mapeo aproximadamente óptimo, es decir, 1.995 bits / nt y este último funciona como la copia de seguridad que entra en juego cuando la primera asignación no produce secuencias de ADN válidas (es decir, secuencias que satisfacen el contenido de GC y las restricciones de ejecución de homopolímeros). En el método de asignación posterior, se construye una tabla de búsqueda auxiliar con baja complejidad de codificación y decodificación. Mientras tanto, este método exhibe un potencial de mapeo de 1.976 bits/nt que es mucho más alto que los códigos de bloque con la complejidad equivalente. La combinación de estas dos estrategias de mapeo da como resultado un potencial de mapeo promedio de alrededor de 1,98 bits/nt con los datos estocásticos. En otras palabras, en el peor de los casos, donde todos los datos se codifican utilizando VLC, aún logramos una estimación de potencial de mapeo alto (1.976 bits/nt). Sin embargo, en el mejor de los casos, cuando todos los datos se mapean utilizando el mapeo intercalado, podríamos lograr un potencial muy alto de 1,995 bits/nt.
Los datos digitales primero pasan por el método de mapeo interleaved para generar las secuencias de ADN. En el método de mapeo intercalado, las secuencias binarias se mapean primero usando mapeo binario a cuaternario. Con el aumento de la longitud oligo, la restricción de contenido GC a menudo se satisface debido a la característica estocástica de los datos binarios. Sin embargo, esta asignación tiende a no satisfacer la restricción de ejecución de homopolímero. Para resolver este problema, introducimos un intercalador después del mapeo binario a cuaternario, que codifica el orden original de las secuencias de nucleótidos. Después del entrelazado, se realiza una prueba de cribado para verificar la ejecución de homopolímeros de la secuencia resultante. Si la secuencia resultante pasa la prueba, esa secuencia se considera una secuencia válida para la síntesis, de lo contrario, el entrelazado se realiza de nuevo en la secuencia original con un patrón de entrelazado diferente. En este trabajo, consideramos 4 patrones de interleaving predefinidos, donde se agrega un nucleótido de bandera (A/T/G/C) al final de la secuencia de ADN intercalado para indicar el patrón de interleaving (Archivo adicional 1: Sección S8). Tenga en cuenta que el nucleótido de bandera adjunto se incluye en la determinación de la ejecución de homopolímeros de la secuencia durante la prueba de detección. Solo usamos un nucleótido adicional (bandera) para mantener una alta densidad de información neta. En consecuencia, el número de ensayos de intercalado se limita a 4. Si la secuencia sigue sin satisfacer la demanda después del número máximo de ensayos, la secuencia se envía al método de mapeo VLC (Fig. 2 B) y expediente adicional 1: Sección S4).
El mapeo VLC está inspirado en la construcción de código de secuencia restringida de longitud variable (VLCS), comúnmente utilizado para codificar datos en códigos que satisfacen restricciones en sistemas restringidos, como los sistemas de grabación óptica donde surgen problemas de límite de longitud de ejecución y sin CC . En el escenario de almacenamiento de ADN donde existen restricciones similares, el código VLCS se puede modificar de manera efectiva a un método de mapeo. Tenga en cuenta que, a medida que usamos el código RA a nivel de paquete para el control de errores, la propagación de errores liderada por el código VLCS está limitada en un paquete y no tiene influencia en la tasa de abandono general de las secuencias codificadas.
Generamos esta regla de asignación en las siguientes cuatro etapas. En primer lugar, considerando la restricción de las corridas máximas de homopolímeros, el almacenamiento basado en ADN se vio como un sistema restringido con límite de longitud de corrida (RLL), denotado por (M,d, k), donde M=4, d=0 y k=2 (Archivo adicional 1: Sección S5). En consecuencia, se generó el diagrama de transición de estados finitos (FSTD) del almacenamiento de datos de ADN restringido por homopolímeros (4,0,2) (Archivo adicional 1: Sección S5 y Fig. 2 C, i)). En la segunda etapa, basándonos en el FSTD generado, deducimos que la capacidad de almacenamiento de ADN restringido por homopolímeros (4, 0, 2) es de 1,982 bits/nt (Archivo adicional 1: Sección S5). También establecimos un conjunto mínimo completo (un conjunto finito de palabras cuyas concatenaciones incluyen todas las posibles secuencias que satisfacen las restricciones), donde enumeramos todas las palabras que se originan y terminan en el estado s0 en la Fig. 2 C, i). Como resultado. obtuvimos un conjunto mínimo {1,2,3,01,02,03,001,002,003}, en el que todos los elementos cumplen con las restricciones y no tienen prefijos. Estas dos propiedades garantizan que cualquier concatenación de los elementos de este conjunto produzca secuencias que satisfagan las restricciones y que sean posibles palabras de código de transición para el sistema restringido. Tenga en cuenta que el conjunto de palabras clave de transición resultante se relaciona con la profundidad y el ancho de la concatenación. Para reducir la complejidad de la codificación, utilizamos directamente el conjunto mínimo completo como conjunto de palabras de código de transición.
En la tercera etapa, utilizamos el árbol de codificación Huffman para generar una asignación óptima del conjunto de palabras fuente binarias de longitud variable al conjunto de palabras de código de transición mencionado anteriormente (Fig. 2 C, ii)). Esta asignación óptima de uno a uno dio una tasa de código promedio de 1.976 bits/nt (Fig. 2 (C, iii) y véase el expediente adicional 1: Sección S5). Mientras tanto, la eficiencia de este mapeo se acerca a \(\sigma = \frac {1.976}{1.982}=99.7\%\), presentando solo un 0,3% de diferencia con respecto a la capacidad del sistema restringido (4,0,2). En términos de potencial de mapeo, este mapeo supera al código restringido por bloques propuesto, en el que se construyó un código restringido (4,0,2) utilizando bloques de ADN de 39 nt como palabras de código, logrando un potencial de mapeo de 1,95 bits/nt. Además, el código de bloque 39nt tampoco es práctico para el almacenamiento de datos de ADN tradicional, donde se consideran secuencias de ADN mucho más largas (palabras de código), es decir, 200nt. En contraste, el enfoque de mapeo de longitud variable tiene una baja complejidad de codificación, independientemente de la longitud total de las secuencias oligo resultantes.
En la última etapa, después de asignar las palabras de origen a las palabras de código de transición en sucesión contra cada secuencia binaria, realizamos la precodificación en las secuencias cuaternarias codificadas de acuerdo con la función de cambio de estado yj=yj-1+xj(mod M), donde yj es el símbolo de precodificación de salida actual, yj-1 es el último símbolo de precodificación de salida, xj es el símbolo de entrada actual, M es el tamaño del alfabeto del sistema. Esta precodificación transferirá el código restringido codificado (M,d,k) al código RLL (M,d+1,k+1). Luego convertimos los símbolos cuaternarios de {0,1,2,3} a {‘A’, ‘T’, ‘C’, ‘G’} y obtuvimos las secuencias oligo finales que satisfacen la restricción de no homopolímero y corren más de 3nt. Un ejemplo de esta estrategia de asignación se puede encontrar en el archivo adicional 1: Sección S6.
A través del esquema de mapeo híbrido, generamos 12,000 secuencias de ADN con una distribución de longitud que varía de 150 a 159 nt (excluyendo 40 nt de sitios de imprimación) para el flujo de datos binarios (Fig. 2 E)). Específicamente, la longitud de las secuencias que se mapearon a través del mapeo intercalado se convirtió en 151nt, mientras que la longitud de las secuencias que se mapearon a través del mapeo VLC osciló entre 150, 152 y 159nt. Tenga en cuenta que no hubo una secuencia con una longitud de 151nt que se originó a partir de la asignación de VLC, ya que se agregó un nucleótido para que esta secuencia mapeada de 151nt fuera de 152nt (Fig. 2 C, iv)). El nucleótido añadido era para distinguir entre los métodos de mapeo. Esto permite el uso de la des-asignación correcta durante la recuperación de los datos almacenados en el decodificador.
Para recuperar datos, las secuencias preparadas del proceso de secuenciación se envían al decodificador para recuperar los datos del usuario (Fig. 2 D)). El decodificador distingue primero el método de asignación. Si la longitud de la secuencia recibida es de 151 nt, el decodificador aplica el reverso de la asignación intercalada basada en el nucleótido de bandera y la regla de asignación binario a cuaternario. De lo contrario, el decodificador aplica el reverso de la asignación VLC donde se realiza el reverso de la precodificación y la asignación. Después de eso, cada secuencia binaria invertida se considera correcta o borrada según la comprobación CRC. Finalmente, con un algoritmo de paso de mensajes, el decodificador RA recupera todos los paquetes de secuencia borrados en función de las conexiones entre paquetes.
Resultados de secuenciación y análisis de recuperación de datos
Después de secuenciar el grupo de oligos sintetizados, recibimos más de 10 millones de lecturas de secuencias sin procesar en un tamaño total de 3,2 Gigabytes de NovogeneAIT. Estas secuencias incluyen lecturas ruidosas generadas durante la secuenciación. En base a los resultados de la secuenciación, en primer lugar, analizamos la fiabilidad de los datos de la secuenciación en términos de examen de la calidad de los datos, distribución del contenido A/T/G/C y distribución de la tasa de error. A partir del resultado del análisis de errores, estudiamos la fiabilidad de nuestro esquema de decodificación para recuperar los datos codificados con diferentes coberturas de muestra.
Resultados de secuenciación
Analizamos el valor de calidad para cada posición de base a lo largo de las lecturas secuenciadas para evaluar la calidad de los datos. La puntuación de calidad es una estimación de la fiabilidad de las lecturas secuenciadas que se relaciona con la tasa de error de cada posición base. Se calcula por Q = – 10log10e, donde e es la tasa de error de la posición base . Las puntuaciones de calidad de cada base de las lecturas de secuenciación varían de 30 a 40 (Fig. 3 A)), de gran calidad. Además, observamos que la tasa de error aumenta con la extensión de las lecturas secuenciadas, mientras que con una tasa promedio de 0,015% en cada base a lo largo de las lecturas (Fig. 3 B)). Esto es probablemente debido al consumo de reactivo de secuenciación, que es un fenómeno común en la plataforma de secuenciación de alto rendimiento Illumina que se basa en la tecnología de secuenciación por síntesis (SBS). Como era de esperar, las primeras bases tienen una tasa de error de secuenciación más alta que otras. Esto podría deberse al enfoque del elemento de detección del sensor de imagen de fluorescencia del secuenciador, que puede no ser lo suficientemente sensible al comienzo de la secuenciación. Como resultado, la calidad de la lectura de fluorescencia adquirida es baja. Recuerde que las secuencias se adjuntaron con un par de sitios de unión de imprimación 20nt en ambos extremos y, por lo tanto, las primeras bases propensas a errores (alrededor de 6nt) no tienen influencia en la decodificación, ya que la prueba CRC y la codificación/decodificación RA se diseñaron excluyendo los sitios de unión. En otras palabras, una secuencia será identificada como borrada por el decodificador CRC debido a los errores en otras posiciones (fuera de los cebadores).