Identificazione di splicing alternativo e poliadenilazione in dati di RNA-seq

Lo splicing alternativo (AS) e la poliadenilazione alternativa (APA) espandono la diversità delle isoforme di trascrizione e dei loro prodotti. Qui, descriviamo i protocolli bioinformatici per analizzare i saggi di sequenziamento di massa dell'RNA-seq e 3' per rilevare e visualizzare AS e APA che variano in base alle condizioni sperimentali.

Oltre all'analisi tipica dell'RNA-Seq per misurare l'espressione genica differenziale (DGE) in condizioni sperimentali/biologiche, i dati di RNA-seq possono anche essere utilizzati per esplorare altri complessi meccanismi regolatori a livello di esone. Lo splicing alternativo e la poliadenilazione svolgono un ruolo cruciale nella diversità funzionale di un gene generando diverse isoforme per regolare l'espressione genica a livello post-trascrizionale e limitando le analisi all'intero livello genico può mancare questo importante strato regolatore. Qui, dimostriamo analisi dettagliate passo dopo passo per l'identificazione e la visualizzazione dell'utilizzo differenziale del sito di esone e poliadenilazione in tutte le condizioni, utilizzando Bioconductor e altri pacchetti e funzioni, tra cui DEXSeq, diffSplice dal pacchetto Limma e rMATS.

L'RNA-seq è stato ampiamente utilizzato nel corso degli anni, tipicamente per stimare l'espressione genica differenziale e la scoperta^{genica 1}. Inoltre, può anche essere utilizzato per stimare l'uso variabile del livello di esone a causa del gene che esprime diverse isoforme, contribuendo così a una migliore comprensione della regolazione genica a livello post-trascrizionale. La maggior parte dei geni eucariotici genera diverse isoforme mediante splicing alternativo (AS) per aumentare la diversità dell'espressione dell'mRNA. Gli eventi AS possono essere suddivisi in diversi modelli: salto di esoni completi (SE) in cui un esone ("cassetta") viene completamente rimosso dalla trascrizione insieme ai suoi introni laterali; selezione alternativa (donatore) del sito di giunzione 5' (A5SS) e selezione alternativa del sito di giunzione 3' (accettore) (A3SS) quando due o più siti di giunzione sono presenti su entrambe le estremità di un esone; ritenzione degli introni (RI) quando un introne viene trattenuto all'interno del trascritto dell'mRNA maturo e mutua esclusione dell'uso dell'esone (MXE) in cui solo uno dei due esoni disponibili può essere trattenuto alla volta ^2,3. La poliadenilazione alternativa (APA) svolge anche un ruolo importante nella regolazione dell'espressione genica utilizzando siti di poli alternativi (A) per generare più isoforme di mRNA da una singola trascrizione⁴. La maggior parte dei siti di poliadenilazione (pAs) si trovano nella regione 3' non tradotta (3' UTR), generando isoforme di mRNA con diverse lunghezze UTR 3'. Poiché l'UTR 3' è l'hub centrale per il riconoscimento degli elementi regolatori, diverse lunghezze UTR 3' possono influenzare la localizzazione, la stabilità e la traduzione dell'mRNA⁵. Esistono saggi di sequenziamento finale di classe 3' ottimizzati per rilevare APA che differiscono nei dettagli del protocollo⁶. La pipeline qui descritta è progettata per PolyA-seq, ma può essere adattata per altri protocolli come descritto.

In questo studio, presentiamo una pipeline di metodi di analisi differenziale degli esoni^7,8 (Figura 1), che possono essere suddivisi in due grandi categorie: basati sull'esone (DEXSeq^9, diffSplice 10) e basati sugli eventi (replicate Multivariate Analysis of Transcript Splicing (rMATS)¹¹). I metodi basati sull'esone confrontano il cambiamento di piega tra le condizioni dei singoli esoni, contro una misura del cambiamento complessivo della piega genica per chiamare l'uso differenziale dell'esone, e da ciò calcolare una misura a livello genetico dell'attività AS. I metodi basati su eventi utilizzano letture di giunzione esone-introne per rilevare e classificare eventi di splicing specifici come il salto dell'esone o la ritenzione di introni e distinguere questi tipi di AS nell'output³. Pertanto, questi metodi forniscono punti di vista complementari per un'analisi completa di AS^12,13. Abbiamo selezionato DEXSeq (basato sul pacchetto DESeq2¹⁴ DGE) e diffSplice (basato sul pacchetto Limma¹⁰ DGE) per lo studio in quanto sono tra i pacchetti più utilizzati per l'analisi di splicing differenziale. rMATS è stato scelto come metodo popolare per l'analisi basata sugli eventi. Un altro metodo popolare basato su eventi è MISO (Mixture of Isoforms)¹. Per l'APA adattiamo l'approccio basato sull'esone.

Figura 1. Pipeline di analisi. Diagramma di flusso dei passaggi utilizzati nell'analisi. I passaggi includono: ottenere i dati, eseguire controlli di qualità e allineamento delle letture seguito dal conteggio delle letture utilizzando annotazioni per esoni, introni e siti pA noti, filtraggio per rimuovere conteggi bassi e normalizzazione. I dati di PolyA-seq sono stati analizzati per siti pA alternativi utilizzando metodi diffSplice/DEXSeq, RNA-Seq di massa sono stati analizzati per lo splicing alternativo a livello di esone con metodi diffSplice/DEXseq e gli eventi AS analizzati con rMATS. Fare clic qui per visualizzare una versione ingrandita di questa figura.

I dati RNA-seq utilizzati in questa indagine sono stati acquisiti da Gene Expression Omnibus (GEO) (GSE138691)¹⁵. Abbiamo utilizzato i dati RNA-seq di topo di questo studio con due gruppi di condizioni: wild-type (WT) e Muscleblind-like type 1 knockout (Mbnl1 KO) con tre repliche ciascuno. Per dimostrare l'analisi differenziale dell'utilizzo del sito di poliadenilazione, abbiamo ottenuto dati PolyA-seq di fibroblasti embrionali di topo (MEF) (GEO Accession GSE60487)¹⁶. I dati hanno quattro gruppi di condizioni: Wild-type (WT), Muscleblind-like type1/type 2 double knockout (Mbnl1/2 DKO), Mbnl 1/2 DKO con Mbnl3 knockdown (KD) e Mbnl1/2 DKO con controllo Mbnl3 (Ctrl). Ogni gruppo di condizioni è costituito da due repliche.

Tabella 1. Riepilogo dei set di dati RNA-Seq e PolyA-seq utilizzati per l'analisi.

1. Installazione di strumenti e pacchetti R utilizzati nell'analisi

1. Conda è un gestore di pacchetti popolare e flessibile che consente una comoda installazione dei pacchetti con le loro dipendenze su tutte le piattaforme. Utilizzare 'Anaconda' (conda package manager) per installare 'conda' che può essere utilizzato per installare gli strumenti/pacchetti necessari per l'analisi.
2. Scarica 'Anaconda' secondo i requisiti di sistema da https://www.anaconda.com/products/individual#Downloads e installalo seguendo le istruzioni nell'installatore grafico. Installa tutti i pacchetti richiesti usando 'conda' digitando quanto segue nella riga di comando di 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. Per scaricare tutti i pacchetti R utilizzati nel protocollo, digitare il seguente codice nella console R (avviata sulla riga di comando di Linux digitando 'R') o nella console 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: in questo protocollo computazionale, i comandi verranno forniti come file R Notebook (file con estensione ". Rmd"), file di codice R (file con estensione ". R"), o script di shell Linux Bash (file con estensione ".sh"). I file R Notebook (Rmd) devono essere aperti in RStudio utilizzando File| Apri File... e singoli blocchi di codice (che possono essere comandi R o comandi della shell Bash) vengono quindi eseguiti in modo interattivo facendo clic sulla freccia verde in alto a destra. I file di codice R possono essere eseguiti aprendo in RStudio o sulla riga di comando di Linux premettendo "Rscript", ad esempio Rscript example.R. Gli script della shell vengono eseguiti sulla riga di comando di Linux anteponendo lo script con il comando "sh", ad esempio .sh example.sh.

2. Analisi di splicing alternativo (AS) mediante RNA-seq

1. Download e pre-elaborazione dei dati
   NOTA: i frammenti di codice annotati di seguito sono disponibili nel file di codice supplementare "AS_analysis_RNASeq.Rmd", per seguire i singoli passaggi in modo interattivo, e sono forniti anche come script bash da eseguire in batch sulla riga di comando di Linux (sh downloading_data_preprocessing.sh).
   1. Download dei dati grezzi.
      1. Scarica i dati grezzi da Sequence Read Archive (SRA) utilizzando il comando 'prefetch' dal toolkit SRA (v2.10.8)¹⁷. Fornire gli ID di adesione SRA in sequenza nel seguente comando per scaricarli in parallelo utilizzando l'utilità parallela GNU¹⁸. Per scaricare i file SRA degli ID di adesione da SRR10261601 a SRR10261606 in parallelo, utilizzare quanto segue nella riga di comando di Linux.
         ​seq 10261601 10261606 | parallel prefetch SRR{}
      2. Estrarre i file fastq dall'archivio utilizzando la funzione 'fastq-dump' dal toolkit SRA. Usa GNU parallelamente e dai i nomi di tutti i file SRA insieme.
         ​parallel -j 3 fastq-dump --gzip --skip-technical --read-filter pass --dumpbase --split-e --clip --origfmt {} :::
      3. Scaricare il genoma di riferimento e le annotazioni per Mouse (Genome assembly GRCm39) da www.ensembl.org utilizzando quanto segue nella riga di comando di 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. Pre-elaborazione e mappatura delle letture nell'assemblaggio del genoma
      1. Controllo qualità. Valuta la qualità delle letture non elaborate con FASTQC (FASTQ Quality Check v0.11.9)¹⁹. Creare una cartella di output ed eseguire fastqc con parallelo su più file fasta di input. Questo passaggio genererà un rapporto sulla qualità per ogni campione. Esaminare i report per assicurarsi che la qualità delle letture sia accettabile prima di eseguire ulteriori analisi. (Fare riferimento al manuale dell'utente per comprendere i report in https://www.bioinformatics.babraham.ac.uk/projects/fastqc/)
         mkdir fastqc_out
parallel "fastqc {} -o fastqc_out" ::: $RAW_DATA/*.fastq.gz
         NOTA: se necessario, eseguire il taglio dell'adattatore con 'cutadapt'20 o '^{trimmomatic'21} per rimuovere il sequenziamento negli adattatori laterali, che varia in base alla dimensione del frammento di RNA e alla lunghezza di lettura. In questa analisi abbiamo saltato questo passaggio poiché la frazione di letture interessate era minima.
      2. Leggi l'allineamento. Il passaggio successivo nella pre-elaborazione include la mappatura delle letture al genoma di riferimento. In primo luogo, costruire l'indice per il genoma di riferimento utilizzando la funzione 'genomeGenerate' di STAR²² e quindi allineare le letture grezze al riferimento (in alternativa gli indici precostruiti sono disponibili sul sito web STAR e possono essere utilizzati direttamente per l'allineamento). Eseguire i seguenti comandi dalla riga di comando di 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: l'allineatore STAR genererà e ordinerà i file BAM (Binary Alignment Map) per ogni campione dopo l'allineamento di lettura. I file Bam devono essere ordinati prima di procedere a ulteriori passaggi.
2. Preparazione delle annotazioni Exon.
   1. Eseguire il file di codice supplementare "prepare_mm_exon_annotation. R" con l'annotazione scaricata in formato GTF (Gene transfer format) per preparare le annotazioni. Per eseguire, digitare quanto segue nella riga di comando di Linux.
      Rscript prepare_mm_exon_annotation.R annotation.gtf
      NOTA: il file GTF contiene più voci di esone per diverse isoforme. Questo file viene utilizzato per "comprimere" gli ID di trascrizione multipli per ciascun esone. È un passo importante definire i contenitori per il conteggio degli esoni.
3. Letture di conteggio. Il passo successivo è contare il numero di letture mappate su diverse trascrizioni / esoni. Vedere File supplementare: "AS_analysis_RNASeq.Rmd".
   1. Caricare le librerie richieste:
      packages<- c("Rsubread","tidyverse", "magrittr", "EnhancedVolcano", "edgeR","openxlsx")
invisible(lapply(packages, library, character.only=TRUE))
   2. Caricare il file di annotazione elaborato ottenuto dal passaggio precedente (2.2).
      load("mm_exon_anno.RData")
   3. Leggere tutti i file bam ottenuti nel passaggio 2.2.2 come input per 'featureCounts' per contare le letture. Leggere la cartella contenente i file bam elencando prima ogni file dalla directory che termina con .bam. Usa 'featureCounts' dal pacchetto Rsubread che prende i file bam e l'annotazione GTF elaborata (riferimento) come input per generare una matrice di conteggi associati a ciascuna funzionalità con righe che rappresentano esoni (caratteristiche) e colonne che rappresentano campioni.
      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. Successivamente, eseguire il filtraggio non specifico per rimuovere gli esoni poco espressi ("non specifico" indica che le informazioni sulla condizione sperimentale non vengono utilizzate nel filtraggio, per evitare distorsioni di selezione). Trasformare i dati da scala grezza a conteggi per milione (cpm) utilizzando la funzione cpm dal pacchetto 'edgeR'²³ e mantenere gli esoni con conteggi superiori a una soglia impostabile (per questo set di dati viene utilizzato un cpm) in almeno tre campioni. Rimuovere anche i geni con un solo esone.
      # 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: controllare i parametri obbligatori per featureCounts quando si utilizzano dati diversi, ad esempio, per letture a estremità singola, impostare 'isPairedEnd = FALSE'. Fare riferimento alla guida utente di RSubread per scegliere le opzioni per i dati e vedere la sezione Discussione di seguito.
4. Splicing differenziale e analisi dell'utilizzo dell'esone. Descriviamo due alternative per questo passaggio: DEXSeq e DiffSplice. Entrambi possono essere utilizzati e dare risultati simili. Per coerenza, selezionare DEXSeq se si preferisce un pacchetto DESeq2 per DGE e utilizzare DiffSplice per un'analisi DGE basata su Limma. Vedi scheda supplementare: "AS_analysis_RNASeq.Rmd".
   1. Utilizzo del pacchetto DEXSeq per l'analisi differenziale degli esoni.
      1. Caricare la libreria e creare una tabella di esempio per definire il progetto sperimentale.
         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 : i nomi di riga devono essere coerenti con i nomi di file bam utilizzati da featureCounts per contare le letture. sampleTable è costituito dai dettagli di ogni esempio che include: tipo di libreria e condizione. Questo è necessario per definire i contrasti o il gruppo di test per rilevare l'utilizzo differenziale.
      2. Preparare il file di informazioni sull'esone. Le informazioni sugli esoni sotto forma di oggetti GRanges (genomic ranges) (https://bioconductor.org/packages/release/bioc/html/GenomicRanges.html) sono necessarie come input per creare l'oggetto DEXSeq nel passaggio successivo. Abbina gli ID genetici con i conteggi di lettura per creare l'oggetto 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. Creare un oggetto DEXSeq utilizzando la funzione DEXSeqDataSet. L'oggetto DEXSeq raccoglie insieme i conteggi di lettura, le informazioni sulle feature dell'esone e le informazioni di esempio. Utilizzare i conteggi di lettura generati nel passaggio 3 e le informazioni sull'esone ottenute dal passaggio precedente per creare l'oggetto DEXSeq dalla matrice di conteggio. L'argomento sampleData accetta un input di frame di dati che definisce i campioni (e i relativi attributi: tipo di libreria e condizione), 'design' utilizza sampleData per generare una matrice di progettazione per il test differenziale utilizzando la notazione della formula del modello. Si noti che un termine di interazione significativo, condizione:esone, indica che la frazione di letture su un gene che cade su un particolare esone dipende dalla condizione sperimentale, cioè c'è AS. Vedere la documentazione di DEXSeq per una descrizione completa dell'impostazione della formula del modello per progetti sperimentali più complessi. Per le informazioni sulle caratteristiche, sono necessari gli ID dell'esone, il gene corrispondente e le trascrizioni.
         ​dxd<-DEXSeqDataSet(countData$counts,sampleData=sampleTable, design=~sample+exon+condition:exon,featureID=exoninfo$ExonID,groupID=exoninfo$GeneID,featureRanges=exoninfo, transcripts=transcripts_l)
      4. Normalizzazione e stima della dispersione. Successivamente, eseguire la normalizzazione tra i campioni e stimare la varianza dei dati, dovuta sia al rumore del conteggio di Poisson dalla natura discreta dell'RNA-seq che alla variabilità biologica, utilizzando i seguenti comandi.
         ​dxd %<>% estimateSizeFactors %>% estimateDispersions %T>% plotDispEsts
      5. Test per l'utilizzo differenziale. Dopo la stima della variazione, testare l'uso differenziale dell'esone per ciascun gene e generare i risultati.
         dxd%<>%testForDEU%>%estimateExonFoldChanges(fitExpToVar=
"condition")#Estimate fold changes
dxr=DEXSeqResults(dxd)
      6. Visualizza gli eventi di splicing per i geni selezionati usando il seguente comando.
         plotDEXSeq(dxr,"Wnk1", displayTranscripts=TRUE, splicing=TRUE,legend
=TRUE,cex.axis=1.2,cex=1.3,lwd=2)
         Esaminare il file R Notebook "AS_analysis_RNASeq.Rmd" per generare grafici aggiuntivi per i geni di interesse e per generare grafici vulcanici a soglie diverse.
   2. Utilizzo di diffSplice di Limma per identificare lo splicing differenziale. Segui il file R Notebook "AS_analysis_RNASeq.Rmd". Assicurarsi che siano stati seguiti i passaggi 2.1-2.3 per preparare i file di input prima di procedere ulteriormente.
      1. Caricare le librerie
         library(limma)
library(edgeR)
      2. Filtraggio non specifico. Estrarre la matrice dei conteggi di lettura ottenuti in 2.3. Crea un elenco di funzionalità usando la funzione 'DGEList' dal pacchetto edgeR, dove le righe rappresentano i geni e le colonne rappresentano i campioni.
         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: come passaggio di filtraggio non specifico, i conteggi vengono filtrati in base a cpm < 1 in x su n campioni, dove x è il numero minimo di repliche in qualsiasi condizione. n = 6 e x = 3 per questi dati di esempio.
      3. Normalizzare i conteggi tra i campioni, con la funzione 'calcNormFactors' dal pacchetto 'edgeR' utilizzando i valori Trimmed Mean of M (metodo di normalizzazione TMM)²⁴ Calcolerà i fattori di scala per regolare le dimensioni della libreria.
         ​dge<-calcNormFactors(dge)
      4. Utilizzare sampleTable come generato nel passaggio 2.4.1.1 e creare la matrice di progettazione. La matrice di progettazione caratterizza il design. Vedere la Guida per l'utente di Limma (https://www.bioconductor.org/packages/devel/bioc/vignettes/limma/inst/doc/usersguide.pdf) capitoli 8 e 9 per i dettagli sulle matrici di progettazione per progetti sperimentali più avanzati.
         Treat<- factor(sampleTable$condition)
design<- model.matrix(~0+Treat)
colnames(design) <- levels(Treat)
      5. Montare un modello lineare per esone. Eseguire la funzione 'voom' del pacchetto 'limma' per elaborare i dati RNA-seq per stimare la varianza e generare pesi di precisione per correggere il rumore del conteggio di Poisson e trasformare i conteggi a livello di esone in log2-counts per milione (logCPM). Quindi eseguire la modellazione lineare utilizzando la funzione 'lmfit' per adattare i modelli lineari ai dati di espressione per ciascun esone. Calcola le statistiche empiriche di Bayes per il modello adattato utilizzando la funzione 'eBayes' per rilevare l'espressione differenziale dell'esone. Quindi, definire una matrice di contrasto per i confronti sperimentali di interesse. Usa 'contrasts.fit' per ottenere coefficienti ed errori standard per ogni coppia di confronti.
         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. Analisi di splicing differenziale. Esegui 'diffSplice' sul modello adattato per testare le differenze nell'uso degli esoni dei geni tra wild-type e knockout ed esplorare i risultati migliori usando la funzione 'topSplice': test="t" fornisce una classificazione degli esoni AS, test="simes" fornisce una classificazione dei geni.
         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. Visualizzazione. Tracciare i risultati con la funzione 'plotSplice', dando il gene di interesse nell'argomento geneide. Salva i risultati principali ordinati per log Piega la modifica in un oggetto e genera un grafico vulcanico per mostrare gli esoni.
         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. Utilizzo di rMATS
      1. Assicurarsi che l'ultima versione di rMATS v4.1.1 (nota anche come rMATS turbo a causa del tempo di elaborazione ridotto e dei minori requisiti di memoria) sia installata utilizzando conda o github (https://github.com/Xinglab/rmats-turbo/releases/download/v4.1.1/rmats_turbo_v4_1_1.tar.gz) nella directory di lavoro. Seguire il paragrafo 4.3 in "AS_analysis_RNASeq.Rmd".
      2. Passare alla cartella contenente i file bam ottenuti dopo il mapping e preparare i file di testo, come richiesto da rMATS, per le due condizioni copiando il nome dei file bam (insieme al percorso) separati da ','. I seguenti comandi devono essere eseguiti dalla riga di comando di 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. Eseguire rmats.py con i due file di input generati nel passaggio precedente, insieme al file GTF ottenuto nella versione 2.1.1.3. Questo genererà una cartella di output 'rmats_out' contenente file di testo che descrivono le statistiche (valori p e livelli di inclusione) per ogni evento di giunzione separatamente.
         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 : è richiesta anche l'annotazione di riferimento sotto forma di file GTF. Controllare i parametri se i dati sono single-end e modificare l'opzione -t di conseguenza.
      4. Esplorazione dei risultati di rMATS. Utilizzare il pacchetto Bioconductor 'maser'25 per esplorare i risultati di rMATS. Carica i file di testo JCEC (Junction and exon counts) nell'oggetto 'maser' e filtra il risultato in base alla copertura includendo almeno cinque letture medie per evento di giunzione.
         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. Visualizzazione dei risultati rMATS. Selezionare gli eventi di splicing significativi al 10% del tasso di falsa scoperta (FDR) e la variazione minima del 10% della percentuale di giunzione in (deltaPSI) utilizzando la funzione 'topEvents' dal pacchetto 'maser'. Quindi, controllare gli eventi genetici per i singoli geni di interesse (gene campione-Wnk1) e tracciare i valori PSI per ogni evento di splicing di quel gene. Generare un grafico a vulcano specificando il tipo di 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. Genera grafici Sashimi per il risultato degli eventi di giunzione ottenuti con rMATS sotto forma di file di testo utilizzando il pacchetto 'rmats2shahimiplot'. Eseguire lo script python dalla riga di comando di 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: questo processo può richiedere molto tempo in quanto genererà il grafico Sashimi per tutti i risultati nel file degli eventi. Scegli i risultati migliori (nomi di geni ed esoni) come visualizzati dalla funzione topEvents da 'maser' e visualizza il grafico del Sashimi corrispondente.

3. Analisi di poliadenilazione alternativa (APA) mediante sequenziamento finale 3'

1. Download e pre-elaborazione dei dati
   Nota : fare riferimento al file supplementare R Notebook "APA_analysis_3PSeq_notebook. Rmd" per i comandi completi per il download dei dati e le fasi di pre-elaborazione, oppure eseguire il file bash supplementare "APA_data_downloading_preprocessing.sh" sulla riga di comando di Linux.
   1. Scaricare i dati dall'SRA con gli ID di adesione (da 1553129 a 1553136).
   2. Adattatori trim e complemento inverso per ottenere la sequenza del filo di senso.
      NOTA: questo passaggio è specifico per il test PolyA-seq utilizzato.
   3. La mappa legge l'assemblaggio del genoma del mouse utilizzando l'allineatore per papillon²⁶.
2. Preparazione delle annotazioni dei siti pA.
   NOTA: l'elaborazione del file di annotazione del sito pA viene eseguita in primo luogo utilizzando il file supplementare R Notebook "APA_analysis_3PSeq_notebook. Rmd" (2.1 - 2.6), quindi utilizzando il file bash "APA_annotation_preparation.sh".
   1. Scarica l'annotazione dei siti pA dal database PolyASite 2.0⁶.
   2. Selezionare le annotazioni del sito pA per conservare i siti pA della regione non tradotta (UTR) 3', annotati come esone terminale (TE) o 1000 nt a valle di un esone terminale annotato (DS) per l'analisi a valle.
   3. Ottieni picchi del sito pA. Ancoraggio in ogni sito di scissione pA e visualizzare la copertura media di lettura utilizzando bedtools e deeptools^27,28. I risultati hanno mostrato che i picchi delle letture mappate erano principalmente dispersi entro ~ 60 bp a monte dei siti di scissione (figura 5 e figura supplementare 5). Pertanto, le coordinate dei siti pA sono state estese dal file di annotazione a 60 bp a monte dei loro siti di scissione. A seconda del protocollo di sequenziamento finale 3' specifico, questo passaggio dovrà essere ottimizzato per saggi diversi da PolyA-seq.
3. Letture di conteggio
   1. Preparare il file di annotazione dei siti 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. Applica 'featureCounts' per acquisire conteggi non elaborati. Salvare la tabella di conteggio come file RData "APA_countData.Rdata" per l'analisi APA utilizzando strumenti diversi.
      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: tenere presente di modificare uno qualsiasi dei parametri elencati nella funzione 'featureCounts'. Modificare il parametro 'strandSpecific' per assicurarsi che sia coerente con la direzione di sequenziamento del saggio di sequenziamento finale 3' utilizzato (empiricamente, visualizzare i dati in un browser del genoma sui geni sui filamenti più e meno chiarirà questo).
   3. Applicare un filtro non specifico di countData. Il filtraggio può migliorare significativamente la robustezza statistica nei test differenziali di utilizzo del sito pA. In primo luogo, abbiamo rimosso quei geni con un solo sito pA, su cui l'utilizzo differenziale del sito pA non può essere definito. In secondo luogo, applichiamo un filtraggio non specifico basato sulla copertura: i conteggi sono filtrati da cpm meno di 1 su x su n campioni, dove x è il numero minimo di repliche in qualsiasi condizione. N = 8 e x = 2 per questi dati di esempio.
      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. Analisi differenziale dell'utilizzo del sito di poliadenilazione mediante pipeline DEXSeq e diffSplice.
   1. Utilizzo del pacchetto DEXSeq
      NOTA: Poiché non è possibile definire una matrice di contrasto per la pipeline DEXSeq, l'analisi APA differenziale di ciascuna due condizioni sperimentali deve essere eseguita separatamente. L'analisi APA differenziale della condizione WT e della condizione DKO viene eseguita come esempio per spiegare la procedura. Fare riferimento al file supplementare "APA_analysis_3PSeq_notebook. RMD ·" per il flusso di lavoro passo-passo di questa sezione e l'analisi APA differenziale di altri contrasti.
      1. Caricare la libreria e creare una tabella di esempio per definire il progetto sperimentale.
         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. Preparare il file di informazioni sui siti pA utilizzando il pacchetto 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. Utilizzare i conteggi di lettura generati nel passaggio 3.3 e le informazioni sul sito pA ottenute dal passaggio precedente per creare l'oggetto 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. Definire la coppia di contrasto definendo i livelli delle condizioni nell'oggetto DEXSeq.
         dxd1$condition<- factor(dxd1$condition, levels= c("WT", "DKO"))
# The contrast pair will be "DKO - WT"
      5. Normalizzazione e stima della dispersione. Analogamente ai dati RNA-seq, per 3' i dati di sequenziamento finale eseguire la normalizzazione tra i campioni (mediana dei rapporti per colonna per ciascun campione) utilizzando la funzione 'estimateSizeFactors' e stimare la variazione dei dati utilizzando la funzione 'estimateDispersions', quindi visualizzare il risultato della stima della dispersione utilizzando la funzione 'plotDispEsts'.
         ​dxd1 %<>% estimateSizeFactors %>% estimateDispersions %T>% plotDispEsts
      6. Test differenziale di utilizzo del sito pA per ciascun gene utilizzando la funzione 'testForDEU', quindi stimare il cambiamento di piega dell'utilizzo del sito pA utilizzando la funzione 'estimateExonFoldChanges'. Controllare i risultati utilizzando la funzione 'DEXSeqResults' e impostare 'FDR < 10%' come criterio per siti pA significativamente differenziali.
         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. Visualizzazione dei risultati dell'utilizzo differenziale del sito pA utilizzando grafici APA differenziali generati dalla funzione 'plotDEXSeq' e diagramma vulcano dalla funzione '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. Utilizzo del pacchetto diffSplice. Fare riferimento al file supplementare R Notebook "APA_analysis_3PSeq_notebook. Rmd" per il flusso di lavoro dettagliato di questa sezione.
      1. Definire i contrasti di interesse per l'analisi dell'utilizzo differenziale di pA.
         Nota : questo passaggio deve essere eseguito dopo la costruzione e l'elaborazione dell'oggetto DGEList, incluso nel file 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. Visualizza il risultato dei contrasti di interesse (qui "DKO - WT") utilizzando grafici APA differenziali con la funzione 'plotSplice' e grafici vulcanologici con la funzione 'EnhancedVolcano'. Fare riferimento al file R Notebook "APA_analysis_3PSeq_notebook. Rmd" 4.2.7 - 4.2.9 per la visualizzazione di altre coppie di contrasto.
         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")

Dopo aver eseguito il flusso di lavoro dettagliato di cui sopra, i risultati dell'analisi AS e APA e i risultati rappresentativi sono sotto forma di tabelle e grafici di dati, generati come segue.

COME:
Il principale risultato dell'analisi AS (Tabella supplementare 1 per diffSplice; La tabella 2 per DEXSeq) è un elenco di esoni che mostrano un uso differenziale tra le condizioni e un elenco di geni che mostrano una significativa attività complessiva di splicing di uno o più dei suoi esoni costituenti, classificati in base alla significatività statistica. La tabella supplementare 1, scheda 2 mostra esoni significativi, con colonne che mostrano FC differenziale di esone rispetto a riposo, valore p non aggiustato per esone e valore p aggiustato (correzione di Benjamini-Hockberg). La soglia sui valori p aggiustati darà un insieme di esoni con FDR definito. La tabella supplementare 1, scheda 3 mostra un elenco classificato di geni che mostrano il significato dell'attività complessiva di splicing, con una colonna che mostra il valore p aggiustato a livello di gene calcolato usando il metodo Simes. Dati simili sono riportati nella tabella 2 per DEXSeq. La Figura supplementare 1 e la Figura supplementare 2 mostrano il modello di splicing differenziale nei geni Mbnl1, Tcf7 e Lef1 che sono stati convalidati sperimentalmente nell'articolo pubblicato presentato con i dati¹⁵. Gli autori hanno mostrato la convalida sperimentale di cinque geni: Mbnl1, Mbnl2, Lef1, Tcf7 e Ncor2. Il nostro approccio ha rilevato pattini di splicing differenziale in tutti questi geni. Qui presentiamo i livelli di FDR per ciascun gene utilizzando rispettivamente DEXSeq, diffSplice e rMATS ottenuti nelle tabelle supplementari 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) e Ncor2 (9.2E-11, 1.6E-06, 0).

La Figura 2 mostra un confronto tra gli output ottenuti da tre diversi strumenti e illustra modelli di splicing alternativi nel gene Wnk1 . I grafici del vulcano sono mostrati in Figura 2A (diffSplice) e Figura 2B (DEXSeq). Altri tre geni altamente classificati sono mostrati nella Figura supplementare 1 (diffSplice) e nella Figura supplementare 2 (DEXSeq). L'asse y mostra la significatività statistica (-log10 P-valori) e l'asse x mostra la dimensione dell'effetto (cambio di piega). I geni situati nei quadranti in alto a sinistra o a destra indicano una FC sostanziale e una forte evidenza statistica di differenze autentiche.

Figura 2. Confronto dei risultati di splicing alternativo ottenuti da diffSplice, DEXSeq e rMATS. |
(A) Grafico del vulcano (a sinistra) di RNA-Seq dall'analisi di Limma diffSplice: L'asse x mostra il cambiamento della piega dell'esone logaritmico; L'asse y mostra il valore p -log₁₀ . Ogni punto corrisponde ad un esone. Linea tratteggiata orizzontale al valore p=1E-5; linee tratteggiate verticali a due cambi di piega (FC). Gli esoni rossi mostrano una sostanziale FC e significatività statistica. Grafico di splicing differenziale (a destra): i modelli di splicing sono esposti per un esempio di gene Wnk1 in cui l'asse x mostra l'id dell'esone per trascrizione; l'asse y mostra la variazione della piega logaritmica relativa dell'esone (la differenza tra logFC dell'esone e logFC complessiva per tutti gli altri esoni). Gli esoni evidenziati in rosso mostrano un'espressione differenziale statisticamente significativa (FDR < 0,1).
(B) Grafico del vulcano (a sinistra) e utilizzo differenziale dell'esone (a destra) dell'RNA-Seq ottenuto dall'analisi DEXSeq. Il gene Wnk1 mostra un uso differenziale significativo degli esoni tra WT e Mbnl1 knock-out evidenziato in rosa, che corrispondono agli stessi esoni differenziali in (A).
(C) Grafico del vulcano (a sinistra) e grafico del sashimi (a destra) per Wnk1 ottenuto dall'analisi rMATS. Grafico del vulcano che raffigura un significativo evento di esone (SE) saltato (cassetta) in wild-type rispetto al knockout al 10% FDR con variazione dei valori percentuali di giunzione (PSI o ΔΨ) > 0,1. L'asse x mostra la variazione dei valori PSI tra le condizioni, mentre l'asse y mostra il valore P del registro. Il grafico del sashimi mostra un evento di esone saltato nel gene Wnk1 , corrispondente a un esone differenziale significativo in (A) e (B). Ogni riga rappresenta un campione di RNA-Seq: tre repliche di wild-type e Mbnl1 knock-out. L'altezza mostra la copertura di lettura in RPKM e gli archi di collegamento raffigurano le letture di giunzione attraverso gli esoni. Le isoforme alternative annotate del modello genetico sono mostrate nella parte inferiore del grafico. Il pannello inferiore di C illustra le letture di giunzione utilizzate per calcolare la statistica PSI.
(D) Diagramma di Venn che confronta il numero di esoni differenziali significativi ottenuti con i diversi metodi. Fare clic qui per visualizzare una versione ingrandita di questa figura.

Figura 2 A (pannello di destra) mostra una visualizzazione schematica delle differenze di esone di uno dei geni classificati in alto, mostrando logFC sull'asse y e il numero di esoni sull'asse x. Questo esempio mostra un esone a cassetta che varia tra le condizioni per il gene Wnk1. Il grafico di utilizzo differenziale dell'esone di DEXSeq mostra prove di splicing differenziale in cinque siti di esoni vicino a Wnk1.6.45. È probabile che gli esoni evidenziati in rosa vengano giuntati nei campioni di Mbnl1 KO rispetto a WT. Questi esoni sono complementari ai risultati ottenuti da diffSplice che mostra un pattern simile nella posizione genomica specifica. Altri esempi sono mostrati nella Figura supplementare 1 e nella Figura supplementare 2. Una visione più dettagliata per confermare risultati interessanti può essere fornita confrontando le tracce di copertura (wiggle) in unità RPM (letture per milione) nei browser genomici UCSC (Università di Santa Cruz) o IGV (Integrative Genomics Viewer) (non mostrato), insieme alla correlazione visiva con altre tracce di interesse, come modelli genici noti, conservazione e altri saggi sull'intero genoma.

La tabella di output rMATS elenca gli eventi di splicing alternativi significativi classificati per tipo (Tabella supplementare 3). La figura 2C mostra un grafico vulcanico di geni che sono giuntati alternativamente, con la dimensione dell'effetto misurata dalla statistica differenziale "percentuale spliced in" (PSI o ΔΨ) di¹¹.

PSI si riferisce alla percentuale di letture coerenti con l'inclusione di un esone a cassetta (cioè, legge la mappatura all'esone della cassetta stesso o le letture di giunzione che si sovrappongono all'esone) rispetto alle letture coerenti con l'esclusione dell'esone, cioè le letture di giunzione attraverso esoni adiacenti a monte e a valle (Il pannello inferiore della Figura 2C). Il pannello di destra della Figura 2C mostra il diagramma del sashimi del gene Wnk1 con evento di splicing differenziale sovrapposto alle tracce di copertura per il gene, con un esone saltato in Mbnl1 KO. Gli archi che uniscono gli esoni mostrano il numero di letture di giunzione (letture che attraversano un introne giuntato). Diverse schede della Tabella supplementare 3 mostrano letture significative di ciascun tipo di evento che si estende su giunzioni con confini di esoni (conteggi di giunzioni e conteggi di esoni (JCEC)). La Figura 2D confronta gli esoni con splicing differenziale significativi rilevati dai tre strumenti.

Figura 3. Eventi di splicing alternativi acquisiti mediante analisi rMATS. A) Tipi di eventi AS. Questa figura è adattata dalla documentazione rMATS¹¹ che spiega l'evento di splicing con esoni costitutivi e alternativamente spliced. Etichettato con il numero di ogni evento a FDR 10%. B) Diagramma di sashimi del gene Add3 che mostra esone saltato (SE). C) Sashimi plot del gene Baz2b che mostra un sito di giunzione alternativo 5' (A5SS). D) Diagramma di sashimi del gene Lsm14b che mostra un sito di giunzione 3' alternativo (A3SS). E) Sashimi plot del gene Mta1 che esibisce esoni mutuamente esclusivi (MXE). F) Sashimi plot del gene Arpp21 che mostra introne trattenuto (RI). Le righe rosse rappresentano tre repliche di wild-type e le righe arancioni rappresentano le repliche knock-out Mbnl1. L'asse x corrisponde alle coordinate genomiche e alle informazioni sul filamento, l'asse y mostra la copertura in RPKM. Fare clic qui per visualizzare una versione ingrandita di questa figura.

La Figura 3 illustra i tipi di eventi di splicing SE, A5SS, A3SS, MXE e RI con l'aiuto dei grafici di Sashimi dei principali geni significativi di tali eventi. Confrontando le tre repliche di WT e Mbnl1_KO, un totale di 1272 eventi SE, 130 A5SS, 116 A3SS, 215 MXE e 313 eventi RI sono stati rilevati a FDR 10%. Il grafico del sashimi illustra il tipo di evento usando i geni principali come esempio.

APA:
L'output dell'analisi APA è simile all'analisi AS a livello di esone. Viene fornita una tabella dei geni principali classificati in base all'attività differenziale dell'APA nel 3'UTR (Tabella supplementare 4 e Tabella supplementare 5). La Figura 4A mostra i grafici vulcanici dei geni per attività APA differenziale in 3'UTR generati utilizzando diffSplice e DEXSeq separatamente. La Figura 4B mostra il grafico di Venn confrontando i risultati di utilizzo del sito pA significativamente differenziali acquisiti da diverse pipeline. Le Figure 4C e 4D mostrano la rappresentazione schematica dell'utilizzo differenziale del sito pA nel 3'UTR del gene Fosl2 (Figura 4C) e Papola (Figura 4D) generato utilizzando sia diffSplice che DEXSeq, che sono convalidati sperimentalmente per mostrare un significativo spostamento da distale a prossimale (Fosl2) e un significativo spostamento da prossimale a distale (Papola) dell'utilizzo del sito pA in DKO, rispettivamente. Altri esempi sono mostrati nella Figura supplementare 3 e nella Figura supplementare 4.

Figura 4. Grafici di poliadenilazione alternativi mediante diffSplice e DEXSeq. A) Grafici vulcanologici di dati PolyA-seq generati utilizzando diffSplice e DEXSeq. L'asse X mostra la modifica della piega del sito pA del registro; L'asse y mostra il valore p -log₁₀ . Ogni punto corrisponde a un sito pA. Linea tratteggiata orizzontale al valore p=1E-5; linee tratteggiate verticali a 2 volte FC. Gli esoni rossi mostrano una sostanziale FC e significatività statistica. B) Grafico di Venn che confronta i risultati di utilizzo del sito pA differenziale significativo acquisiti da diverse pipeline. C-D) Grafici APA differenziali generati utilizzando diffSplice e DEXSeq che mostrano i siti pA prossimale, interno e distale per i geni Fosl2 e Papola . Le figure sono generate dalla stessa funzione della Figura 2 (B) ma con siti pA che sostituiscono gli esoni. Fare clic qui per visualizzare una versione ingrandita di questa figura.

La Figura 5 è un grafico diagnostico per confermare la distribuzione di lettura prevista attorno ai siti di scissione pA annotati per il test PolyA-seq utilizzato. Mostra la copertura media nelle regioni fiancheggianti ancorate a siti di scissione pA noti a livello di genoma. In questo caso, viene visualizzata l'accumulo previsto di letture a monte dei siti. Le distribuzioni di lettura ancorate ai siti pA per tutti i campioni PolyA-seq sono mostrate nella Figura supplementare 5.

Figura 5. Grafico di copertura media intorno ai siti di scissione pA. Il sito di scissione per i dati PolyA-seq è mostrato da una linea tratteggiata verticale. L'asse X mostra la posizione della base rispetto ai siti di scissione pA, fino a 100 nucleotidi a monte e a valle; L'asse y mostra la copertura media di lettura su tutti i siti di scissione pA, normalizzata in base alle dimensioni della libreria in CPM. Fare clic qui per visualizzare una versione ingrandita di questa figura.

I risultati APA differenziali di diversi contrasti generati dalla stessa pipeline possono essere confrontati e verificati visualizzando la copertura di lettura di siti pA rappresentativi significativamente differenziali nel browser del genoma. La Figura 6A è il grafico di Venn che confronta l'utilizzo significativamente differenziale del sito pA di diversi contrasti acquisiti da diffSplice. La Figura 6B-D sono le istantanee IGV della copertura letta nei siti pA per diversi geni, che mostrano i modelli coerenti con quelli scoperti nell'analisi APA utilizzando diffSplice. La Figura 6B convalida il significativo spostamento da prossimale a distale dell'utilizzo del sito pA per il gene Paip2, che è rilevato in modo univoco nel contrasto DKO vs WT, ma non in altri due contrasti KD vs WT e Ctr vs WT. La Figura 6C convalida il significativo spostamento da distale a prossimale dell'utilizzo del sito pA per il gene Ccl2 rilevato in modo univoco nel contrasto KD vs WT, mentre la Figura 6D convalida l'uso differenziale significativo di pA di tutti i contrasti per il gene Cacna2d1.

Figura 6. Confronto del contrasto e verifica dei risultati di diffSplice. A) Diagramma di Venn che confronta i risultati significativi dell'utilizzo differenziale del sito pA di diversi contrasti acquisiti da diffSplice. B-D) Istantanea IGV che visualizza la copertura dei picchi pA dei geni Paip2, Ccl2 e Cacna2d1 in tutte le condizioni. Ogni traccia rappresenta la copertura di lettura in una condizione specifica. Fare clic qui per visualizzare una versione ingrandita di questa figura.

Figura supplementare 1. Analisi RNA-Seq dello splicing differenziale con Limma diffSplice. (A) Diagramma vulcanologico di RNA-Seq dall'analisi di Limma diffSplice: l'asse x mostra il cambiamento della piega dell'esone logaritmico; L'asse y mostra il valore p -log₁₀ . Ogni punto corrisponde ad un esone. Linea tratteggiata orizzontale al valore p=1E-5; linee tratteggiate verticali a doppio cambiamento (FC). Gli esoni rossi mostrano una sostanziale FC e significatività statistica. (B-D) Grafici di giunzione differenziale: I modelli di splicing sono esposti ad esempio i geni Mbnl1, Tcf7 e Lef1, rispettivamente, dove l'asse x mostra l'id dell'esone per trascritto; l'asse y mostra la variazione della piega logaritmica relativa dell'esone (la differenza tra logFC dell'esone e logFC complessiva per tutti gli altri esoni). Gli esoni evidenziati in rosso mostrano un'espressione differenziale statisticamente significativa (FDR < 0,1). Clicca qui per scaricare questo file.

Figura supplementare 2. Analisi RNA-Seq dell'utilizzo differenziale dell'esone con DEXSeq. (A) Trama del vulcano. (B-D) Uso differenziale dell'esone di RNA-Seq ottenuto dall'analisi DEXSeq. I geni Mbnl1, Tcf7 e Lef1, rispettivamente, mostrano un uso differenziale significativo degli esoni tra WT e Mbnl1 knock-out evidenziato in rosa, che corrispondono agli stessi esoni differenziali nella figura supplementare 1. Clicca qui per scaricare questo file.  

Figura supplementare 3. Grafici alternativi di poliadenilazione mediante diffSplice. A) Grafici vulcanologici di dati PolyA-seq generati utilizzando diffSplice in tre coppie di contrasto ottenuti dai dati PolyA-seq del mouse, tra cui double knockout (DKO) vs wild-type (WT), knock-down (KD) vs WT e control (Ctrl) vs WT. L'asse X mostra il cambiamento di piegatura del sito pA del log; L'asse y mostra il valore p -log₁₀. Ogni punto corrisponde a un sito pA. Linea tratteggiata orizzontale al valore p=1E-5; linee tratteggiate verticali a 2 volte FC. I siti pA rossi mostrano una sostanziale FC e significatività statistica. B) Grafici APA differenziali generati utilizzando diffSplice che mostrano i siti pA prossimali, interni e distali per i geni altamente classificati S100a7a, Tpm1 e Smc6. Clicca qui per scaricare questo file.  

Figura supplementare 4. Analisi dell'utilizzo differenziale di pA tramite pipeline DEXSeq. A) Grafici vulcano di dati PolyA-seq generati utilizzando DEXSeq in tre coppie ottenute dai dati PolyA-seq del mouse, tra cui doppio knockout (DKO) vs wild-type (WT), knock-down (KD) vs WT e controllo (Ctrl) vs WT. L'asse X mostra il cambiamento di piega del sito pA del registro; L'asse y mostra il valore p -log₁₀. Ogni punto corrisponde a un sito pA. Linea tratteggiata orizzontale al valore p=1E-5; linee tratteggiate verticali a 2 volte FC. I siti pA rossi mostrano una sostanziale FC e significatività statistica. B) Grafici APA differenziali generati utilizzando DEXSeq che mostrano i siti pA prossimali, interni e distali per i geni altamente classificati S100a7a, Tpm1 e Smc6. Clicca qui per scaricare questo file.  

Figura supplementare 5. Trama di copertura media e mappe di calore intorno ai siti di scissione pA.  Copertura mostrata per quattro condizioni, con geni su filamenti avanti e indietro mostrati separatamente. L'asse X mostra la posizione della base rispetto ai siti di scissione pA, fino a 100 nucleotidi a monte e a valle; L'asse y si riferisce alla copertura media nelle corrispondenti posizioni di base relative in tutti i siti di scissione pA. Le mappe di calore forniscono una vista alternativa, con ogni sito di scissione pA mostrato come una riga separata sull'asse x, ordinata per copertura. L'intensità mostra la copertura di lettura (vedi legenda). Clicca qui per scaricare questo file.

Tabella supplementare 1. diffSplice output dell'analisi AS. La prima scheda definisce i nomi delle colonne per gli output originali presentati nella seconda (livello esone) e nella terza (livello genetico). Clicca qui per scaricare questa tabella.

Tabella supplementare 2. Output DEXSeq dell'analisi AS. La prima scheda definisce i nomi delle colonne per l'output originale presentato nella seconda (livello esone) e nella terza (livello gene). Clicca qui per scaricare questa tabella.

Tabella supplementare 3. Output rMATS dell'analisi AS. La prima scheda definisce i nomi delle colonne per il file di riepilogo (scheda 2) e i file JCEC per ciascun evento (scheda 3-7). Clicca qui per scaricare questa tabella.

Tabella supplementare 4. diffSplice output dell'analisi APA. La prima scheda definisce i nomi delle colonne per l'output originale presentato nella seconda (a livello di sito pA) e nella terza (a livello di gene). Clicca qui per scaricare questa tabella.

Tabella supplementare 5. Output DEXSeq dell'analisi APA. La prima scheda definisce i nomi delle colonne per l'output originale presentato nella seconda (a livello di sito pA) e nella terza (a livello di gene). Clicca qui per scaricare questa tabella.

Tabella supplementare 6. Un riepilogo del numero di esoni significativamente modificati per i siti AS e pA per APA. Clicca qui per scaricare questa tabella.

Tabella supplementare 7. Riepilogo degli strumenti e dei pacchetti utilizzati nell'analisi AS/APA. Clicca qui per scaricare questa tabella.

In questo studio, abbiamo valutato approcci basati su esoni e basati su eventi per rilevare AS e APA in massa RNA-Seq e 3' dati di sequenziamento finale. Gli approcci AS basati sugli esoni producono sia un elenco di esoni differenzialmente espressi sia una classificazione a livello genico ordinata in base alla significatività statistica dell'attività complessiva di splicing differenziale a livello genetico (Tabelle 1-2, 4-5). Per il pacchetto diffSplice, l'uso differenziale è determinato adattando modelli lineari ponderati a livello di esone per stimare il cambiamento di piega del log differenziale di un esone rispetto al cambiamento medio di piega logaritmico degli altri esoni all'interno dello stesso gene (chiamato per esone FC). La significatività statistica a livello genetico viene calcolata combinando i singoli test di significatività a livello di esone in un test genetico con il metodo Simes¹⁰. Questa classificazione in base all'attività di splicing differenziale a livello genico può essere successivamente utilizzata per eseguire un'analisi di arricchimento del set genico (GSEA) dei percorsi chiave coinvolti¹⁰. DEXSeq utilizza una strategia simile, adattando un modello lineare generalizzato per misurare l'utilizzo differenziale dell'esone, sebbene differisca in alcuni passaggi come il filtraggio, la normalizzazione e la stima della dispersione. Confrontando i primi 500 esoni classificati che mostrano attività AS e APA utilizzando DEXSeq e DiffSplice, abbiamo trovato una sovrapposizione di 310 esoni e 300 siti pA, rispettivamente, dimostrando la concordanza dei due approcci basati sugli esoni, che è stata dimostrata anche in uno studio precedente ²⁹. Si consiglia di utilizzare una combinazione di un approccio basato sull'esone (DEXSeq o diffSplice) e basato su eventi per il rilevamento e la classificazione completi dell'AS. Per APA, gli utenti possono scegliere DEXSeq o diffSplice: entrambi i metodi hanno dimostrato di funzionare bene in una vasta gamma di esperimenti di trascrittomica²⁹.

Nel preparare la libreria RNA-seq per un'analisi AS, è importante utilizzare un protocollo 8 di RNA-seq di massa specifico per il filamento, poiché una grande frazione di geni nei genomi dei vertebrati si sovrappone su diversi filamenti e^{un protocollo non} specifico del filamento non è in grado di distinguere queste regioni sovrapposte, confondendo la rilevazione finale dell'esone. Un'altra considerazione è la profondità di lettura, con analisi di splicing che richiedono un sequenziamento più approfondito rispetto al DGE, ad esempio 30-60 milioni di letture per campione, contro 5-25 milioni di letture per campione per DGE (https://sapac.support.illumina.com/bulletins/2017/04/considerations-for-rna-seq-read-length-and-coverage-.html). Tutti gli strumenti illustrati nel protocollo supportano sia dati di sequenziamento single-end che paired-end. Se vengono utilizzate solo annotazioni genetiche note per rilevare le letture di giunzione, è possibile utilizzare letture più brevi a terminazione singola (≥ 50 bp), sebbene il rilevamento de novo di nuove giunzioni di giunzione benefici da letture a estremità accoppiata e più lunghe (≥ 100 bp) letture^30,31. La scelta del protocollo di estrazione dell'RNA - selezione di poliA o deplezione di rRNA - dipende dalla qualità dell'RNA e dalla domanda sperimentale - vedi³¹ per una discussione. A seconda dei dettagli della costruzione della libreria, saranno necessarie modifiche agli script di esempio forniti qui per i parametri di allineamento alla lettura, conteggio delle funzionalità e rMATS. Nel calcolare i conteggi di lettura iniziali del livello di esone utilizzando featureCounts, o metodi simili, è necessario prestare attenzione a configurare correttamente le opzioni di funzione per conteggi e strandedness: in featureCounts, impostiamo l'argomento "strandSpecific" in modo appropriato per il protocollo RNA-seq strand-specific utilizzato; e per la quantificazione a livello di esone ci si aspetta che una lettura venga mappata sugli esoni adiacenti, quindi impostiamo il parametro allowMultiOverlap su TRUE. Per l'APA, ci sono diversi protocolli di sequenziamento finale 3'⁶ che variano nella posizione precisa dei picchi rispetto al sito pA. Per i nostri dati di esempio determiniamo che il picco è 60 bp a monte del sito pA, come mostrato dalla Figura 5, e questa analisi dovrà essere adattata per altri protocolli di sequenziamento finale 3'.

In questo protocollo limitiamo l'ambito alla discussione delle analisi differenziali a livello dei singoli esoni e degli eventi di splicing costituiti da combinazioni esone-introne adiacenti. Non discutiamo la classe di analisi basate sulla ricostruzione de novo delle isoforme come Cufflinks, Cuffdiff³², RSEM 33, Kallisto³⁴ che mirano a rilevare e quantificare l'espressione assoluta e relativa di intere isoforme alternative. I metodi basati sugli esoni e sugli eventi sono più sensibili per rilevare singoli eventi di splicing³⁰ e in molti casi forniscono tutte le informazioni necessarie per ulteriori analisi, senza la necessità di quantificazione a livello di isoforma.

L'ultima versione dei file sorgente in questo protocollo è disponibile all'indirizzo https://github.com/jiayuwen/AS_APA_JoVE

Gli autori non hanno nulla da rivelare.

Questo studio è stato supportato da una Future Fellowship dell'Australian Research Council (ARC) (FT16010043) e da ANU Futures Scheme.