Información

Cómo seleccionar genes antes de la relación log2 en una matriz de expresión génica RNASeq, según la mediana de la señal


Quiero transformar una matriz de expresión de ARNm TCGA (en formato de datos lineales) para registrar2-ratios y luego ejecute una selección de características (gen), seleccionando los 1000 genes más variantes (genes con una desviación estándar más alta entre las muestras). El flujo de trabajo es el siguiente:

  1. Seleccione genes "buenos" antes del registro2 proporción (genes cada uno con una señal mediana al menos t en pag% de muestras);
  2. En genes seleccionados, ejecute log2 proporción, dividiendo cada gen por su señal mediana y luego log2-transformar la matriz de resultados;
  3. Seleccione los 1000 genes más variantes de todas las muestras.

Como selecciono t y pag?


No hay una regla para arreglar t y pag. Depende del nivel de rigor que espere. Valor de t depende de lo que se considere una concentración activa; esto no tiene por qué ser el mismo para todos los genes.

Este es un dato de RNAseq; No entiendo cuál es la señal "mediana" de la que estás hablando. Para cada muestra, un gen tendría un valor de expresión normalizado que normalmente es RPKM (lecturas por kilobase por millón de lecturas mapeadas). Si tiene réplicas para cada muestra, tome la media, no la mediana.

Con respecto al cálculo de las relaciones logarítmicas: Siempre tenga cuidado con esto, especialmente en caso de ceros. En lugar de las relaciones de registro, puede utilizar algún tipo de métrica de ganancia:

si relación = x / y entonces ganancia = (x-y) / y

También puede hacer un análisis de componentes principales en los datos y seleccionar primeronortecomponentes principales.


Cómo seleccionar genes antes de la relación log2 en una matriz de expresión génica RNASeq, según la mediana de la señal - Biología

No existe una solución general para seleccionar "t" y "p". Tales elecciones son en gran parte arbitrarias. Además, para una plataforma de matriz, si se supone que "t" tiene algo que ver con "expresado", el valor de "t" será diferente para cada sonda en la matriz.

Dado que, en última instancia, va a filtrar en función de la varianza, le sugiero que comience con sus datos transformados logarítmicamente centrados en la mediana y simplemente elija los 1000 genes más variables.

Los datos a los que se refiere OP son RNAseq, por lo que no hay sondas. Se puede realizar una corrección de sesgo de secuenciación para ellos.

Gracias Sean. Con respecto a su última sugerencia, pensé que podría presentar problemas, ya que los genes con una señal mediana muy baja pueden mostrar una alta varianza cuando se transforma la logratio. ¿Qué piensas?

La transformación del registro de pozo reducirá la varianza. Imagínese 4 muestras [0,5, 2, 8, 32]. Sin la transformación logarítmica, la varianza es 213,5625, pero cuando transforma los datos en log2, la varianza se reduce a 6,67

En cualquier caso, si la expresión es constantemente baja, la varianza será baja. Debe tener cuidado con las transformaciones logarítmicas, especialmente cuando se realizan estudios de expresión diferencial. Sugeriría que hagas la transformación logarítmica después de seleccionar la mediana y la varianza.

Seguí su primer consejo utilizando umbrales personalizados.

Con la esperanza de que esto sea útil, el código de la canalización está disponible en https://gist.github.com/fbrundu/adf45f7ce817572bba10

Filtre los genes que estaban por debajo del percentil 5 en general en más del 5% de las muestras. Creo que podría ser un umbral razonable. Por si acaso dime.


Archivos de datos GSEA

¿Cómo creo un archivo de conjunto de datos de expresión? ¿Qué tipos de datos de expresión puedo analizar?

GSEA requiere que los datos de expresión estén en un archivo RES, GCT, PCL o TXT. Los cuatro formatos de archivo son archivos de texto delimitados por tabulaciones. Para obtener detalles de cada formato de archivo, consulte Formatos de datos.

GenePattern proporciona varios módulos para convertir datos de expresión en archivos gct y / o res:

ExpressionFileCreator convierte datos de expresión sin procesar de archivos Affymetrix CEL.

GEOImporter y caArrayImportViewer crean un archivo GCT basado en datos de expresión extraídos del repositorio de datos de expresión de microarrays de GEO o caArray, respectivamente.

El módulo MAGEImportViewer convierte datos en formato MAGE-ML. MAGE-ML es el formato estándar para almacenar datos de microarrays de ADNc y Affymetrix en el repositorio de ArrayExpress.

Para usar datos de expresión almacenados en cualquier otro formato (como datos de microarrays de ADNc), primero convierta los datos en un archivo de texto delimitado por tabuladores que contenga mediciones de expresión con genes como filas y muestras como columnas y luego modifique ese archivo de texto para cumplir con la gct requisitos de formato de archivo como se describe en Expression Datasets en el Guía del usuario de GSEA.

Si utiliza datos de relación de dos colores, consulte también Datos de microarrays de ADNc.

Errores de análisis: Si ve el siguiente error de análisis cuando carga su archivo de datos, verifique la extensión del archivo:
Hubo errores: ERRORES #: 1 Problemas de reproducción…

La extensión del archivo del conjunto de datos de expresión identifica el formato del archivo. Si un archivo gct, res o pcl tiene una extensión de archivo .txt, verá el error de análisis cuando cargue el archivo en GSEA. Compruebe que la extensión del archivo coincida con el formato del archivo. Tenga en cuenta que algunos sistemas operativos (como Windows) se pueden configurar para ocultar extensiones de archivo conocidas. Si su sistema operativo está configurado para ocultar extensiones conocidas, un archivo llamado test.gct.txt será listado como test.gct. Mire el tipo de archivo del archivo: debe ser GCT (o RES o PCL), no Documento de texto.

¿Cómo filtro o preproceso mi conjunto de datos para GSEA?

La forma de filtrar o preprocesar sus datos depende de su estudio. Aquí hay algunas pautas a considerar:

  • Identificadores de sonda versus identificadores de genes. Normalmente, su conjunto de datos contiene los identificadores de sonda nativos de su chip de ADN de plataforma de microarrays. GSEA puede analizar los identificadores de sondas o colapsar cada conjunto de sondas en un vector genético, donde el gen se identifica mediante el símbolo del gen. El colapso de los conjuntos de sondas evita que varias sondas por gen inflen las puntuaciones de enriquecimiento y facilita la interpretación biológica de los resultados del análisis.
  • Filtros de llamadas AP. Puede ejecutar GSEA en datos filtrados o sin filtrar. Normalmente, el equipo de GSEA ejecuta el análisis en datos sin filtrar. Un enfoque sugerido es ejecutar GSEA en los datos sin filtrar. Si los resultados parecen dominados por conjuntos de genes con genes pobremente expresados, puede obtener información sobre qué umbrales utilizar para los filtros de llamadas.
  • Valores de expresión. El algoritmo GSEA examina las diferencias en los valores de expresión en lugar de los valores en sí. Por ejemplo, puede tener datos de escala natural o niveles de expresión registrados, puede tener datos Affymetrix o datos de proporción de dos colores. & Lta ​​name = "_ Toc120959112" & gt & lt / a & gt Como en la mayoría de las metodologías de análisis de datos, los mismos datos de expresión representados en diferentes formatos pueden generar diferentes resultados de análisis. Se esperan las diferencias. GSEA no puede determinar qué resultados son & quotcorrectos & quot & lta name = "_ Toc120959112" & gt & lt / a & gt

Para obtener más información, consulte Preparación de archivos de datos en la Guía del usuario de GSEA.

¿Debo usar datos de escala natural o logarítmica para GSEA?

Recomendamos utilizar datos a escala natural. Lo usamos cuando calibramos el método GSEA y parece funcionar bien en casos generales.

Las técnicas de modelado tradicionales, como la agrupación en clústeres, a menudo se benefician del preprocesamiento de datos. Por ejemplo, se pueden filtrar los datos de expresión para eliminar los genes que tienen baja varianza en el conjunto de datos y / o transformar los datos en logaritmos para hacer que la distribución sea más simétrica. El algoritmo GSEA no se beneficia de tal procesamiento previo de los datos.

¿Cuántas muestras necesito para GSEA?

Esto depende de su problema específico y de las características de los datos; sin embargo, como regla general, normalmente desea analizar al menos diez muestras.

Si tiene réplicas técnicas, generalmente desea eliminarlas promediando o alguna otra técnica de reducción de datos. Por ejemplo, suponga que tiene cinco muestras de tumores y cinco muestras de control cada una de las cuales se procesan tres veces (tres columnas de réplica) para un total de 30 columnas de datos. Promediaría las tres columnas replicadas para cada muestra y crearía un conjunto de datos que contenga 10 columnas de datos (cinco de tumor y cinco de control).

¿Cómo creo un archivo de etiqueta de fenotipo? ¿Qué tipo de experimentos puedo analizar?

GSEA se puede utilizar para analizar experimentos de cualquier tipo (incluidas series de tiempo, tres o más clases, etc.). El archivo ASCII de etiquetas de fenotipo (cls) define los fenotipos experimentales y asocia cada muestra de tu conjunto de datos con uno de esos fenotipos. El archivo cls es un archivo delimitado por tabulaciones ASCII, que puede crear fácilmente con un editor de texto. Para obtener más información, consulte Preparación de archivos de datos en la Guía del usuario de GSEA.

¿Qué conjuntos de genes están disponibles? ¿Puedo crear mis propios conjuntos de genes?

Puede utilizar los conjuntos de genes en la base de datos de firmas moleculares (MSigDB) o crear el suyo propio. Para obtener más información sobre los conjuntos de genes MSigDB, consulte la página MSigDB de este sitio web. Para obtener más información sobre la creación de conjuntos de genes o el uso de conjuntos de genes con GSEA, consulte Preparación de archivos de datos en el Guía del usuario de GSEA.

¿Cuántos genes debería haber en un conjunto de genes?

GSEA ajusta automáticamente las estadísticas de enriquecimiento para tener en cuenta los diferentes tamaños de conjuntos de genes, como se describe en el Información suplementaria para el documento GSEA 2005 PNAS.

¿Puede GSEA analizar un conjunto de genes que contiene genes duplicados? conjuntos de genes duplicados?

Los genes duplicados en un conjunto de genes y los conjuntos de genes duplicados afectan los resultados de GSEA. GSEA elimina automáticamente los genes duplicados de cada conjunto de genes, pero no busca conjuntos de genes duplicados. Para obtener más información, consulte Conjuntos de genes en el Guía del usuario de GSEA.

¿Puede GSEA analizar un conjunto de genes que contiene genes que no están en mi conjunto de datos de expresión?

El análisis de enriquecimiento de conjuntos de genes restringe automáticamente los conjuntos de genes a los genes en el conjunto de datos de expresión. El informe de análisis enumera los conjuntos de genes y el número de genes que se incluyeron y excluyeron del análisis.

¿Qué plataformas de matriz y especies de organismos admite GSEA?

GSEA trabaja en alguna datos, siempre que los identificadores de genes en sus datos de expresión coincidan con los del archivo de conjuntos de genes.

Normalmente, GSEA utiliza conjuntos de genes de MSigDB. Todos los conjuntos de genes en MSigDB consisten en símbolos de genes humanos. GSEA tiene herramientas integradas para la conversión entre una variedad de otros identificadores de genes en símbolos de genes humanos por medio de archivos CHIP con formato especial. Los archivos CHIP proporcionan el mapeo entre los identificadores de genes en sus datos de expresión y los identificadores de genes en los conjuntos de genes. Específicamente, nuestros archivos CHIP proporcionan las asignaciones de todo tipo de plataformas diferentes (por ejemplo, ID de conjuntos de sondas de Affymetrix de ratón, ID de conjuntos de sondas de Affymetrix humanas, etc.) a símbolos de genes humanos.

Si sus datos se generaron a partir de muestras no humanas, debe decidir si el uso de MSigDB satisface sus necesidades. Las opciones son:

  1. La especie no humana sirve como modelo para estudiar las condiciones relevantes para la biología humana. En este caso, desea conjuntos de genes que se conserven entre los humanos y su organismo modelo. Entonces, MSigDB es la opción correcta y solo necesitará proporcionar el archivo CHIP apropiado para el análisis.
  2. La especie no humana es el tema de su investigación y no tiene planes de compararla con conjuntos de genes humanos. En este caso, aún puede usar MSigDB si su organismo se encuentra entre las fuentes de algunos conjuntos de genes de MSigDB (por ejemplo, ratón o rata) y solo necesitará proporcionar el archivo CHIP apropiado para el análisis.
  3. La especie no humana es el tema de su investigación y no desea utilizar conjuntos de genes MSigDB por otras razones. En este caso, debe proporcionar su propia base de datos de conjuntos de genes como un archivo GMT o GMX. Los formatos de archivo se describen aquí. Por supuesto, aún debe asegurarse de que los identificadores de genes en sus datos coincidan con los de su base de datos de conjuntos de genes. Si los identificadores no coinciden entre sí, también debe proporcionar un archivo CHIP con las asignaciones adecuadas. El formato de archivo CHIP se describe aquí.

Para ver qué archivos CHIP están disponibles en nuestra distribución (nota: nuestros archivos CHIP solo proporcionan asignaciones a símbolos de genes humanos): inicie la aplicación de escritorio GSEA y haga clic en [. ] en "Plataforma (s) de chip" en la página "Ejecutar GSEA".

Si su plataforma no está en esta lista, tiene las siguientes opciones:

  1. Cree su propio archivo CHIP para asignar los identificadores de genes específicos de su plataforma a los símbolos de genes humanos y luego use su archivo CHIP para colapsar el conjunto de datos en GSEA. El formato de archivo CHIP se describe aquí.
  2. Convierta los identificadores de su plataforma en símbolos de genes humanos fuera de GSEA, luego ejecute GSEA con 'Collapse dataset' = FALSE.

Asegúrese de que los símbolos genéticos en el conjunto de datos contraído aparezcan solo una vez. El simple hecho de reemplazar los identificadores con símbolos de genes humanos generalmente no es suficiente porque algunos de los identificadores pueden corresponder a los mismos símbolos de genes humanos, lo que da como resultado filas duplicadas con diferentes valores de expresión. En este caso, GSEA elegirá arbitrariamente una de las filas con los mismos símbolos genéticos para el análisis, lo que no recomendamos.

¿Puede GSEA analizar los datos de expresión de miARN?

La única forma de que GSEA analice los datos de expresión con identificadores de miARN es proporcionar conjuntos de genes hechos de identificadores de miARN coincidentes. Esto no es posible con los conjuntos de genes MSigDB, que consisten predominantemente en genes que codifican proteínas en forma de símbolos de genes humanos.


El análisis de secuencia de ARN de los patrones de expresión de genes espaciotemporales durante el desarrollo de la fruta reveló genes de referencia para la normalización de la transcripción en ciruelas

El análisis transcripcional que descubre las redes reguladoras de genes relacionados con la maduración de la fruta es cada vez más importante para maximizar la calidad y minimizar las pérdidas de frutas económicamente importantes como las ciruelas. La secuenciación de ARN (RNA-Seq) y la reacción en cadena de la polimerasa con transcripción inversa cuantitativa en tiempo real (qRT-PCR) son herramientas importantes para realizar transcriptómica de alto rendimiento. El éxito de la transcriptómica depende de las transcripciones de alta calidad de frutos de ciruela enriquecidos con polifenólicos y polisacáridos, mientras que la fiabilidad de los datos de cuantificación se basa en una normalización precisa utilizando genes de referencia adecuados. Optimizamos un procedimiento para el aislamiento de ARN de alta calidad a partir de tejidos vegetativos y reproductivos de cultivares de ciruela climatéricos y no climatéricos y realizamos transcriptómica de alto rendimiento. Identificamos 20 genes de referencia candidatos a partir de transcripciones significativamente no expresadas diferencialmente de datos de RNA-Seq y verificamos su estabilidad de expresión utilizando qRT-PCR en un total de 141 muestras de ciruela que incluían tejidos de pulpa, piel y hojas de varios cultivares recolectados en tres ubicaciones. durante un período de 3 años. Análisis de estabilidad del ciclo umbral (C T) valores usando BestKeeper, delta (Δ) CT, Se revela el software NormFinder, geNorm y RefFinder S Y proteína de tráfico relacionada con proteínas (LUN), factor de alargamiento 1 alfa (EF1α), y factor de iniciación 5A (IF5A) como los mejores genes de referencia para la normalización precisa de la transcripción en diferentes muestras de tejido. Supervisamos los patrones de expresión espacio-temporal de las transcripciones expresadas diferencialmente durante el proceso de desarrollo después de la normalización precisa de los datos de qRT-PCR utilizando la combinación de los dos mejores genes de referencia. Este estudio también ofrece una guía para seleccionar los mejores genes de referencia para futuros estudios de expresión génica en otros cultivares de ciruela.

Esta es una vista previa del contenido de la suscripción, acceda a través de su institución.


Los genes que tienen recuentos muy bajos en todas las bibliotecas deben eliminarse antes del análisis posterior. Esto se justifica tanto por motivos biológicos como estadísticos. Desde el punto de vista biológico, un gen debe expresarse a un nivel mínimo antes de que sea probable que se traduzca en una proteína o se considere biológicamente importante. Desde un punto de vista estadístico, es muy poco probable que los genes con recuentos consistentemente bajos se evalúen como DE significativamente porque los recuentos bajos no proporcionan suficiente evidencia estadística para realizar un juicio confiable. Por lo tanto, estos genes pueden eliminarse del análisis sin pérdida de información.

Como regla general, requerimos que un gen tenga un recuento de al menos 10-15 en al menos algunas bibliotecas antes de que se considere que se expresa en el estudio. Podríamos seleccionar explícitamente genes que tengan al menos un par de recuentos de 10 o más, pero es un poco mejor basar el filtrado en valores de recuento por millón (CPM) para evitar favorecer genes que se expresan en bibliotecas más grandes. sobre los expresados ​​en bibliotecas más pequeñas. Para el análisis actual, mantenemos los genes que tienen valores de CPM por encima de 0,5 en al menos dos bibliotecas:

Aquí se ha elegido el límite de 0,5 para el CPM porque es aproximadamente igual a (10 ​​/ L ) donde (L ) es el tamaño mínimo de biblioteca en millones. El tamaño de las bibliotecas aquí es de 20 a 25 millones. Usamos un valor redondo de 0.5 solo por simplicidad, el valor exacto no es importante porque el análisis de expresión diferencial descendente no es sensible a los pequeños cambios en este parámetro. El requisito de las bibliotecas ( ge 2 ) es que cada grupo contiene dos réplicas. Esto asegura que se retendrá un gen si se expresa en ambas bibliotecas que pertenecen a cualquiera de los seis grupos.

La regla de filtrado anterior intenta mantener el número máximo de genes interesantes en el análisis, pero también son posibles otros criterios de filtrado sensibles. Por ejemplo, keep & lt- rowSums (y $ count) & gt 50 es un criterio muy simple que mantendría genes con un recuento total de lecturas de más de 50. Esto daría resultados posteriores similares para este conjunto de datos al filtrado realmente utilizado. Cualquiera que sea la regla de filtrado, debe ser independiente de la información del archivo de destinos. No debe hacer ninguna referencia a qué bibliotecas de ARN pertenecen a qué grupo, porque hacerlo sesgaría el análisis de expresión diferencial posterior.

El objeto DGEList está subconjunto para retener solo los genes no filtrados:

La opción keep.lib.sizes = FALSE hace que los tamaños de la biblioteca se vuelvan a calcular después del filtrado. Esto se recomienda generalmente, aunque el efecto en el análisis posterior suele ser pequeño.


Conclusiones

En estudios de simulación y datos reales, limma con el l, l2, y r2 transformaciones funcionaron mejor que limma con el voom transformación para datos con tamaño de muestra pequeño (nCases = nControls = 3) o grande (nCases = nControls = 100). Para un tamaño de muestra moderado (nCases = nControls = 30 o 50), limma con el rv y rv2 La transformación funcionó mejor que limma con el voom transformación. Esperamos que estas nuevas transformaciones de datos puedan proporcionar a los investigadores un análisis de expresión diferencial más potente utilizando datos de RNA-seq.


Inspección de los resultados del mapeo

El archivo BAM contiene información sobre dónde se asignan las lecturas en el genoma de referencia. Pero como es un archivo binario que contiene información para muchas lecturas (varios millones para estas muestras), es difícil inspeccionar y explorar el archivo.

Una poderosa herramienta para visualizar el contenido de los archivos BAM es el Visor de Genómica Integrada (IGV).

Hands_on Hands-on: Inspección de los resultados del mapeo

  1. Instale IGV (si aún no está instalado)
  2. Inicie IGV localmente
  3. Expanda el archivo param-archivo mapped.bam (salida de ESTRELLA DE ARN herramienta) para GSM461177
  4. Haga clic en el local en pantalla con IGV local D. melanogaster (dm6) para cargar las lecturas en el navegador IGV

Comentario Comentarios

Para que este paso funcione, necesitará tener instalado IGV o Java web start en su máquina. Sin embargo, las preguntas de esta sección también se pueden responder inspeccionando las capturas de pantalla de IGV a continuación.

Consulte la documentación de IGV para obtener más información.

IGV herramienta: Zoom a chr4: 540,000-560,000 (Cromosoma 4 entre 540 kb y 560 kb)

Pregunta pregunta

  1. ¿Qué información aparece en la parte superior como picos grises?
  2. ¿Qué indican las líneas de conexión entre algunas de las lecturas alineadas?

Solución solución

  1. La gráfica de cobertura: la suma de las lecturas mapeadas en cada posición.
  2. Indican eventos de unión (o sitios de empalme), es decir. lecturas que se asignan a través de un intrón

IGV herramienta: Inspeccione las uniones de empalme con una Parcela de sashimi

Comentar Creación de una trama de Sashimi

Pregunta pregunta

  1. ¿Qué representa el gráfico de barras verticales rojas? ¿Y los arcos con números?
  2. ¿Qué significan los números en los arcos?
  3. ¿Por qué observamos diferentes grupos apilados de cajas vinculadas en azul en la parte inferior?

Solución solución

  1. La cobertura de cada pista de alineación se traza como un gráfico de barras rojo. Los arcos representan uniones de empalme observadas, es decir., lee intrones que abarcan
  2. Los números se refieren al número de lecturas de unión observadas.
  3. Los diferentes grupos de cuadros enlazados en la parte inferior representan las diferentes transcripciones de los genes en esta ubicación, que están presentes en el archivo GTF.

Comentario comentario

Después del mapeo, tenemos la información sobre dónde se encuentran las lecturas en el genoma de referencia. También sabemos qué tan bien se mapearon. El siguiente paso en el análisis de datos de RNA -Seq es la cuantificación del número de lecturas asignadas a características genómicas (genes, transcripciones, exones,…).

Comentario comentario

La cuantificación depende tanto del genoma de referencia (el archivo FASTA) como de sus anotaciones asociadas (el archivo GTF). Es extremadamente importante usar un archivo de anotación que corresponda a la misma versión del genoma de referencia que usó para el mapeo (por ejemplo, dm6 aquí), ya que las coordenadas cromosómicas de los genes suelen ser diferentes entre las diferentes versiones del genoma de referencia.

Para identificar exones que están regulados por la Pasilla gen, necesitamos identificar genes y exones que se expresan diferencialmente entre muestras con agotamiento del gen PS (tratadas) y muestras de control (no tratadas). Luego analizaremos la expresión génica diferencial y también el uso diferencial del exón.


Genómica computacional con R

Con el advenimiento de las tecnologías de secuenciación de segunda generación (también conocida como próxima generación o de alto rendimiento), la cantidad de genes que pueden perfilarse para niveles de expresión con un solo experimento ha aumentado en el orden de decenas de miles de genes. Por lo tanto, el cuello de botella en este proceso se ha convertido en el análisis de datos más que en la generación de datos. Se requieren muchos métodos estadísticos y herramientas computacionales para obtener resultados significativos de los datos, que vienen con mucha información valiosa junto con muchas fuentes de ruido. Afortunadamente, la mayoría de los pasos del análisis de RNA-seq se han vuelto bastante maduros a lo largo de los años. A continuación, describiremos primero cómo llegar a una tabla de recuento de lecturas a partir de lecturas de fastq sin procesar obtenidas de una ejecución de secuenciación de Illumina. Luego demostraremos en R cómo procesar la tabla de conteo, hacer un análisis de expresión diferencial de casos y controles y hacer un análisis de enriquecimiento funcional posterior.

8.3.1 Procesamiento de datos brutos

8.3.1.1 Control de calidad y procesamiento de lectura

El primer paso en cualquier experimento que implique secuenciación de lectura corta de alto rendimiento debe ser verificar la calidad de secuenciación de las lecturas antes de comenzar a hacer cualquier análisis posterior. La calidad de las secuencias de entrada tiene una importancia fundamental en la confianza de las conclusiones biológicas extraídas del experimento. Hemos introducido el control de calidad y el procesamiento en el Capítulo 7, y esas herramientas y flujos de trabajo también se aplican en el análisis de RNA-seq.

8.3.1.2 Mejora de la calidad

El segundo paso en el flujo de trabajo de análisis de RNA-seq es mejorar la calidad de las lecturas de entrada. Este paso podría considerarse un paso opcional cuando la calidad de secuenciación es muy buena. Sin embargo, incluso con los conjuntos de datos de secuenciación de la más alta calidad, este paso aún puede mejorar la calidad de las secuencias de entrada. Los artefactos técnicos más comunes que se pueden filtrar son las secuencias adaptadoras que contaminan las lecturas secuenciadas y las bases de baja calidad que generalmente se encuentran al final de las secuencias. Las herramientas de uso común en el campo (trimmomatic (Bolger, Lohse y Usadel 2014), trimGalore (Andrews 2010)) nuevamente no están escritas en R, sin embargo, existen bibliotecas R alternativas para llevar a cabo la misma funcionalidad, por ejemplo, QuasR (Gaidatzis , Lerch, Hahne, et al.2015) (ver función QuasR :: preprocessReads) y ShortRead (Morgan, Anders, Lawrence, et al. 2009) (ver función ShortRead :: filterFastq). Algunos de estos enfoques se presentan en el Capítulo 7.

Los pasos de control de calidad de secuenciación y preprocesamiento de lectura se pueden visitar varias veces hasta lograr un nivel satisfactorio de calidad en los datos de secuencia antes de pasar a los pasos de análisis posteriores.

8.3.2 Alineación

Una vez que se alcanza un nivel decente de calidad en las secuencias, el nivel de expresión de los genes se puede cuantificar mapeando primero las secuencias a un genoma de referencia y, en segundo lugar, haciendo coincidir las lecturas alineadas con las anotaciones del gen, para contar el número de lecturas. mapeo de cada gen. Si la especie en estudio tiene un transcriptoma bien anotado, las lecturas se pueden alinear con las secuencias de la transcripción en lugar del genoma de referencia. En los casos en los que no existe un genoma o transcriptoma de referencia de buena calidad, es posible ensamblar de novo el transcriptoma a partir de las secuencias y luego cuantificar los niveles de expresión de genes / transcripciones.

Para las alineaciones de lectura de RNA-seq, además de la disponibilidad de genomas de referencia y anotaciones, probablemente el factor más importante a considerar al elegir una herramienta de alineación es si el método de alineación considera la ausencia de regiones intrónicas en las lecturas secuenciadas, mientras que el genoma objetivo puede contienen intrones. Por lo tanto, es importante elegir herramientas de alineación que tengan en cuenta empalmes alternativos. En el entorno básico donde una lectura, que se origina a partir de una secuencia de ADNc correspondiente a una unión exón-exón, debe dividirse en dos partes cuando se alinea contra el genoma. Existen varias herramientas que consideran este factor como STAR (Dobin, Davis, Schlesinger, et al. 2013), Tophat2 (Kim, Pertea, Trapnell, et al. 2013), Hisat2 (Kim, Langmead y Salzberg 2015), y GSNAP (Wu, Reeder, Lawrence, et al. 2016). La mayoría de las herramientas de alineación están escritas en lenguajes C / C ++ debido a problemas de rendimiento. También hay bibliotecas de R que pueden realizar alineaciones de lectura breves, que se describen en el Capítulo 7.

8.3.3 Cuantificación

Una vez que las lecturas se alinean con el objetivo, se debería haber obtenido un archivo SAM / BAM ordenado por coordenadas. El archivo BAM contiene toda la información relacionada con la alineación de todas las lecturas que se han intentado alinear con la secuencia objetivo. Esta información consiste, básicamente, en las coordenadas genómicas (cromosoma, inicio, final, hebra) de dónde se emparejó una secuencia (si es que lo hizo) en el objetivo, inserciones / deleciones / desajustes específicos que describen las diferencias entre la entrada y el objetivo secuencias. Estas piezas de información se utilizan junto con las coordenadas genómicas de las anotaciones del genoma, como los modelos de genes / transcripciones, para contar cuántas lecturas se han secuenciado a partir de un gen / transcripción. Por simple que parezca, no es una tarea trivial asignar lecturas a un gen / transcripción simplemente comparando las coordenadas genómicas de las anotaciones y las secuencias, debido a factores de confusión como anotaciones de genes superpuestas, anotaciones de exones superpuestas de diferentes transcripciones. isoformas de un gen y anotaciones superpuestas de cadenas de ADN opuestas en ausencia de un protocolo de secuenciación específico de cadena. Por lo tanto, para el recuento de lecturas, es importante tener en cuenta:

  1. Especificidad de hebra del protocolo de secuenciación: ¿Se espera que las lecturas se originen en la hebra directa, la hebra inversa o inespecíficas?
  2. Modo de conteo: - cuando se cuenta a nivel de gen: cuando hay anotaciones superpuestas, ¿a qué características se debe asignar la lectura? Las herramientas suelen tener un parámetro que permite al usuario seleccionar un modo de conteo. - al contar a nivel de transcripción: cuando hay múltiples isoformas de un gen, ¿a qué isoforma debe asignarse la lectura? Esta consideración suele ser una consideración algorítmica que el usuario final no puede modificar.

Algunas herramientas pueden acoplar la alineación a la cuantificación (por ejemplo, STAR), mientras que algunas asumen que las alineaciones ya están calculadas y requieren archivos BAM como entrada. Por otro lado, en presencia de buenas anotaciones de transcriptoma, los métodos sin alineación (Salmon (Patro, Duggal, Love, et al. 2017), Kallisto (Bray, Pimentel, Melsted, et al. 2016), Sailfish (Patro, Mount y Kingsford 2014)) también se pueden utilizar para estimar los niveles de expresión de transcripciones / genes. También existen métodos de cuantificación sin referencia que primero pueden ensamblar de novo el transcriptoma y estimar los niveles de expresión basados ​​en este ensamblaje. Esta estrategia puede ser útil para descubrir transcripciones novedosas o puede ser necesaria en los casos en que no exista una buena referencia. Si existe un transcriptoma de referencia pero de baja calidad, se puede utilizar un ensamblador de transcriptoma basado en referencias como Cufflinks (Trapnell, Williams, Pertea, et al. 2010) para mejorar el transcriptoma. En caso de que no haya una anotación de transcriptoma disponible, se puede utilizar un ensamblador de novo como Trinity (Haas, Papanicolaou, Yassour, et al. 2013) o Trans-ABySS (Robertson, Schein, Chiu, et al. 2010) para ensamblar el transcriptoma desde cero.

Dentro de R, la cuantificación se puede realizar usando: - Rsubread :: featureCounts - QuasR :: qCount - GenomicAlignments :: summaryizeOverlaps

8.3.4 Normalización dentro de la muestra de los recuentos de lectura

La aplicación más común después de cuantificar la expresión de un gen (como el número de lecturas alineadas con el gen), es comparar la expresión del gen en diferentes condiciones, por ejemplo, en un entorno de casos y controles (p. Ej., Enfermedad versus normal) o en un series de tiempo (por ejemplo, a lo largo de diferentes etapas de desarrollo). Hacer tales comparaciones ayuda a identificar los genes que podrían ser responsables de una enfermedad o una trayectoria de desarrollo alterada. Sin embargo, hay varias advertencias que deben abordarse antes de hacer una comparación entre los recuentos de lectura de un gen en diferentes condiciones (Maza, Frasse, Senin, et al. 2013).

  • El tamaño de la biblioteca (es decir, la profundidad de secuenciación) varía entre muestras procedentes de diferentes carriles de la celda de flujo de la máquina de secuenciación.
  • Los genes más largos tendrán un mayor número de lecturas.
  • La composición de la biblioteca (es decir, el tamaño relativo del transcriptoma estudiado) puede ser diferente en dos condiciones biológicas diferentes.
  • Los sesgos de contenido de GC en diferentes muestras pueden conducir a un muestreo sesgado de genes (Risso, Schwartz, Sherlock, et al. 2011).
  • La cobertura de lectura de una transcripción puede estar sesgada y distribuida de manera no uniforme a lo largo de la transcripción (Mortazavi, Williams, McCue, et al. 2008).

Por lo tanto, estos factores deben tenerse en cuenta antes de realizar comparaciones.

Los enfoques de normalización más básicos abordan el sesgo de profundidad de secuenciación. Dichos procedimientos normalizan el recuento de lecturas por gen dividiendo el recuento de lecturas de cada gen por un valor determinado y multiplicándolo por 10 ^ 6. Estos valores normalizados suelen denominarse CPM (recuentos por millón de lecturas):

  • Normalización de recuentos totales (divida los recuentos por suma de todos los conteos)
  • Normalización del cuartil superior (dividir los recuentos por Cuartilla superior valor de los recuentos)
  • Normalización mediana (dividir los recuentos por mediana de todos los conteos)

Las métricas populares que mejoran el CPM son RPKM / FPKM (lecturas / fragmentos por kilobase de millón de lecturas) y TPM (transcripciones por millón). RPKM se obtiene dividiendo el valor de CPM por otro factor, que es la longitud del gen por kilobase. FPKM es lo mismo que RPKM, pero se utiliza para lecturas de extremo emparejado. Por tanto, los métodos RPKM / FPKM dan cuenta, en primer lugar, de la tamaño de la biblioteca, y en segundo lugar, el longitudes de genes.

TPM also controls for both the library size and the gene lengths, however, with the TPM method, the read counts are first normalized by the gene length (per kilobase), and then gene-length normalized values are divided by the sum of the gene-length normalized values and multiplied by 10^6. Thus, the sum of normalized values for TPM will always be equal to 10^6 for each library, while the sum of RPKM/FPKM values do not sum to 10^6. Therefore, it is easier to interpret TPM values than RPKM/FPKM values.

8.3.5 Computing different normalization schemes in R

Here we will assume that there is an RNA-seq count table comprising raw counts, meaning the number of reads counted for each gene has not been exposed to any kind of normalization and consists of integers. The rows of the count table correspond to the genes and the columns represent different samples. Here we will use a subset of the RNA-seq count table from a colorectal cancer study. We have filtered the original count table for only protein-coding genes (to improve the speed of calculation) and also selected only five metastasized colorectal cancer samples along with five normal colon samples. There is an additional column width that contains the length of the corresponding gene in the unit of base pairs. The length of the genes are important to compute RPKM and TPM values. The original count tables can be found from the recount2 database (https://jhubiostatistics.shinyapps.io/recount/) using the SRA project code SRP029880, and the experimental setup along with other accessory information can be found from the NCBI Trace archive using the SRA project code SRP029880`.

8.3.5.1 Computing CPM

Let’s do a summary of the counts table. Due to space limitations, the summary for only the first three columns is displayed.

To compute the CPM values for each sample (excluding the width column):

Check that the sum of each column after normalization equals to 10^6 (except the width column).


Fondo

The remarkable diversity of flower colors, especially in wild plants has fascinated botanists, ecologists, and horticulturists for centuries [1,2,3]. The coloring of floral organs, a remarkable character of flowering plants, is a striking feature of the angiosperm radiation [4, 5]. Flower color diversity is recognized to be one of key adaptive traits correlated predominantly with pollinators (e.g. insects, birds) and animals for seed dispersal [6, 7]. Moreover, the flower color phenotype is an important feature for plants used for their classification by taxonomists. However, flower color appears evolutionarily to be one of the most labile traits, down to populations in the same species [7, 8].

The cellular compounds of flowers that contribute to the color profile and visually perceived by humans are generally referred to as “pigments”. A group of secondary metabolites belonging to flavonoids are the main determinants of pigments for coloration in plants, where anthocyanins are responsible for red orange to red, purple to violet pigments found in flowers, leaves, fruits, seeds and other tissues [9, 10]. Anthocyanins are the predominant compounds of floral coloration, existing in over 90% of angiosperms [11]. The flavonoid biosynthetic pathway leading to accumulation of anthocyanins is highly conserved and well characterized, and has been extensively studied in many species, most of which are in model plants or agriculturally and horticulturally important plants [12,13,14,15]. Few studies have examined the molecular basis underlying the formation and accumulation of anthocyanin in wild species [16, 17]. Based on these studies, three major associated factors have been proposed to be involved in anthocyanin accumulation, including transcription regulatory genes (MYB-bHLH-WD40 complex) that occur in the nucleus, structural genes (CHS, FLS, DFR, ANS) acting in the biosynthetic pathway, and transporter genes (GST) transferring anthocyanin from the cytosol into the vacuole [10, 18, 19]. The expression of these genes could also be affected by natural variation in sequences and cis-regulatory elements as well as epigenetic modifications (such as DNA methylation) in the promoter regions [18, 20]. Moreover, the color of flowers can be stabilized and enhanced by co-pigmentation of anthocyanins by flavonols, where it is observed as hyperchromic effect, in which the intensity of an anthocyanin content is fortified [21]. For instance, the DFR gene along with the FLS gene can compete for a substrate leading to the production of different anthocyanins and flavanols through two primary branches [22, 23], thus resulting in co-pigmentation. In contrast to the biosynthesis pathways, knowledge of anthocyanin catabolism in plants is limited. Some catabolic genes like BGLU y PER have been shown to be responsible for anthocyanin degradation [24, 25]. Nevertheless, the molecular mechanism regulating anthocyanin synthesis has been shown to vary among plant species resulting in structural diversity of anthocyanins, because the biosynthesis pathway is regulated by multiple factors through regulatory networks [26].

Color is a form of electromagnetic radiation in the range of the visible spectrum. The wavelengths reflected by pigments determine the color of a flower [27]. Color can be defined and classified in terms of Brightness (the intensity of a signal, B), Saturation (the purity of a color, S) and Hue (the spectral descriptor of color, H), and those features are commonly used to distinguish colors [27, 28]. Brightness refers to the color intensity that is determined by the amount of anthocyanin [29, 30], and different color component combinations such as B/H, S/H were found to be significantly correlated with anthocyanin content as well [31]. Liu et al. [32] proposed that the color brightness decreased as the total anthocyanin content increased. It was also demonstrated that a correlation exists between the saturation/hue ratio (S/H) and anthocyanin content [31]. With these parameters, the anthocyanin content can be rapidly and non-destructively determined.

In evergreen azaleas (Rhododendron), anthocyanins and flavanols are the main flower pigments, and especially the composition of anthocyanin constituents (i.e. cyanidin, delphinidin, malvidin, pelargonidin, peonidin, and petunidin), and their quantities determine their flower color that ranges from light pink to violet [11, 33]. Some studies have reported that R. kiusianum with purple-colored flowers contain derivatives of both anthocyanidins cyanidin and delphinidin, whereas the red-colored flowers of R. kaempferi contain only cyanidin derivatives [34]. Le Maitre et al. [35, 36] studying Erica species, belonging to the same family Ericaceae as Rhododendron, used qRT-PCR and UPLC-MS, unraveled the anthocyanin genetic network of floral color shifts between red or pink and white or yellow flowered species and found losses of single pathway gene expression, abrogation of the entire pathway due to loss of the expression of a transcription factor or loss of function mutations in pathway genes resulted in striking floral color shifts.

Here, we investigated the genetic basis of flower coloration using a highly color polymorphic Rhododendron sanguineum complex. The complex (R. subgen. Hymenanthes) includes plants with yellow to pink or crimson to blackish crimson flowers that are classified into six varieties mainly based on their flower color differences [37]. Members of this complex are basically located at high elevations (> 3000 m) associated with snow cover [37]. They are endemic to northwest Yunnan and southeast Tibet, one of the global biodiversity hotspots [38]. This region is also recognized as one of the centers for diversification and differentiation of Rhododendron [37, 39]. The flower color polymorphisms of this genus have been traditionally viewed as an ecologically adaptive trait that is essential in attracting specific pollinators [40,41,42], and may also be the response to environmental variation, such as UV radiations at different elevations, temperatures, and soil conditions [32]. Although there are studies published on the anthocyanin components and contents in Rhododendron flowers, most were solely dedicated to the identification of the pigment constituents in the petals of some wild and cultivated azaleas using thin-layer chromatography (TLC) and high-performance liquid chromatography (HPLC) [11, 33]. No study so far focused on the molecular mechanisms underlying infraspecific color polymorphisms in Rhododendron. The study of closely related entities such as a species complex has the advantage of a fairly homogeneous genetic background where flower color genes vary and cases of homoplasy are limited. Previous studies mainly focused on color shifts at different developmental stages of a single species [14, 18], or covered a number of related species [26, 35].

In the present study, we combined transcriptome sequencing (RNA-seq) and genome resequencing with reflectance spectra analyses to elucidate molecular and anthocyanin content differences among three differently colored naturally occurring varieties of the R. sanguineum complex, with yellow flushed pink to deep blackish crimson colored flowers. We aimed at studying the correlation between infraspecific flower color variation and the expression of candidate genes of the anthocyanin / flavonoid biosynthesis pathway. Our findings may allow the proposal of a hypothesis for the genetic mechanism of the expression of flower color variation and a representative case of pollinator-mediated incipient sympatric speciation in the R. sanguineum complex. In addition, it is the first study to compare transcriptome profiles in a natural system of a non-model species of Rhododendron.


3 RESULTS

3.1 Especially progression genes show elevated expression at the onset of tooth patterning

For a robust readout of gene expression profiles, we first obtained gene expression levels using both microarray and RNAseq techniques from E13 (bud stage) and E14 (cap stage) mouse molars (Section 2). From dissected tooth germs we obtained five microarray and seven RNAseq replicates for both developmental stages. The results show that especially the progression category genes (genes required for the progression of tooth development) are highly expressed during E13 compared to the control gene sets (tissue, dispensable, and developmental-process categories, pag values range from .0003 to .0426 for RNAseq and microarray experiments, tested using random resampling, for details and all the tests, see Section 2, Figure 2, and Tables S2 and S3). Comparable differences are observed in E14 molars (pag values range from .0000 to .0466 Figure 2, Tables S2 and S3).

In general, the expression differences between progression and tissue categories appear greater than between progression and dispensable categories (pag values range from .0028 to .0379 and .0059 to .0466, respectively Table S3), suggesting that some of the genes in the dispensable category may still play a functional role in tooth development. In our data we have 11 genes that cause a developmental arrest of the tooth when double mutated (Appendix S1). The expression level of this double-mutant category shows incipient upregulation compared to that of the developmental-process category (pag values range from .0322 to .1637 Table S3), but not when compared to the tissue or dispensable categories (pag values range from .0978 to .5010 Table S3). Therefore, it is plausible, based on the comparable expression levels between double and some of the dispensable category genes, that several of the genes in the dispensable category may cause phenotypic effects when mutated in pairs.

Even though expression levels of the shape category genes (genes required for normal shape development) are lower than that of the progression category (Figure 2), at least the E14 microarray data suggests elevated expression levels relative to all the other control categories (pag values range from .0001 to .0901 Table S3). The moderately elevated levels of expression by the shape category genes could indicate that they are required slightly later in development, or that the most robust upregulation happens only for genes that are essential for the progression of the development. The latter option seems to be supported by a RNAseq analysis of E16 molar, showing only slight upregulation of shape category genes in the bell stage molars (Table S3).

3.2 Transcriptomes of developing rat molars show elevated expression of the progression genes

Because our gene categories were based on experimental evidence from the mouse, we also tested whether comparable expression levels can be detected for the same genes in the rat. Evolutionary divergence of Mus-Rattus dates back to the Middle Miocene (Kimura et al., 2015 ), allowing a modest approximation of conservation in the expression levels. Examination of bud (E15) and cap (E17) stage RNAseq of rat molars shows comparable upregulation of progression and shape category genes as in the mouse (Figure 3, Table S2 and S3). Considering also that many of the null mutations in keystone gene in the mouse are known to have comparable phenotypic effects in humans (Nieminen, 2009 ), our keystone gene categories and analyses are likely to apply to mammalian teeth in general.

One complication of our expression level analyses is that these have been done at the whole organ level. Because many of the genes regulating tooth development are known to have spatially complex expression patterns within the tooth (Nieminen at al. 1998 ), cell-level examinations are required to decompose the patterns within the tissue.

3.3 Single-cell RNAseq reveals cell-level patterns of keystone genes

Tooth development is punctuated by iteratively forming epithelial signaling centers, the enamel knots. The first, primary enamel knot, is active in E14 mouse molar and at this stage many genes are known to have complex expression patterns. Some progression category genes have been reported to be expressed in the enamel knot, whereas others have mesenchymal or complex combinatorial expression patterns (Jernvall & Thesleff, 2012 Nieminen et al., 1998 ). To quantify these expression levels at the cell-level, we performed a single-cell RNAseq (scRNAseq) on E14 mouse molars (Section 2). We focused on capturing a representative sample of cells by dissociating each tooth germ without cell sorting (norte = 4). After data filtering, 7000–8811 cells per tooth were retained for the analyses, providing 30,930 aggregated cells for a relatively good proxy of the E14 mouse molar (Section 2).

First we examined whether the scRNAseq produces comparable expression levels to our previous analyses. For the comparisons, the gene count values from the cells were summed up and treated as bulk RNAseq data (Figure 4a and Section 2). We analyzed the expression levels of different gene categories as in the mouse bulk data (Figure 2) and the results show a general agreement between the experiments (Figures 2 and 4b). As in the previous analyses (Table S3), the progression category shows the highest expression levels compared to the control gene sets (pag values range from .0071 to .0310 Table S3). Although the mean expression of the shape category is intermediate between progression and control gene sets, scRNAseq shape category is not significantly upregulated in the randomization tests (pag values range from .7788 to .9968). This pattern reflects the bulk RNAseq analyses (for both mouse and rat) while the microarray analysis showed slightly stronger upregulation, suggesting subtle differences between the methodologies (the used mouse strain was also different in the microarray experiment).

Unlike the bulk transcriptome data, the scRNAseq data can be used to quantify the effect of expression domain size on the overall expression level of a gene. The importance of expression domain size is well evident in the scRNAseq data when we calculated the number of cells that express each gene (Section 2). The data shows that the overall tissue level gene expression is highly correlated with the cell population size (Figure 5a). In other words, the size of the expression domain is the key driver of expression levels measured at the whole tissue level.

To examine the cell level patterns further, we calculated the mean transcript abundances for each gene for the cells that express that gene (see Section 2). This metric approximates the cell-level upregulation of a particular gene, and is thus independent of the size of the expression domain. We calculated the transcript abundance values for progression, shape, tissue, double, and dispensable category genes in each cell that expresses any of those genes. The resulting mean transcript abundances were contrasted to that of the dispensable category (Section 2). The results show that the average transcript abundance is high in the progression category whereas the other categories show roughly comparable transcript abundances (Figure 5b). Considering that the progression category genes have highly heterogeneous expression patterns (e.g., Nieminen at al. 1998 Figure 5c), their high cell-level transcript abundance (Figure 5b) is suggestive about their critical role at the cell level. That is, progression category genes are not only highly expressed at the tissue level because they have broad expression domains, but rather because they are upregulated in individual cells irrespective of domain identity or size. These results suggest that high cell-level transcript abundance is a system-level feature of genes essential for the progression of tooth development, a pattern that seems to be shared with essential genes of single cell organisms (Dong et al., 2018 ). We note that although the dispensable category has several genes showing comparable expression levels with that of the progression category genes at the tissue level (Figure 2), their cell-level transcript abundances are predominantly low (Figure 5b).

Next we examined more closely the differences between progression and shape category genes, and to what extent the upregulation of the keystone genes reflects the overall expression of the corresponding pathways.

3.4 Keystone gene upregulation in the context of their pathways

In our data the developmental-process genes appear to have slightly elevated expression levels compared to the other protein coding genes (Figures 2, 3, and 4b), suggesting an expected and general recruitment of the pathways required for organogenesis. To place the progression and shape category genes into the specific context of their corresponding pathways, we investigated in E14 mouse bulk RNAseq whether the pathways implicated in tooth development show elevated expression levels. Six pathways, Fgf, Wnt, Tgfβ, Hedgehog (Hh), Notch, and Ectodysplasin (Eda), contain the majority of progression and shape genes (Section 2). First we used the RNAseq of E14 stage molars to test whether these pathways show elevated expression levels. We manually identified 272 genes belonging to the six pathways (Section 2 and Table S4). Comparison of the median expression levels of the six-pathway genes with the developmental-process genes shows that the pathway genes are a highly upregulated set of genes (Figure 6a pag < .0001, random resampling). This difference suggests that the experimentally identified progression and shape genes might be highly expressed partly because they belong to the developmentally upregulated pathways. To specifically test this possibility, we contrasted the expression levels of the progression and shape genes to the genes of their corresponding signaling families.

The 15 progression category genes belong to four signaling families (Wnt, Tgfβ, Fgf, Hh) with 221 genes in our tabulations. Even though these pathways are generally upregulated in the E14 tooth, the median expression level of the progression category is still further elevated (Figure 6b pag < .0001). In contrast, the analyses for the 28 shape category genes and their corresponding pathways (272 genes from Wnt, Tgfβ, Fgf, Hh, Eda, Notch) show comparable expression levels (Figure 6c pag = .5919). Whereas this contrasting pattern between progression and shape genes within their pathways may explain the subtle upregulation of the shape category (Figure 2), the difference warrants a closer look. Examination of the two gene categories reveals that compared to the progression category genes, relatively large proportion of the shape category genes are ligands (36% shape genes compared to 20% progression genes, Appendix S1). In our E14 scRNAseq data, ligands show generally smaller expression domains than other genes (roughly by half, Figure 6d,e), and the low expression of the shape category genes seems to be at least in part driven by the ligands (Figure 6c and Table S5).

Overall, the upregulation of the keystone genes within their pathways appears to be influenced by the kind of proteins they encode. In this context it is noteworthy that patterning of tooth shape requires spatial regulation of secondary enamel knots and cusps, providing a plausible explanation for the high proportion of genes encoding diffusing ligands in the shape category.


  1. Trapnell, C., L. Pachter, and S. L. Salzberg, 2009 TopHat: discovering splice junctions with RNA-Seq. Bioinformatics 25: 1105–1111. https://academic.oup.com/bioinformatics/article/25/9/1105/203994
  2. Levin, J. Z., M. Yassour, X. Adiconis, C. Nusbaum, D. A. Thompson et al., 2010 Comprehensive comparative analysis of strand-specific RNA sequencing methods. Nature Methods 7: 709. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3005310/
  3. Young, M. D., M. J. Wakefield, G. K. Smyth, and A. Oshlack, 2010 Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biology 11: R14. https://genomebiology.biomedcentral.com/articles/10.1186/gb-2010-11-2-r14
  4. Brooks, A. N., L. Yang, M. O. Duff, K. D. Hansen, J. W. Park et al., 2011 Conservation of an RNA regulatory map between Drosophila and mammals. Genome Research 21: 193–202. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3032923/
  5. Marcel, M., 2011 Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.journal 17: 10–12. http://journal.embnet.org/index.php/embnetjournal/article/view/200
  6. Robinson, J. T., H. Thorvaldsdóttir, W. Winckler, M. Guttman, E. S. Lander et al., 2011 Integrative genomics viewer. Nature Biotechnology 29: 24. https://www.nature.com/nbt/journal/v29/n1/abs/nbt.1754.html
  7. Wang, L., S. Wang, and W. Li, 2012 RSeQC: quality control of RNA-seq experiments. Bioinformatics 28: 2184–2185. https://www.ncbi.nlm.nih.gov/pubmed/22743226
  8. Dobin, A., C. A. Davis, F. Schlesinger, J. Drenkow, C. Zaleski et al., 2013 STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29: 15–21. https://academic.oup.com/bioinformatics/article/29/1/15/272537
  9. Kim, D., G. Pertea, C. Trapnell, H. Pimentel, R. Kelley et al., 2013 TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biology 14: R36. https://genomebiology.biomedcentral.com/articles/10.1186/gb-2013-14-4-r36
  10. Liao, Y., G. K. Smyth, and W. Shi, 2013 featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 30: 923–930. https://academic.oup.com/bioinformatics/article/31/2/166/2366196
  11. Luo, W., and C. Brouwer, 2013 Pathview: an R/Bioconductor package for pathway-based data integration and visualization. Bioinformatics 29: 1830–1831. https://academic.oup.com/bioinformatics/article-abstract/29/14/1830/232698
  12. Love, M. I., W. Huber, and S. Anders, 2014 Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology 15: 550. https://genomebiology.biomedcentral.com/articles/10.1186/s13059-014-0550-8
  13. Anders, S., P. T. Pyl, and W. Huber, 2015 HTSeq—a Python framework to work with high-throughput sequencing data. Bioinformatics 31: 166–169. https://academic.oup.com/bioinformatics/article/31/2/166/2366196
  14. Kim, D., B. Langmead, and S. L. Salzberg, 2015 HISAT: a fast spliced aligner with low memory requirements. Nature Methods 12: 357. https://www.nature.com/articles/nmeth.3317
  15. Ewels, P., M. Magnusson, S. Lundin, and M. Käller, 2016 MultiQC: summarize analysis results for multiple tools and samples in a single report. Bioinformatics 32: 3047–3048. https://academic.oup.com/bioinformatics/article/32/19/3047/2196507
  16. Thurmond, J., J. L. Goodman, V. B. Strelets, H. Attrill, L. S. Gramates et al., 2018 FlyBase 2.0: the next generation. Nucleic Acids Research 47: D759–D765. https://academic.oup.com/nar/article-abstract/47/D1/D759/5144957
  17. Kim, D., J. M. Paggi, C. Park, C. Bennett, and S. L. Salzberg, 2019 Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nature Biotechnology 37: 907–915. https://www.nature.com/articles/s41587-019-0201-4

Did you use this material as an instructor? Feel free to give us feedback on how it went.