Identificación de empalme alternativo y poliadenilación en datos de RNA-seq

Summary

El empalme alternativo (AS) y la poliadenilación alternativa (APA) amplían la diversidad de isoformas de transcripción y sus productos. Aquí, describimos protocolos bioinformáticos para analizar RNA-seq a granel y ensayos de secuenciación final 3' para detectar y visualizar AS y APA que varían según las condiciones experimentales.

Abstract

Además del análisis típico de RNA-Seq para medir la expresión génica diferencial (DGE) en condiciones experimentales / biológicas, los datos de RNA-seq también se pueden utilizar para explorar otros mecanismos reguladores complejos a nivel de exón. El empalme alternativo y la poliadenilación juegan un papel crucial en la diversidad funcional de un gen al generar diferentes isoformas para regular la expresión génica a nivel post-transcripcional, y limitar los análisis a todo el nivel del gen puede pasar por alto esta importante capa reguladora. Aquí, demostramos análisis detallados paso a paso para la identificación y visualización del uso diferencial del exón y el sitio de poliadenilación en todas las condiciones, utilizando Bioconductor y otros paquetes y funciones, incluidos DEXSeq, diffSplice del paquete Limma y rMATS.

Introduction

RNA-seq ha sido ampliamente utilizado a lo largo de los años, generalmente para estimar la expresión génica diferencial y el descubrimiento de genes¹. Además, también se puede utilizar para estimar el uso variable del nivel de exón debido a que los genes expresan diferentes isoformas, lo que contribuye a una mejor comprensión de la regulación génica a nivel post-transcripcional. La mayoría de los genes eucariotas generan diferentes isoformas mediante empalme alternativo (AS) para aumentar la diversidad de la expresión del ARNm. Los eventos AS se pueden dividir en diferentes patrones: salto de exones completos (SE) donde un exón ("casete") se elimina completamente de la transcripción junto con sus intrones flanqueantes; selección alternativa (donante) del sitio de empalme de 5' (A5SS) y selección alternativa 3' (aceptor) del sitio de empalme (A3SS) cuando dos o más sitios de empalme están presentes en cada extremo de un exón; retención de intrones (RI) cuando un intrón se retiene dentro de la transcripción de ARNm maduro y exclusión mutua del uso de exones (MXE) donde solo uno de los dos exones disponibles puede ser retenido a la vez ^2,3. La poliadenilación alternativa (APA) también juega un papel importante en la regulación de la expresión génica utilizando sitios alternativos de poli (A) para generar múltiples isoformas de ARNm a partir de una sola transcripción⁴. La mayoría de los sitios de poliadenilación (pAs) se encuentran en la región 3' no traducida (3' UTRs), generando isoformas de ARNm con diversas longitudes 3' UTR. Como el 3' UTR es el eje central para reconocer elementos reguladores, diferentes longitudes 3' UTR pueden afectar la localización, estabilidad y traducción del ARNm⁵. Hay una clase de ensayos de secuenciación final 3' optimizados para detectar APA que difieren en los detalles del protocolo⁶. La canalización descrita aquí está diseñada para PolyA-seq, pero se puede adaptar para otros protocolos como se describe.

En este estudio, presentamos una tubería de métodos de análisis de exones diferenciales^7,8 (Figura 1), que se pueden dividir en dos grandes categorías: basados en exones (DEXSeq^9, diffSplice¹⁰) y basados en eventos (replicate Multivariate Analysis of Transcript Splicing (rMATS)¹¹). Los métodos basados en exones comparan el cambio de pliegue a través de las condiciones de los exones individuales, contra una medida del cambio general del pliegue genético para llamar al uso de exones expresado diferencialmente, y a partir de eso calcular una medida a nivel de gen de la actividad de AS. Los métodos basados en eventos utilizan lecturas de unión de expansión exón-intrón para detectar y clasificar eventos de empalme específicos, como la omisión de exón o la retención de intrones, y distinguir estos tipos de AS en la salida³. Por lo tanto, estos métodos proporcionan vistas complementarias para un análisis completo de la EA^12,13. Se seleccionaron DEXSeq (basado en el paquete DESeq2¹⁴ DGE) y diffSplice (basado en el paquete Limma¹⁰ DGE) para el estudio, ya que se encuentran entre los paquetes más utilizados para el análisis de empalme diferencial. rMATS fue elegido como un método popular para el análisis basado en eventos. Otro método popular basado en eventos es MISO (mezcla de isoformas)¹. Para APA adaptamos el enfoque basado en exones.

Figura 1. Pipeline de análisis. Diagrama de flujo de los pasos utilizados en el análisis. Los pasos incluyen: obtener los datos, realizar controles de calidad y alineación de lecturas seguidas de contar lecturas utilizando anotaciones para exones conocidos, intrones y sitios pA, filtrado para eliminar recuentos bajos y normalización. Los datos de PolyA-seq se analizaron para sitios de pA alternativos utilizando métodos diffSplice/DEXSeq, RNA-Seq a granel se analizaron para splicing alternativo a nivel de exón con métodos diffSplice/DEXseq, y los eventos de AS se analizaron con rMATS. Haga clic aquí para ver una versión más grande de esta figura.

Los datos de RNA-seq utilizados en este estudio fueron adquiridos de Gene Expression Omnibus (GEO) (GSE138691)¹⁵. Utilizamos datos de ARN-seq de ratón de este estudio con dos grupos de condiciones: tipo salvaje (WT) y knockout tipo 1 similar a Muscleblind (Mbnl1 KO) con tres réplicas cada uno. Para demostrar el análisis diferencial de uso del sitio de poliadenilación, obtuvimos datos de PolyA-seq de fibroblastos embrionarios de ratón (MEF) (GEO Accesion GSE60487)¹⁶. Los datos tienen cuatro grupos de condiciones: tipo salvaje (WT), tipo ciego muscular tipo 1 / tipo 2 doble knockout (Mbnl1/2 DKO), Mbnl 1/2 DKO con derribo Mbnl3 (KD) y Mbnl1/2 DKO con control Mbnl3 (Ctrl). Cada grupo de condición consta de dos réplicas.

	Adhesión al GEO	Número de ejecución de SRA	Nombre de la muestra	Condición	Replicar	Tejido	Secuenciación	Longitud de lectura
RNA-Seq	GSM4116218	SRR10261601	Mbnl1KO_Thymus_1	Mbnl1 nocaut	Rep 1	Timo	Extremo emparejado	100 pb
	GSM4116219	SRR10261602	Mbnl1KO_Thymus_2	Mbnl1 nocaut	Rep 2	Timo	Extremo emparejado	100 pb
	GSM4116220	SRR10261603	Mbnl1KO_Thymus_3	Mbnl1 nocaut	Rep 3	Timo	Extremo emparejado	100 pb
	GSM4116221	SRR10261604	WT_Thymus_1	Tipo salvaje	Rep 1	Timo	Extremo emparejado	100 pb
	GSM4116222	SRR10261605	WT_Thymus_2	Tipo salvaje	Rep 2	Timo	Extremo emparejado	100 pb
	GSM4116223	SRR10261606	WT_Thymus_3	Tipo salvaje	Rep 3	Timo	Extremo emparejado	100 pb
3P-seq	GSM1480973	SRR1553129	WT_1	Tipo salvaje (WT)	Rep 1	Fibroblastos embrionarios de ratón (MEF)	Extremo único	40 pb
	GSM1480974	SRR1553130	WT_2	Tipo salvaje (WT)	Rep 2	Fibroblastos embrionarios de ratón (MEF)	Extremo único	40 pb
	GSM1480975	SRR1553131	DKO_1	Mbnl 1/2 doble knockout (DKO)	Rep 1	Fibroblastos embrionarios de ratón (MEF)	Extremo único	40 pb
	GSM1480976	SRR1553132	DKO_2	Mbnl 1/2 doble knockout (DKO)	Rep 2	Fibroblastos embrionarios de ratón (MEF)	Extremo único	40 pb
	GSM1480977	SRR1553133	DKOsiRNA_1	Mbnl 1/2 doble knockout con Mbnl 3 siRNA (KD)	Rep 1	Fibroblastos embrionarios de ratón (MEF)	Extremo único	40 pb
	GSM1480978	SRR1553134	DKOsiRNA_2	Mbnl 1/2 doble knockout con Mbnl 3 siRNA (KD)	Rep 2	Fibroblastos embrionarios de ratón (MEF)	Extremo único	36 pb
	GSM1480979	SRR1553135	DKONTsiRNA_1	Mbnl 1/2 doble knockout con siRNA no dirigido (Ctrl)	Rep 1	Fibroblastos embrionarios de ratón (MEF)	Extremo único	40 pb
	GSM1480980	SRR1553136	DKONTsiRNA_2	Mbnl 1/2 doble knockout con siRNA no dirigido (Ctrl)	Rep 2	Fibroblastos embrionarios de ratón (MEF)	Extremo único	40 pb

Tabla 1. Resumen de los conjuntos de datos RNA-Seq y PolyA-seq utilizados para el análisis.

Protocol

1. Instalación de herramientas y paquetes R utilizados en el análisis

1. Conda es un administrador de paquetes popular y flexible que permite la instalación conveniente de paquetes con sus dependencias en todas las plataformas. Use 'Anaconda' (administrador de paquetes conda) para instalar 'conda' que se puede usar para instalar las herramientas / paquetes necesarios para el análisis.
2. Descargue 'Anaconda' de acuerdo con los requisitos del sistema de https://www.anaconda.com/products/individual#Downloads e instálelo siguiendo las instrucciones en el instalador gráfico. Instale todos los paquetes necesarios usando 'conda' escribiendo lo siguiente en la línea de comandos de Linux.
   conda install -c daler sratoolkit
conda install -c conda-forge parallel
conda install -c bioconda star bowtie fastqc rmats rmats2sashimiplot samtools fasterq-dump cutadapt bedtools deeptools
3. Para descargar todos los paquetes de R utilizados en el protocolo, escriba el código siguiente en la consola de R (iniciada en la línea de comandos de Linux escribiendo 'R') o en la consola de Rstudio.
   bioc_packages<- c("DEXSeq", "Rsubread", "EnhancedVolcano", "edgeR", "limma", "maser","GenomicRanges")
packages<- c("magrittr", "rtracklayer", "tidyverse", "openxlsx", "BiocManager")
#Install if not already installed
installed_packages<-packages%in% rownames(installed.packages())
installed_bioc_packages<-bioc_packages%in% rownames(installed.packages())
if(any(installed_packages==FALSE)) {
install.packages(packages[!installed_packages],dependencies=TRUE)
BiocManager::install(packages[!installed_bioc_packages], dependencies=TRUE)
}
   NOTA: En este protocolo computacional, los comandos se darán como archivos de R Notebook (archivos con extensión ". Rmd"), archivos de código R (archivos con extensión ". R"), o scripts de shell Bash de Linux (archivos con extensión ".sh"). Los archivos de R Notebook (Rmd) deben abrirse en RStudio usando File| Abra Archivo..., y los fragmentos de código individuales (que pueden ser comandos de R o comandos de shell Bash) se ejecutan interactivamente haciendo clic en la flecha verde en la parte superior derecha. Los archivos de código R se pueden ejecutar abriendo en RStudio, o en la línea de comandos de Linux anteponiendo "Rscript", por ejemplo, Rscript example.R. Los scripts de shell se ejecutan en la línea de comandos de Linux anteponiendo el script con el comando "sh", por ejemplo.sh example.sh.

2. Análisis de empalme alternativo (AS) utilizando RNA-seq

1. Descarga y preprocesamiento de datos
   NOTA: Los fragmentos de código anotados a continuación están disponibles en el archivo de código suplementario "AS_analysis_RNASeq.rmd", para seguir los pasos individuales de forma interactiva, y también se proporcionan como un script bash para ejecutarse por lotes en la línea de comandos de Linux (sh downloading_data_preprocessing.sh).
   1. Descarga de los datos sin procesar.
      1. Descargue los datos sin procesar de Sequence Read Archive (SRA) utilizando el comando 'prefetch' del kit de herramientas SRA (v2.10.8)¹⁷. Proporcione los ID de acceso de SRA en secuencia en el siguiente comando para descargarlos en paralelo utilizando la utilidad paralela^{GNU 18}. Para descargar archivos SRA de identificadores de acceso de SRR10261601 a SRR10261606 en paralelo, use lo siguiente en la línea de comandos de Linux.
         ​seq 10261601 10261606 | parallel prefetch SRR{}
      2. Extraiga los archivos fastq del archivo utilizando la función 'fastq-dump' del kit de herramientas SRA. Utilice GNU paralelo y dé los nombres de todos los archivos SRA juntos.
         ​parallel -j 3 fastq-dump --gzip --skip-technical --read-filter pass --dumpbase --split-e --clip --origfmt {} :::
      3. Descargue el genoma de referencia y las anotaciones para Mouse (ensamblaje del genoma GRCm39) desde www.ensembl.org utilizando lo siguiente en la línea de comandos de Linux.
         wget -nv -O annotation.gtf.gz http://ftp.ensembl.org/pub/release-103/gtf/mus_musculus/Mus_musculus.GRCm39.103.gtf.gz \ && gunzip -f annotation.gtf.gz
wget -nv -O genome.fa.gz http://ftp.ensembl.org/pub/release-103/fasta/mus_musculus/dna/Mus_musculus.GRCm39.dna.primary_assembly.fa.gz \ && gunzip -f genome.fa.gz
GTF=$(readlink -f annotation.gtf)
GENOME=$(readlink -f genome.fa)
   2. Preprocesamiento y mapeo de lecturas al ensamblaje del genoma
      1. Control de calidad. Evalúe la calidad de las lecturas sin procesar con FASTQC (FASTQ Quality Check v0.11.9)¹⁹. Cree una carpeta de salida y ejecute fastqc con paralelo en varios archivos fasta de entrada. Este paso generará un informe de calidad para cada muestra. Examine los informes para asegurarse de que la calidad de las lecturas es aceptable antes de realizar un análisis adicional. (Consulte el manual del usuario para comprender los informes en https://www.bioinformatics.babraham.ac.uk/projects/fastqc/)
         mkdir fastqc_out
parallel "fastqc {} -o fastqc_out" ::: $RAW_DATA/*.fastq.gz
         NOTA: Si es necesario, realice el recorte del adaptador con 'cutadapt'20 o 'trimmomatic'21 para eliminar la secuenciación en los adaptadores de flanqueo, que varía según el tamaño del fragmento de ARN y la longitud de lectura. En este análisis nos saltamos este paso ya que la fracción de lecturas afectadas fue mínima.
      2. Leer alineación. El siguiente paso en el preprocesamiento incluye el mapeo de las lecturas al genoma de referencia. En primer lugar, construya el índice para el genoma de referencia utilizando la función 'genomeGenerate' de STAR²² y luego alinee las lecturas sin procesar con la referencia (alternativamente, los índices preconstruidos están disponibles en el sitio web de STAR y se pueden usar directamente para la alineación). Ejecute los siguientes comandos en la línea de comandos de Linux.
         #Build STAR index
GDIR=STAR_indices
mkdir $GDIR
STAR --runMode genomeGenerate --genomeFastaFiles $GENOME --sjdbGTFfile $GTF --runThreadN 8 --genomeDir $GDIR
ODIR=results/mapping
mkdir -p $ODIR
#Align reads to the genome
for fq1 in $RAW_DATA/*R1.fastq.gz;
do
fq2=$(echo $fq1| sed 's/1.fastq.gz/2.fastq.gz/g');
OUTPUT=$(basename ${fq1}| sed 's/R1.fastq.gz//g');
STAR --genomeDir $GDIR \
--runThreadN 12 \
--readFilesCommand zcat \
--readFilesIn ${fq1}${fq2}\
--outFileNamePrefix $ODIR\/${OUTPUT} \
--outSAMtype BAM SortedByCoordinate \
--outSAMunmapped Within \
--outSAMattributes Standard
Done
         NOTA: El alineador STAR generará y ordenará archivos BAM (mapa de alineación binaria) para cada muestra después de la alineación de lectura. Los archivos Bam deben ordenarse antes de continuar con los pasos adicionales.
2. Preparación de anotaciones Exon.
   1. Ejecute el archivo de código complementario "prepare_mm_exon_annotation. R" con la anotación descargada en formato GTF (formato de transferencia génica) para preparar las anotaciones. Para ejecutarlo, escriba lo siguiente en la línea de comandos de Linux.
      Rscript prepare_mm_exon_annotation.R annotation.gtf
      NOTA: El archivo GTF contiene varias entradas de exón para diferentes isoformas. Este archivo se utiliza para "contraer" los múltiples ID de transcripción para cada exón. Es un paso importante definir los contenedores de conteo de exones.
3. El conteo se lee. El siguiente paso es contar el número de lecturas asignadas a diferentes transcripciones/exones. Ver archivo complementario: "AS_analysis_RNASeq.Rmd".
   1. Cargar las bibliotecas necesarias:
      packages<- c("Rsubread","tidyverse", "magrittr", "EnhancedVolcano", "edgeR","openxlsx")
invisible(lapply(packages, library, character.only=TRUE))
   2. Cargue el archivo de anotación procesado obtenido del paso anterior (2.2).
      load("mm_exon_anno.RData")
   3. Lea todos los archivos bam obtenidos en el paso 2.2.2 como entrada para que 'featureCounts' cuente las lecturas. Lea la carpeta que contiene los archivos bam enumerando primero cada archivo del directorio que termina con .bam. Utilice 'featureCounts' del paquete Rsubread que toma archivos bam y anotación GTF procesada (referencia) como entrada para generar una matriz de recuentos asociados con cada característica con filas que representan exones (características) y columnas que representan muestras.
      countData <- dir("bams", pattern=".bam$", full.names=T) %>%
featureCounts(annot.ext=anno,
isGTFAnnotationFile=FALSE,
minMQS=0,useMetaFeatures=FALSE,
allowMultiOverlap=TRUE,
largestOverlap=TRUE,
countMultiMappingReads=FALSE,
primaryOnly=TRUE,
isPairedEnd=TRUE,
nthreads=12)
   4. A continuación, realice un filtrado no específico para eliminar exones poco expresados ("no específico" indica que la información de la condición experimental no se utiliza en el filtrado, para evitar sesgos de selección). Transforme los datos de escala sin procesar a recuentos por millón (cpm) utilizando la función cpm del paquete 'edgeR'²³ y mantenga exones con recuentos superiores a un umbral configurable (para este conjunto de datos se usa un cpm) en al menos tres muestras. También eliminan genes con un solo exón.
      # Non-specific filtering: Remove the exons with low counts
isexpr<- rownames(countData$counts)[rowSums(cpm(countData$counts)>1) >=3]
countData$counts<-countData$counts[rownames(countData$counts) %in%isexpr, ]
anno<-anno%>% filter(GeneID%in% rownames(countData$counts))
# Remove genes with only 1 site and NA in geneIDs
dn<-anno%>%group_by(GeneID)%>%summarise(nsites=n())%>% filter(nsites>1&!is.na(GeneID))
anno<-anno%>% filter(GeneID%in%dn$GeneID)
countData$counts<-countData$counts[rownames(countData$counts) %in%anno$GeneID, ]
      NOTA: Compruebe los parámetros necesarios para featureCounts cuando utilice datos diferentes, por ejemplo, para lecturas de un solo extremo, establezca 'isPairedEnd = FALSE'. Consulte la guía del usuario de RSubread para elegir las opciones para sus datos y consulte la sección Discusión a continuación.
4. Empalme diferencial y análisis de uso de exones. Describimos dos alternativas para este paso: DEXSeq y DiffSplice. Cualquiera de los dos puede ser utilizado y dar resultados similares. Para mantener la coherencia, seleccione DEXSeq si prefiere un paquete DESeq2 para DGE y utilice DiffSplice para un análisis DGE basado en Limma. Ver ficha complementaria: "AS_analysis_RNASeq.rmd".
   1. Uso del paquete DEXSeq para el análisis diferencial de exones.
      1. Cargue la biblioteca y cree una tabla de ejemplo para definir el diseño experimental.
         library(DEXSeq)
sampleTable<-data.frame(row.names= c("Mbnl1KO_Thymus_1", "Mbnl1KO_Thymus_2", "Mbnl1KO_Thymus_3", "WT_Thymus_1", "WT_Thymus_2", "WT_Thymus_3"), condition= rep(c("Mbnl1_KO", "WT"),c(3,3)), libType= rep(c("paired-end")))
         Nota : los nombres de fila deben ser coherentes con los nombres de archivo bam utilizados por featureCounts para contar las lecturas. sampleTable consta de detalles de cada muestra que incluye: tipo de biblioteca y condición. Esto es necesario para definir los contrastes o el grupo de prueba para detectar el uso diferencial.
      2. Prepare el archivo de información de exón. La información de exón en forma de objetos GRanges (rangos genómicos) (https://bioconductor.org/packages/release/bioc/html/GenomicRanges.html) se requiere como entrada para crear el objeto DEXSeq en el siguiente paso. Haga coincidir los identificadores de genes con los recuentos de lectura para crear el objeto exoninfo.
         exoninfo<-anno[anno$GeneID%in% rownames(countData$counts),]
exoninfo<-GRanges(seqnames=anno$Chr,
ranges=IRanges(start=anno$Start, end=anno$End, width=anno$Width),strand=Rle(anno$Strand))
mcols(exoninfo)$TranscriptIDs<-anno$TranscriptIDs
mcols(exoninfo)$Ticker<-anno$Ticker
mcols(exoninfo)$ExonID<-anno$ExonID
mcols(exoninfo)$n<-anno$n
mcols(exoninfo)$GeneID<-anno$GeneID
transcripts_l= strsplit(exoninfo$TranscriptIDs, "\\,")
save(countData, sampleTable, exoninfo, transcripts_l, file="AS_countdata.RData")
      3. Crear objeto DEXSeq mediante la función DEXSeqDataSet. El objeto DEXSeq recopila recuentos de lectura, información de características de exón e información de ejemplo. Utilice los recuentos de lectura generados en el paso 3 y la información de exón obtenida del paso anterior para crear el objeto DEXSeq a partir de la matriz de recuento. El argumento sampleData toma una entrada de marco de datos que define las muestras (y sus atributos: tipo de biblioteca y condición), 'design' usa sampleData para generar una matriz de diseño para la prueba diferencial utilizando la notación de fórmula de modelo. Tenga en cuenta que un término de interacción significativo, condición: exón, indica que la fracción de lecturas sobre un gen que cae en un exón particular depende de la condición experimental, es decir, hay AS. Consulte la documentación de DEXSeq para obtener una descripción completa de la configuración de la fórmula del modelo para diseños experimentales más complejos. Para obtener información sobre las características, se requieren identificadores de exón, genes correspondientes y transcripciones.
         ​dxd<-DEXSeqDataSet(countData$counts,sampleData=sampleTable, design=~sample+exon+condition:exon,featureID=exoninfo$ExonID,groupID=exoninfo$GeneID,featureRanges=exoninfo, transcripts=transcripts_l)
      4. Normalización y estimación de dispersión. A continuación, realice la normalización entre muestras y estime la varianza de los datos, debido tanto al ruido de conteo de Poisson de la naturaleza discreta de RNA-seq como a la variabilidad biológica, utilizando los siguientes comandos.
         ​dxd %<>% estimateSizeFactors %>% estimateDispersions %T>% plotDispEsts
      5. Prueba de uso diferencial. Después de la estimación de la variación, pruebe el uso diferencial del exón para cada gen y genere los resultados.
         dxd%<>%testForDEU%>%estimateExonFoldChanges(fitExpToVar=
"condition")#Estimate fold changes
dxr=DEXSeqResults(dxd)
      6. Visualice los eventos de empalme para los genes seleccionados mediante el siguiente comando.
         plotDEXSeq(dxr,"Wnk1", displayTranscripts=TRUE, splicing=TRUE,legend
=TRUE,cex.axis=1.2,cex=1.3,lwd=2)
         Examine el archivo de R Notebook "AS_analysis_RNASeq.Rmd" para generar gráficos adicionales para genes de interés y para generar gráficos de volcanes en diferentes umbrales.
   2. Uso de diffSplice de Limma para identificar el empalme diferencial. Siga el archivo de R Notebook "AS_analysis_RNASeq.rmd". Asegúrese de que se han seguido los pasos 2.1-2.3 para preparar los archivos de entrada antes de continuar.
      1. Cargar bibliotecas
         library(limma)
library(edgeR)
      2. Filtrado no específico. Extraer la matriz de recuentos de lectura obtenidos en 2.3. Cree una lista de entidades utilizando la función 'DGEList' del paquete edgeR, donde las filas representan genes y las columnas representan muestras.
         mycounts=countData$counts
#Change the rownames of the countdata to exon Ids instead of genes for unique rownames.
rownames(mycounts) = exoninfo$ExonID
dge<-DGEList(counts=mycounts)
#Filtering
isexpr<- rowSums(cpm(dge)>1) >=3
dge<-dge[isexpr,,keep.lib.sizes=FALSE]
#Extract the exon annotations for only transcripts meeting non-specific filter
exoninfo=anno%>% filter(ExonID%in% rownames(dge$counts))
#Convert the exoninfo into GRanges object
exoninfo1<-GRanges(seqnames=exoninfo$Chr,
ranges=IRanges(start=exoninfo$Start, end=exoninfo$End, width=exoninfo$Width),strand=Rle(exoninfo$Strand))
mcols(exoninfo1)$TranscriptIDs<-exoninfo$TranscriptIDs
mcols(exoninfo1)$Ticker<-exoninfo$Ticker
mcols(exoninfo1)$ExonID<-exoninfo$ExonID
mcols(exoninfo1)$n<-exoninfo$n
mcols(exoninfo1)$GeneID<-exoninfo$GeneID
transcripts_l= strsplit(exoninfo1$TranscriptIDs, "\\,")
         NOTA: Como paso de filtrado no específico, los recuentos se filtran por cpm < 1 en x de n muestras, donde x es el número mínimo de réplicas en cualquier condición. n = 6 y x = 3 para estos datos de ejemplo.
      3. Normalice los recuentos entre muestras, con la función 'calcNormFactors' del paquete 'edgeR' utilizando la media recortada de los valores M (método de normalización TMM)²⁴ Calculará factores de escala para ajustar los tamaños de la biblioteca.
         ​dge<-calcNormFactors(dge)
      4. Utilice sampleTable como se generó en el paso 2.4.1.1 y cree la matriz de diseño. La matriz de diseño caracteriza el diseño. Consulte los capítulos 8 y 9 de la Guía del usuario de Limma (https://www.bioconductor.org/packages/devel/bioc/vignettes/limma/inst/doc/usersguide.pdf) para obtener detalles sobre las matrices de diseño para diseños experimentales más avanzados.
         Treat<- factor(sampleTable$condition)
design<- model.matrix(~0+Treat)
colnames(design) <- levels(Treat)
      5. Ajuste un modelo lineal por exón. Ejecute la función 'voom' del paquete 'limma' para procesar datos de RNA-seq para estimar la varianza y generar pesos de precisión para corregir el ruido del recuento de Poisson, y transformar los recuentos a nivel de exón en log2 recuentos por millón (logCPM). A continuación, ejecute el modelado lineal utilizando la función 'lmfit' para ajustar los modelos lineales a los datos de expresión de cada exón. Calcule estadísticas empíricas de Bayes para el modelo ajustado utilizando la función 'eBayes' para detectar la expresión diferencial del exón. A continuación, defina una matriz de contraste para las comparaciones experimentales de interés. Utilice 'contrasts.fit' para obtener coeficientes y errores estándar para cada par de comparación.
         v<-voom(dge,design,plot=FALSE)
fit<-lmFit(v,design)
fit<-eBayes(fit)
colnames(fit)
cont.matrix<-makeContrasts(
Mbnl1_KO_WT=Mbnl1_KO-WT,
levels=design)
fit2<-contrasts.fit(fit,cont.matrix)
      6. Análisis de empalme diferencial. Ejecute 'diffSplice' en el modelo ajustado para probar las diferencias en el uso de exones de genes entre wild-type y knockout y explore los resultados mejor clasificados usando la función 'topSplice': test = "t" da una clasificación de los exones AS, test = "simes" da una clasificación de genes.
         ex<-diffSplice(fit2,geneid=exoninfo$GeneID,exonid=exoninfo$ExonID)
ts<-topSplice(ex,n=Inf,FDR=0.1, test="t", sort.by="logFC")
​tg<-topSplice(ex,n=Inf,FDR=0.1, test="simes")
      7. Visualización. Trazar los resultados con la función 'plotSplice', dando al gen de interés en el argumento geneid. Guarde los resultados principales ordenados por registro Pliegue el cambio en un objeto y genere un diagrama de volcán para exhibir los exones.
         plotSplice(ex,geneid="Wnk1", FDR=0.1)
#Volcano plot
EnhancedVolcano(ts,lab=ts$ExonID,selectLab= head((ts$ExonID),2000), xlab= bquote(~Log[2]~'fold change'), x='logFC', y='P.Value', title='Volcano Plot', subtitle='Mbnl1_KO vs WT (Limma_diffSplice)', FCcutoff=2, labSize=4,legendPosition="right", caption= bquote(~Log[2]~"Fold change cutoff, 2; FDR 10%"))
   3. Uso de rMATS
      1. Asegúrese de que la última versión de rMATS v4.1.1 (también conocida como rMATS turbo debido al tiempo de procesamiento reducido y menos requisitos de memoria) se instala utilizando conda o github (https://github.com/Xinglab/rmats-turbo/releases/download/v4.1.1/rmats_turbo_v4_1_1.tar.gz) en el directorio de trabajo. Siga la sección 4.3 en "AS_analysis_RNASeq.Rmd".
      2. Vaya a la carpeta que contiene los archivos bam obtenidos después de asignar y prepare los archivos de texto, según lo requiera rMATS, para las dos condiciones copiando el nombre de los archivos bam (junto con la ruta) separados por ','. Los siguientes comandos deben ejecutarse en la línea de comandos de Linux:
         mkdir rMATS_analysis
cd bams/
ls -pd "$PWD"/*| grep "WT"| tr '\n'','> Wt.txt
ls -pd "$PWD"/*| grep "Mb"| tr '\n'','> KO.txt
mv *.txt ../rMATS_analysis
      3. Ejecute rmats.py con los dos archivos de entrada generados en el paso anterior, junto con el archivo GTF obtenido en 2.1.1.3. Esto generará una carpeta de salida 'rmats_out' que contiene archivos de texto que describen estadísticas (valores p y niveles de inclusión) para cada evento de empalme por separado.
         python rmats-turbo/rmats.py --b1 KO.txt --b2 Wt.txt --gtf annotation.gtf -t paired --readLength 50 --nthread 8 --od rmats_out/ --tmp rmats_tmp --task pos
         NOTA: La anotación de referencia en forma de archivo GTF también es necesaria. Compruebe los parámetros si los datos son de un solo extremo y cambie la opción -t en consecuencia.
      4. Explorando los resultados de rMATS. Utilice el paquete Bioconductor '^maser'25 para explorar los resultados de rMATS. Cargue los archivos de texto de conteo de uniones y exones (JCEC) en el objeto 'maser' y filtre el resultado en función de la cobertura incluyendo al menos cinco lecturas promedio por evento de empalme.
         library(maser)
mbnl1<-maser("/rmats_out/", c("WT","Mbnl1_KO"), ftype="JCEC")
#Filtering out events by coverage
mbnl1_filt<-filterByCoverage(mbnl1,avg_reads=5)
      5. Visualización de los resultados de rMATS. Seleccione los eventos de empalme significativos a la tasa de descubrimiento falso (FDR) del 10% y el cambio mínimo del 10% en el porcentaje empalmado (deltaPSI) utilizando la función 'topEvents' del paquete 'maser'. A continuación, verifique los eventos genéticos para genes individuales de interés (gen de muestra-Wnk1) y trace los valores de PSI para cada evento de empalme de ese gen. Genere un diagrama de volcán especificando el tipo de evento.
         #Top splicing events at 10% FDR
mbnl1_top<-topEvents(mbnl1_filt,fdr=0.1, deltaPSI=0.1)
mbnl1_top
#Check the gene events for a particular gene
mbnl1_wnk1<-geneEvents(mbnl1_filt,geneS="Wnk1", fdr=0.1, deltaPSI=0.1)
maser::display(mbnl1_wnk1,"SE")
plotGenePSI(mbnl1_wnk1,type="SE", show_replicates
=TRUE)
​volcano(mbnl1_filt,fdr=0.1, deltaPSI=0.1,type="SE")
+xlab("deltaPSI")+ylab("Log10 Adj. Pvalue")+ggtitle("Volcano Plot of exon skipping events")
      6. Genere gráficos de Sashimi para el resultado de eventos de empalme obtenidos con rMATS en forma de archivos de texto utilizando el paquete 'rmats2shahimiplot'. Ejecute el script de Python en la línea de comandos de Linux.
         python ./src/rmats2sashimiplot/rmats2sashimiplot.py --b1 ../bams/WT_Thymus_1.bam,../bams/WT_Thymus_2.bam,../bams/WT_Thymus_3.bam --b2 ../bams/Mbnl1KO_Thymus_1.bam,../bams/Mbnl1KO_Thymus_2.bam,../bams/Mbnl1KO_Thymus_3.bam -t SE -e ../rMATS_analysis/rmats_out/SE.MATS.JC.txt --l1 WT --l2 Mbnl1_KO --exon_s 1 --intron_s 5 -o ../rMATS_analysis/rmats2shasmi_output
         NOTA: Este proceso puede llevar mucho tiempo, ya que generará el gráfico de Sashimi para todos los resultados en el archivo de eventos. Elija los resultados principales (nombres de genes y exones) como se muestra en la función topEvents de 'maser' y visualice el gráfico de Sashimi correspondiente.

3. Análisis alternativo de poliadenilación (APA) utilizando secuenciación de extremo 3'

1. Descarga y preprocesamiento de datos
   NOTA: Consulte el archivo suplementario de R Notebook "APA_analysis_3PSeq_notebook. Rmd" para los comandos completos para los pasos de descarga y preprocesamiento de datos, o ejecute el archivo bash suplementario "APA_data_downloading_preprocessing.sh" en la línea de comandos de Linux.
   1. Descargar datos de SRA con los ID de Acceso (1553129 a 1553136).
   2. Adaptadores de recorte y complemento inverso para obtener la secuencia de hebras de detección.
      NOTA: Este paso es específico para el ensayo PolyA-seq utilizado.
   3. Lecturas de mapa para el ensamblaje del genoma del ratón usando el alineador de pajarita²⁶.
2. Preparación de anotaciones de sitios pA.
   NOTA: El procesamiento del archivo de anotación del sitio pA se realiza en primer lugar utilizando el archivo suplementario de R Notebook "APA_analysis_3PSeq_notebook. Rmd" (2.1 - 2.6), y luego usando el archivo bash "APA_annotation_preparation.sh".
   1. Descargar anotación de sitios pA de la base de datos PolyASite 2.0⁶.
   2. Seleccione anotaciones de sitio pA para retener sitios pA de región no traducida (UTR) de 3', que se anotan como Exón terminal (TE) o 1000 nt aguas abajo de un exón terminal anotado (DS) para el análisis descendente.
   3. Obtener picos de sitio pA. Ancla en cada sitio de escisión pA y visualiza la cobertura de lectura promedio usando herramientas de cama y herramientas profundas^27,28. Los resultados mostraron que los picos de las lecturas mapeadas se dispersaron principalmente dentro de ~ 60 pb aguas arriba de los sitios de escisión (figura 5 y figura complementaria 5). Por lo tanto, las coordenadas de los sitios pA se extendieron desde el archivo de anotación a 60 pb aguas arriba de sus sitios de escisión. Dependiendo del protocolo específico de secuenciación final de 3' utilizado, este paso deberá optimizarse para ensayos que no sean PolyA-seq.
3. Lecturas de conteo
   1. Prepare el archivo de anotación de sitios pA.
      anno<- read.table(file= "flanking60added.pA_annotation.bed",
stringsAsFactors=FALSE, check.names=FALSE, header=FALSE, sep="")
colnames(anno) <- c("chrom", "chromStart", "chromEnd", "name", "score", "strand", "rep", "annotation", "gene_name", "gene_id")
anno<- dplyr::select(anno,name,chrom, chromStart,chromEnd, strand,gene_id,gene_name,rep)
colnames(anno) <- c("GeneID", "Chr", "Start", "End", "Strand", "Ensembl", "Symbol", "repID")
   2. Aplique 'featureCounts' para adquirir recuentos sin procesar. Guarde la tabla de recuento como el archivo RData "APA_countData.Rdata" para el análisis APA utilizando diferentes herramientas.
      countData<- dir("bamfiles", pattern="sorted.bam$", full.names=TRUE) %>%
# Read all bam files as input for featureCounts
featureCounts(annot.ext=anno, isGTFAnnotationFile= FALSE,minMQS=0,useMetaFeatures= TRUE,allowMultiOverlap=TRUE, largestOverlap= TRUE,strandSpecific=1, countMultiMappingReads =TRUE,primaryOnly= TRUE,isPairedEnd= FALSE,nthreads=12)%T>%
save(file="APA_countData.Rdata")
      NOTA: Tenga en cuenta que debe cambiar cualquiera de los parámetros enumerados en la función 'featureCounts'. Modifique el parámetro 'strandSpecific' para asegurarse de que sea consistente con la dirección de secuenciación del ensayo de secuenciación final 3' utilizado (empíricamente, visualizar los datos en un navegador del genoma sobre los genes en las hebras más y menos aclarará esto).
   3. Aplicar filtrado no específico de countData. El filtrado puede mejorar significativamente la solidez estadística en las pruebas de uso diferencial del sitio pA. Primero, eliminamos esos genes con un solo sitio pA, en el que no se puede definir diferencialmente el uso del sitio pA. En segundo lugar, aplicamos un filtrado no específico basado en la cobertura: los recuentos se filtran por cpm menos de 1 en x de n muestras, donde x es el número mínimo de réplicas en cualquier condición. N = 8 y x = 2 para estos datos de ejemplo.
      load(file= "APA_countData.Rdata")# Skip this step if already loaded
# Non-specific filtering: Remove the pA sites not differentially expressed in the samples
countData<-countData$counts%>%as.data.frame%>% .[rowSums(edgeR::cpm(.)>1) >=2, ]
anno%<>% .[.$GeneID%in% rownames(countData), ]
# Remove genes with only 1 site and NA in geneIDs
dnsites<-anno%>%group_by(Symbol)%>%summarise(nsites=n())%>% filter(nsites>1&!is.na(Symbol))
anno<-anno%>% filter(Symbol%in%dnsites$Symbol)
countData<-countData[rownames(countData) %in%anno$GeneID, ]
4. Análisis diferencial de uso de sitios de poliadenilación utilizando tuberías DEXSeq y diffSplice.
   1. Uso del paquete DEXSeq
      NOTA: Como no se puede definir una matriz de contraste para la tubería DEXSeq, el análisis APA diferencial de cada dos condiciones experimentales debe realizarse por separado. El análisis diferencial APA de la condición WT y la condición DKO se realiza como ejemplo para explicar el procedimiento. Consulte el archivo complementario "APA_analysis_3PSeq_notebook. RMD" para el flujo de trabajo paso a paso de esta sección y el análisis diferencial APA de otros contrastes.
      1. Cargue la biblioteca y cree una tabla de ejemplo para definir el diseño experimental.
         c("DEXSeq", "GenomicRanges") %>% lapply(library, character.only=TRUE) %>%invisible
sampleTable1<- data.frame(row.names= c("WT_1","WT_2","DKO_1","DKO_2"),
condition= c(rep("WT", 2), rep("DKO", 2)),
​libType= rep("single-end", 4))
      2. Prepare el archivo de información de sitios pA utilizando el paquete Bioconductor GRanges.
         # Prepare the GRanges object for DEXSeqDataSet object construction
PASinfo <- GRanges(seqnames = anno$Chr,
ranges = IRanges(start = anno$Start, end = anno$End),strand = Rle(anno$Strand))
mcols(PASinfo)$PASID<-anno$repID
mcols(PASinfo)$GeneEns<-anno$Ensembl
mcols(PASinfo)$GeneID<-anno$Symbol
# Prepare the new feature IDs, replace the strand information with letters to match the current pA site clusterID
new.featureID <- anno$Strand %>% as.character %>% replace(. %in% "+", "F") %>% replace(. %in% "-", "R") %>% paste0(as.character(anno$repID), .)
      3. Utilice los recuentos de lectura generados en el paso 3.3 y la información del sitio pA obtenida del paso anterior para crear el objeto DEXSeq.
         # Select the read counts of the condition WT and DKO
countData1<- dplyr::select(countData, SRR1553129.sorted.bam, SRR1553130.sorted.bam, SRR1553131.sorted.bam, SRR1553132.sorted.bam)
# Rename the columns of countData using sample names in sampleTable
colnames(countData1) <- rownames(sampleTable1)
dxd1<-DEXSeqDataSet(countData=countData1,
sampleData=sampleTable1,
design=~sample+exon+condition:exon,
featureID=new.featureID,
groupID=anno$Symbol,
featureRanges=PASinfo)
      4. Defina el par de contraste definiendo los niveles de condiciones en el objeto DEXSeq.
         dxd1$condition<- factor(dxd1$condition, levels= c("WT", "DKO"))
# The contrast pair will be "DKO - WT"
      5. Normalización y estimación de dispersión. Similar a los datos de RNA-seq, para los datos de secuenciación final de 3' realice la normalización entre muestras (mediana de proporciones en columna para cada muestra) utilizando la función 'estimateSizeFactors', y estime la variación de los datos usando con la función 'estimateDispersions', luego visualice el resultado de la estimación de dispersión usando la función 'plotDispEsts'.
         ​dxd1 %<>% estimateSizeFactors %>% estimateDispersions %T>% plotDispEsts
      6. Prueba diferencial de uso del sitio pA para cada gen usando la función 'testForDEU', luego estime el cambio de pliegue del uso del sitio pA usando la función 'estimateExonFoldChanges'. Verifique los resultados utilizando la función 'DEXSeqResults' y establezca 'FDR < 10%' como criterio para sitios de pA significativamente diferenciales.
         dxd1 %<>% testForDEU %>% estimateExonFoldChanges(fitExpToVar = "condition")
dxr1 <- DEXSeqResults(dxd1)
dxr1
mcols(dxr1)$description
table(dxr1$padj<0.1) # Check the number of differential pA sites (FDR < 0.1)
table(tapply(dxr1$padj<0.1, dxr1$groupID, any)) # Check the number of gene overlapped with differential pA site
      7. La visualización de los resultados de uso diferencial del sitio pA utilizando gráficos APA diferenciales generados por la función 'plotDEXSeq' y gráficos de volcán por la función 'EnhancedVolcano'.
         # Select the top 100 significant differential pA sites ranked by FDR
topdiff.PAS<- dxr1%>%as.data.frame%>%rownames_to_column%>%arrange(padj)%$%groupID[1:100]

# Apply plotDEXSeq for the visualization of differential polyA usage
plotDEXSeq(dxr1,"S100a7a", legend=TRUE, expression=FALSE,splicing=TRUE, cex.axis=1.2, cex=1.3,lwd=2)

# Apply perGeneQValue to check the top genes with differential polyA site usage
dxr1%<>% .[!is.na(.$padj), ]
dgene<- data.frame(perGeneQValue= perGeneQValue (dxr1)) %>%rownames_to_column("groupID")

dePAS_sig1<-dxr1%>% data.frame() %>%
dplyr::select(-matches("dispersion|stat|countData|genomicData"))%>%
inner_join(dgene)%>%arrange(perGeneQValue)%>%distinct()%>%
filter(padj<0.1)

# Apply EnhancedVolcano package to visualise differential polyA site usage
"EnhancedVolcano"%>% lapply(library, character.only=TRUE) %>%invisible
EnhancedVolcano(dePAS_sig1, lab=dePAS_sig1$groupID, x='log2fold_DKO_WT',
y='pvalue',title='Volcano Plot',subtitle='DKO vs WT',
FCcutoff=1,labSize=4, legendPosition="right",
caption= bquote(~Log[2]~"Fold change cutoff, 1; FDR 10%"))
   2. Usando el paquete diffSplice. Consulte el archivo complementario de R Notebook "APA_analysis_3PSeq_notebook. Rmd" para el flujo de trabajo paso a paso de esta sección.
      1. Definir contrastes de interés para el análisis diferencial de uso de pA.
         Nota : este paso debe realizarse después de la construcción y el procesamiento del objeto DGEList, que se incluye en el archivo de R Notebook "APA_analysis_3PSeq_notebook. Rmd".
         contrast.matrix<-makeContrasts(DKO_vs_WT=DKO-
WT,Ctrl_vs_DKO=Ctrl-DKO,
KD_vs_Ctrl=KD-Ctrl,KD_vs_DKO=KD-DKO,levels=design)
fit2<-fit%>%contrasts.fit(contrast.matrix)%>%eBayes
summary(decideTests(fit2))
ex<-diffSplice(fit2,geneid=anno$Symbol,exonid=new.featureID)
topSplice(ex) #Check the top significant results with topSplice
      2. Visualice el resultado de contrastes de interés (aquí "DKO - WT") utilizando gráficos APA diferenciales mediante la función 'plotSplice' y gráficos de volcanes con la función 'EnhancedVolcano'. Consulte el archivo de R Notebook "APA_analysis_3PSeq_notebook. Rmd" 4.2.7 - 4.2.9 para la visualización de otros pares de contraste.
         sig1<-topSplice(ex,n=Inf,FDR=0.1,coef=1, test="t", sort.by="logFC")
sig1.genes<-topSplice(ex,n=Inf,FDR=0.1,coef=1, test="simes")
plotSplice(ex, coef=1,geneid="S100a7a", FDR = 0.1)
plotSplice(ex,coef=1,geneid="Tpm1", FDR = 0.1)
plotSplice(ex,coef=1,geneid="Smc6", FDR = 0.1)
EnhancedVolcano(sig1, lab=sig1$GeneID,xlab= bquote(~Log[2]~'fold change'),
x='logFC', y='P.Value', title='Volcano Plot', subtitle='DKO vs WT',
FCcutoff=1, labSize=6, legendPosition="right")

Representative Results

Después de ejecutar el flujo de trabajo paso a paso anterior, los resultados del análisis AS y APA y los resultados representativos se presentan en forma de tablas y gráficos de datos, generados de la siguiente manera.

COMO:
El resultado principal del análisis AS (Tabla suplementaria 1 para diffSplice; Tabla 2 para DEXSeq) es una lista de exones que muestran el uso diferencial entre condiciones, y una lista de genes que muestran una actividad de empalme general significativa de uno o más de sus exones constituyentes, clasificados por significación estadística. La tabla complementaria 1, pestaña 2 muestra exones significativos, con columnas que muestran CF diferencial de exón versus reposo, valor p no ajustado por exón y valor p ajustado (corrección de Benjamini-Hockberg). El umbral en los valores p ajustados dará un conjunto de exones con FDR definido. La Tabla complementaria 1, pestaña 3 muestra una lista clasificada de genes que muestran la importancia de la actividad general de empalme, con una columna que muestra el valor p ajustado a nivel de gen calculado utilizando el método Simes. Datos similares se muestran en la Tabla 2 para DEXSeq. La Figura Suplementaria 1 y la Figura Suplementaria 2 muestran un patrón de empalme diferencial en los genes Mbnl1, Tcf7 y Lef1 que han sido validados experimentalmente en el artículo publicado presentado con los datos¹⁵. Los autores han demostrado la validación experimental de cinco genes: Mbnl1, Mbnl2, Lef1, Tcf7 y Ncor2. Nuestro enfoque detectó patones de empalme diferencial en todos estos genes. Aquí presentamos los niveles de FDR para cada gen utilizando DEXSeq, diffSplice y rMATS respectivamente obtenidos en las Tablas Suplementarias 1-3: Mbnl1 (0, 6.6E-61,0), Mbnl2 (0,0.18,0), Lef1 (1.4E-10, 1.3E-04, 0), Tcf7 (0, 1.1E-6, 0) y Ncor2 (9.2E-11, 1.6E-06, 0).

La Figura 2 muestra una comparación entre los resultados obtenidos de tres herramientas diferentes e ilustra patrones de empalme alternativos en el gen Wnk1. Los diagramas de volcanes se muestran en la Figura 2A (diffSplice) y la Figura 2B (DEXSeq). Otros tres genes altamente clasificados se muestran en la Figura Suplementaria 1 (diffSplice) y la Figura Suplementaria 2 (DEXSeq). El eje y muestra la significación estadística (valores p -log10) y el eje x muestra el tamaño del efecto (cambio de pliegue). Los genes ubicados en los cuadrantes superiores izquierdo o derecho indican una CF sustancial y una fuerte evidencia estadística de diferencias genuinas.

Figura 2. Comparación de resultados de empalme alternativos obtenidos de diffSplice, DEXSeq y rMATS. |
(A) Diagrama de volcán (izquierda) de RNA-Seq de Limma diffAnálisis de empalme: El eje x muestra el cambio de pliegue del exón logarítmico; El eje Y muestra el valor p -log₁₀ . Cada punto corresponde a un exón. Línea discontinua horizontal en valor p=1E-5; líneas discontinuas verticales con doble cambio (FC). Los exones rojos muestran una FC sustancial y significación estadística. Gráfico de empalme diferencial (derecha): Los patrones de empalme se exhiben para un ejemplo del gen Wnk1 donde el eje x muestra el exón id por transcripción; el eje y muestra el cambio relativo del pliegue logarítmico del exón (la diferencia entre el logFC del exón y el logFC general para todos los demás exones). Los exones resaltados en rojo muestran una expresión diferencial estadísticamente significativa (FDR < 0,1).
(B) Diagrama de volcán (izquierda) y uso diferencial de exones (derecha) de RNA-Seq obtenido del análisis DEXSeq. El gen Wnk1 muestra un uso diferencial significativo de exones entre WT y Mbnl1 knock-out resaltados en rosa, que corresponden a los mismos exones diferenciales en (A).
(C) Diagrama de volcán (izquierda) y gráfica de Sashimi (derecha) para Wnk1 obtenida del análisis rMATS. Diagrama de volcán que representa un evento significativo de exón (SE) omitido (casete) en tipo salvaje en comparación con knockout al 10% FDR con cambio en los valores de porcentaje empalmado (PSI o ΔΨ) > 0.1. El eje x muestra el cambio en los valores de PSI en todas las condiciones, y el eje y muestra el valor p del registro. La gráfica de Sashimi muestra un evento de exón omitido en el gen Wnk1 , correspondiente a un exón diferencial significativo en (A) y (B). Cada fila representa una muestra de RNA-Seq: tres réplicas de tipo salvaje y knock-out de Mbnl1. La altura muestra la cobertura de lectura en RPKM y los arcos de conexión representan lecturas de unión a través de exones. Las isoformas alternativas del modelo genético anotado se muestran en la parte inferior de la gráfica. El panel inferior de C ilustra las lecturas de unión utilizadas para calcular la estadística PSI.
(D) Diagrama de Venn comparando el número de exones diferenciales significativos obtenidos por los diferentes métodos. Haga clic aquí para ver una versión más grande de esta figura.

Figura 2 Un (panel derecho) muestra una visualización diagramática de las diferencias de exón de uno de los genes mejor clasificados, mostrando logFC en el eje y y el número de exón en el eje x. Este ejemplo muestra un exón de casete que varía entre las condiciones para el gen Wnk1. La gráfica de uso de exones diferenciales de DEXSeq muestra evidencia de empalme diferencial en cinco sitios de exón cerca de Wnk1.6.45. Es probable que los exones resaltados en rosa se empalmen en muestras de Mbnl1 KO en comparación con WT. Estos exones son complementarios a los resultados obtenidos por diffSplice que muestra un patrón similar en la posición genómica específica. Más ejemplos se muestran en la Figura Suplementaria 1 y la Figura Suplementaria 2. Se puede dar una visión más detallada para confirmar resultados interesantes comparando las pistas de cobertura (ondulación) en unidades RPM (lecturas por millón) en los navegadores del genoma de la UCSC (Universidad de Santa Cruz) o IGV (Integrative Genomics Viewer) (no mostrado), junto con la correlación visual con otras pistas de interés, como modelos genéticos conocidos, conservación y otros ensayos de todo el genoma.

La tabla de salida de rMATS enumera los eventos de empalme alternativos significativos categorizados por tipo (Tabla complementaria 3). La Figura 2C muestra un diagrama de volcán de genes que son empalmados alternativamente, con el tamaño del efecto medido por la estadística diferencial "porcentaje empalmado" (PSI o ΔΨ) de¹¹.

PSI se refiere al porcentaje de lecturas consistentes con la inclusión de un exón de casete (es decir, lecturas que se asignan al exón del casete o lecturas de unión que se superponen al exón) en comparación con lecturas consistentes con la exclusión del exón, es decir, lecturas de unión a través de exones adyacentes aguas arriba y aguas abajo (El panel inferior de la Figura 2C). El panel derecho de la Figura 2C muestra el diagrama de sashimi del gen Wnk1 con un evento de empalme diferencial superpuesto en las pistas de cobertura para el gen, con un exón omitido en Mbnl1 KO. Los arcos que unen exones muestran el número de lecturas de unión (lecturas que cruzan un intrón empalmado). Las diferentes pestañas de la Tabla Suplementaria 3 muestran lecturas significativas de cada tipo de evento que abarca uniones con límites de exón (recuentos de uniones y recuentos de exones (JCEC)). La Figura 2D compara los exones empalmados diferencialmente significativos detectados por las tres herramientas.

Figura 3. Eventos de empalme alternativos adquiridos por análisis rMATS. A) Tipos de eventos AS. Esta figura está adaptada de la documentación de rMATS¹¹ que explica el evento de empalme con exones constitutivos y, alternativamente, empalmados. Etiquetado con el número de cada evento en FDR 10%. B) Diagrama de sashimi del gen Add3 que exhibe exón omitido (SE). C) Sashimi plot del gen Baz2b que exhibe un sitio alternativo de empalme 5' (A5SS). D) Sashimi plot del gen Lsm14b que exhibe un sitio de empalme alternativo 3' (A3SS). E) Diagrama de sashimi del gen Mta1 que exhibe exones mutuamente excluyentes (MXE). F) Diagrama de sashimi del gen Arpp21 que exhibe intrón retenido (RI). Las filas rojas representan tres réplicas de tipo salvaje y las filas naranjas representan réplicas knock-out de Mbnl1. El eje x corresponde a coordenadas genómicas e información de hebras, el eje y muestra cobertura en RPKM. Haga clic aquí para ver una versión más grande de esta figura.

La Figura 3 ilustra los tipos de eventos de empalme SE, A5SS, A3SS, MXE y RI con la ayuda de gráficos de Sashimi de los principales genes significativos de esos eventos. Al comparar las tres réplicas de WT y Mbnl1_KO, se detectaron un total de 1272 eventos SE, 130 eventos A5SS, 116 A3SS, 215 MXE y 313 eventos RI en FDR 10%. El diagrama de Sashimi ilustra el tipo de evento usando los genes principales como ejemplo.

APA:
El resultado del análisis APA es similar al análisis AS a nivel de exón. Se proporciona una tabla de los principales genes clasificados por actividad diferencial APA en el 3'UTR (Tabla suplementaria 4 y Tabla complementaria 5). La Figura 4A muestra los gráficos volcánicos de genes por actividad APA diferencial en 3'UTRs generados usando diffSplice y DEXSeq por separado. La Figura 4B muestra el diagrama de Venn comparando los resultados de uso del sitio de pA significativamente diferenciales adquiridos de diferentes tuberías. Las Figuras 4C y 4D muestran la representación esquemática del uso diferencial del sitio pA en el 3'UTR del gen Fosl2 (Figura 4C) y Papola (Figura 4D) generado utilizando tanto diffSplice como DEXSeq, que se validan experimentalmente para mostrar un desplazamiento distal a proximal significativo (Fosl2) y un desplazamiento proximal a distal significativo (Papola) del uso del sitio pA en DKO, respectivamente. Se muestran más ejemplos en la Figura Suplementaria 3 y la Figura Suplementaria 4.

Figura 4. Parcelas alternativas de poliadenilación por diffSplice y DEXSeq. A) Gráficos volcánicos de datos PolyA-seq generados usando diffSplice y DEXSeq. El eje X muestra el cambio de pliegue del sitio log pA; El eje Y muestra el valor p -log₁₀ . Cada punto corresponde a un sitio pA. Línea discontinua horizontal en valor p=1E-5; líneas discontinuas verticales en 2 veces FC. Los exones rojos muestran una FC sustancial y significación estadística. B) Diagrama de Venn que compara los resultados diferenciales significativos de uso del sitio de pA adquiridos de diferentes tuberías. C-D) Gráficos APA diferenciales generados usando diffSplice y DEXSeq que muestran los sitios pA proximal, interno y distal para los genes Fosl2 y Papola . Las figuras se generan con la misma función que la Figura 2 (B) pero con sitios pA reemplazando exones. Haga clic aquí para ver una versión más grande de esta figura.

La Figura 5 es una gráfica de diagnóstico para confirmar la distribución de lectura esperada alrededor de los sitios de escisión de pA anotados para el ensayo PolyA-seq utilizado. Muestra la cobertura media en regiones flanqueantes ancladas en sitios de escisión pA conocidos a nivel de todo el genoma. En este caso, se visualiza la acumulación esperada de lecturas aguas arriba de los sitios. Las distribuciones de lectura ancladas en sitios pA para todas las muestras de PolyA-seq se muestran en la Figura complementaria 5.

Figura 5. Diagrama de cobertura media alrededor de los sitios de escisión de pA. El sitio de escisión para los datos PolyA-seq se muestra mediante una línea discontinua vertical. El eje X muestra la posición de la base en relación con los sitios de escisión de pA, hasta 100 nucleótidos aguas arriba y aguas abajo; El eje y muestra la cobertura media de lectura en todos los sitios de escisión de pA, normalizada por el tamaño de la biblioteca en CPM. Haga clic aquí para ver una versión más grande de esta figura.

Los resultados diferenciales APA de diferentes contrastes generados por la misma canalización se pueden comparar y verificar visualizando la cobertura de lectura de sitios pA representativos significativamente diferenciales en el navegador del genoma. La Figura 6A es el diagrama de Venn que compara el uso significativamente diferencial del sitio pA de diferentes contrastes adquiridos de diffSplice. La Figura 6B-D son las instantáneas IGV de la cobertura de lectura en los sitios pA para diferentes genes, que muestran los patrones consistentes con los descubiertos en el análisis APA usando diffSplice. La Figura 6B valida el desplazamiento significativo proximal a distal del uso del sitio pA para el gen Paip2, que se detecta de forma única en el contraste DKO vs WT, pero no en otros dos contrastes KD vs WT, y Ctr vs WT. La Figura 6C valida el cambio significativo distal a proximal del uso del sitio pA para el gen Ccl2 detectado únicamente en el contraste KD vs WT, mientras que la Figura 6D valida el uso significativo de pA diferencial de todos los contrastes para el gen Cacna2d1.

Figura 6. Comparación de contraste y verificación de los resultados de diffSplice. A) Diagrama de Venn que compara los resultados diferenciales significativos de uso del sitio de pA de diferentes contrastes adquiridos de diffSplice. B-D) Instantánea IGV que visualiza picos de cobertura de pA de los genes Paip2, Ccl2 y Cacna2d1 a través de condiciones. Cada pista representa la cobertura de lectura en una condición específica. Haga clic aquí para ver una versión más grande de esta figura.

Figura complementaria 1. Análisis RNA-Seq de empalme diferencial con Limma diffSplice. (A) Diagrama volcánica de RNA-Seq de Limma diffSplice analysis: El eje x muestra el cambio de pliegue del exón logarítmico; El eje Y muestra el valor p -log₁₀ . Cada punto corresponde a un exón. Línea discontinua horizontal en valor p=1E-5; líneas discontinuas verticales con doble cambio (FC). Los exones rojos muestran una FC sustancial y significación estadística. (B-D) Gráficos de empalme diferencial: Los patrones de empalme se exhiben, por ejemplo, los genes Mbnl1, Tcf7 y Lef1, respectivamente, donde el eje x muestra el exón id por transcripción; el eje y muestra el cambio relativo del pliegue logarítmico del exón (la diferencia entre el logFC del exón y el logFC general para todos los demás exones). Los exones resaltados en rojo muestran una expresión diferencial estadísticamente significativa (FDR < 0,1). Haga clic aquí para descargar este archivo.

Figura complementaria 2. Análisis RNA-Seq del uso diferencial de exones con DEXSeq. (A) Parcela de volcán. (B-D) Uso diferencial del exón de RNA-Seq obtenido del análisis DEXSeq. Los genes Mbnl1, Tcf7 y Lef1, respectivamente, muestran un uso diferencial significativo de exones entre WT y Mbnl1 knock-out resaltados en rosa, que corresponden a los mismos exones diferenciales en la figura suplementaria 1. Haga clic aquí para descargar este archivo.  

Figura complementaria 3. Parcelas alternativas de poliadenilación por diffSplice. A) Gráficos volcánicos de datos PolyA-seq generados usando diffSplice en tres pares de contraste obtenidos de los datos PolyA-seq del ratón, incluyendo doble knockout (DKO) vs wild-type (WT), knock-down (KD) vs WT y control (Ctrl) vs WT. El eje X muestra el cambio de pliegue del sitio log pA; El eje Y muestra el valor p -log₁₀. Cada punto corresponde a un sitio pA. Línea discontinua horizontal en valor p=1E-5; líneas discontinuas verticales en 2 veces FC. Los sitios rojos de pA muestran una FC sustancial y significación estadística. B) Gráficos diferenciales APA generados mediante diffSplice que muestran los sitios pA proximal, interno y distal para los genes altamente clasificados S100a7a, Tpm1 y Smc6. Haga clic aquí para descargar este archivo.  

Figura complementaria 4. Análisis diferencial de uso de pA por pipeline DEXSeq. A) Gráficos volcánicos de datos PolyA-seq generados usando DEXSeq en tres pares obtenidos de los datos PolyA-seq del ratón, incluyendo doble knockout (DKO) vs wild-type (WT), knock-down (KD) vs WT, y control (Ctrl) vs WT. El eje X muestra el cambio de pliegue del sitio log pA; El eje Y muestra el valor p -log₁₀. Cada punto corresponde a un sitio pA. Línea discontinua horizontal en valor p=1E-5; líneas discontinuas verticales en 2 veces FC. Los sitios rojos de pA muestran una FC sustancial y significación estadística. B) Gráficos diferenciales APA generados utilizando DEXSeq que muestran los sitios pA proximal, interno y distal para los genes altamente clasificados S100a7a, Tpm1 y Smc6. Haga clic aquí para descargar este archivo.  

Figura complementaria 5. Diagrama de cobertura media y mapas de calor alrededor de los sitios de escisión de pA.  La cobertura se muestra para cuatro afecciones, con genes en las hebras hacia adelante y hacia atrás que se muestran por separado. El eje X muestra la posición de la base en relación con los sitios de escisión de pA, hasta 100 nucleótidos aguas arriba y aguas abajo; El eje y se refiere a la cobertura media en las posiciones de base relativas correspondientes en todos los sitios de escisión de pA. Los mapas de calor proporcionan una vista alternativa, con cada sitio de escisión pA mostrado como una fila separada en el eje x, ordenada por cobertura. La intensidad muestra la cobertura de lectura (ver leyenda). Haga clic aquí para descargar este archivo.

Cuadro complementario 1. Salida diffSplice del análisis AS. La primera pestaña define los nombres de columna para las salidas originales presentadas en la segunda pestaña (nivel de exón) y tercera (nivel de gen). Haga clic aquí para descargar esta tabla.

Cuadro complementario 2. Salida DEXSeq del análisis AS. La primera pestaña define los nombres de columna para la salida original presentada en la segunda pestaña (nivel de exón) y tercera (nivel de gen). Haga clic aquí para descargar esta tabla.

Cuadro complementario 3. Salida rMATS del análisis AS. La primera pestaña define los nombres de columna para el archivo de resumen (pestaña 2) y los archivos JCEC para cada evento (pestaña 3-7). Haga clic aquí para descargar esta tabla.

Cuadro complementario 4. Salida diffSplice del análisis APA. La primera pestaña define los nombres de columna para la salida original presentada en la segunda pestaña (nivel de sitio pA) y tercera (nivel de gen). Haga clic aquí para descargar esta tabla.

Cuadro complementario 5. Salida DEXSeq del análisis APA. La primera pestaña define los nombres de columna para la salida original presentada en la segunda pestaña (nivel de sitio pA) y tercera (nivel de gen). Haga clic aquí para descargar esta tabla.

Cuadro complementario 6. Un resumen del número de exones significativamente cambiados para los sitios AS y pA para APA. Haga clic aquí para descargar esta tabla.

Cuadro complementario 7. Un resumen de las herramientas y paquetes utilizados en el análisis AS/APA. Haga clic aquí para descargar esta tabla.

Discussion

En este estudio, evaluamos enfoques basados en exones y eventos para detectar AS y APA en datos masivos de RNA-Seq y secuenciación final 3'. Los enfoques de EA basados en exones producen tanto una lista de exones expresados diferencialmente como una clasificación a nivel de gen ordenada por la significación estadística de la actividad general de empalme diferencial a nivel de gen (Tablas 1-2, 4-5). Para el paquete diffSplice, el uso diferencial se determina ajustando modelos lineales ponderados a nivel de exón para estimar el cambio de pliegue logarítmico diferencial de un exón contra el cambio promedio de pliegue logarítmico de los otros exones dentro del mismo gen (llamado FC por exón). La significación estadística a nivel genético se calcula combinando pruebas individuales de significación a nivel de exón en una prueba genética mediante el método de Simes¹⁰. Esta clasificación por actividad de empalme diferencial a nivel de gen se puede utilizar posteriormente para realizar un análisis de enriquecimiento de conjuntos de genes (GSEA) de las vías clave involucradas¹⁰. DEXSeq utiliza una estrategia similar, ajustando un modelo lineal generalizado para medir el uso diferencial de exones, aunque difiere en ciertos pasos, como el filtrado, la normalización y la estimación de dispersión. Al comparar los 500 principales exones clasificados que muestran actividad de AS y APA usando DEXSeq y DiffSplice, encontramos una superposición de 310 exones y 300 sitios de pA, respectivamente, lo que demuestra la concordancia de los dos enfoques basados en exones, que también se demostró en un estudio anterior ²⁹. Se recomienda utilizar una combinación de un enfoque basado en exones (DEXSeq o diffSplice) y un enfoque basado en eventos para la detección y clasificación integrales de EA. Para APA, los usuarios pueden elegir DEXSeq o diffSplice: ambos métodos han demostrado funcionar bien en una amplia gama de experimentos transcriptómicos²⁹.

Al preparar la biblioteca RNA-seq para un análisis AS, es importante utilizar un protocolo RNA-seq a granel específico de hebras⁸, ya que una gran fracción de los genes en los genomas de vertebrados se superponen en diferentes hebras, y un protocolo no específico de hebras no puede distinguir estas regiones superpuestas, confundiendo la detección final del exón. Otra consideración es la profundidad de lectura, con análisis de empalme que requieren una secuenciación más profunda que DGE, por ejemplo, 30-60 millones de lecturas por muestra, frente a 5-25 millones de lecturas por muestra para DGE (https://sapac.support.illumina.com/bulletins/2017/04/considerations-for-rna-seq-read-length-and-coverage-.html). Todas las herramientas demostradas en el protocolo admiten datos de secuenciación de extremo único y pareado. Si solo se usan anotaciones genéticas conocidas para detectar lecturas de unión, se pueden usar lecturas más cortas de un solo extremo (≥ 50 pb), aunque la detección de novo de nuevas uniones de empalme se beneficia de lecturas de extremo pareado y más largas (≥ 100 pb)^30,31. La elección del protocolo de extracción de ARN, ya sea la selección de poliA o el agotamiento del ARNr, depende de la calidad del ARN y de la pregunta experimental, consulte³¹ para una discusión. Dependiendo de los detalles de la construcción de la biblioteca, se requerirán modificaciones a los scripts de ejemplo dados aquí para los parámetros de alineación de lectura, recuento de características y rMATS. Al calcular los recuentos iniciales de lectura a nivel de exón utilizando featureCounts, o métodos similares, se debe tener cuidado de configurar las opciones de función correctamente para los recuentos y el strandedness: en featureCounts, establecemos el argumento "strandSpecific" apropiadamente para el protocolo RNA-seq específico de hebra utilizado; y para la cuantificación a nivel de exón, se espera que una lectura se asigne sobre exones adyacentes, por lo que establecemos el parámetro allowMultiOverlap en TRUE. Para APA, existen diferentes protocolos de secuenciación de extremo 3'⁶ que varían en la ubicación precisa de los picos en relación con el sitio de pA. Para nuestros datos de ejemplo, determinamos que el pico es 60 pb aguas arriba del sitio pA como se muestra en la Figura 5, y este análisis deberá adaptarse para otros protocolos de secuenciación final 3'.

En este protocolo limitamos el alcance a la discusión de análisis diferenciales a nivel de exones individuales, y eventos de empalme que consisten en combinaciones adyacentes exón-intrón. No discutimos la clase de análisis basados en la reconstrucción de isoformas de novo como Cufflinks, Cuffdiff³², RSEM 33, Kallisto³⁴ que tienen como objetivo detectar y cuantificar la expresión absoluta y relativa de isoformas alternativas completas. Los métodos basados en exones y eventos son más sensibles para detectar eventos de empalme individuales³⁰ y en muchos casos proporcionan toda la información necesaria para un análisis posterior, sin necesidad de cuantificación a nivel de isoforma.

La última versión de los archivos de origen de este protocolo está disponible en https://github.com/jiayuwen/AS_APA_JoVE

Materials

Name	Company	Catalog Number	Comments
Not relevent for computational study