Identificatie van alternatieve splicing en polyadenylation in RNA-seq-gegevens

Summary

Alternatieve splicing (AS) en alternatieve polyadenylation (APA) vergroten de diversiteit van transcriptisovormen en hun producten. Hier beschrijven we bioinformatische protocollen om bulk RNA-seq en 3'-eindsequencingtests te analyseren om AS en APA te detecteren en te visualiseren, variërend tussen experimentele omstandigheden.

Abstract

Naast de typische analyse van RNA-Seq om differentiële genexpressie (DGE) te meten over experimentele / biologische omstandigheden, kunnen RNA-seq-gegevens ook worden gebruikt om andere complexe regulerende mechanismen op exon-niveau te verkennen. Alternatieve splicing en polyadenylation spelen een cruciale rol in de functionele diversiteit van een gen door verschillende isovormen te genereren om genexpressie op post-transcriptioneel niveau te reguleren, en het beperken van analyses tot het hele genniveau kan deze belangrijke regulerende laag missen. Hier demonstreren we gedetailleerde stapsgewijze analyses voor identificatie en visualisatie van differentiële exon- en polyadenylation-sitegebruik onder omstandigheden, met behulp van Bioconductor en andere pakketten en functies, waaronder DEXSeq, diffSplice uit het Limma-pakket en rMATS.

Introduction

RNA-seq is in de loop der jaren op grote schaal gebruikt, meestal voor het schatten van differentiële genexpressie en genontdekking¹. Bovendien kan het ook worden gebruikt om het gebruik op exonniveau te schatten als gevolg van gen dat verschillende isovormen tot expressie brengt, waardoor het bijdraagt aan een beter begrip van genregulatie op post-transcriptioneel niveau. De meerderheid van de eukaryote genen genereren verschillende isovormen door alternatieve splicing (AS) om de diversiteit van mRNA-expressie te vergroten. AS-gebeurtenissen kunnen worden onderverdeeld in verschillende patronen: het overslaan van complete exonen (SE) waarbij een ("cassette") exon volledig uit het transcript wordt verwijderd samen met de flankerende introns; alternatieve (donor) 5' splice site selection (A5SS) en alternatieve 3' (acceptor) splice site selection (A3SS) wanneer twee of meer splice sites aanwezig zijn aan beide uiteinden van een exon; retentie van introns (RI) wanneer een intron wordt bewaard binnen het volwassen mRNA-transcript en wederzijdse uitsluiting van exongebruik (MXE) waarbij slechts één van de twee beschikbare exonen tegelijk kan worden bewaard ^2,3. Alternatieve polyadenylering (APA) speelt ook een belangrijke rol bij het reguleren van genexpressie met behulp van alternatieve poly (A) -sites om meerdere mRNA-isovormen te genereren uit een enkel transcript⁴. De meeste polyadenylation sites (pa's) bevinden zich in het 3' onvertaalde gebied (3' UTR's), waardoor mRNA-isovormen met diverse 3' UTR-lengtes worden gegenereerd. Aangezien de 3' UTR de centrale hub is voor het herkennen van regulerende elementen, kunnen verschillende 3' UTR-lengtes van invloed zijn op mRNA-lokalisatie, stabiliteit en translatie⁵. Er zijn een klasse van 3'-eindsequencingtests die zijn geoptimaliseerd om APA te detecteren die verschillen in de details van het protocol⁶. De hier beschreven pijplijn is ontworpen voor PolyA-seq, maar kan worden aangepast voor andere protocollen zoals beschreven.

In deze studie presenteren we een pijplijn van differentiële exon-analysemethoden ^7,8 (figuur 1), die kunnen worden onderverdeeld in twee brede categorieën: exon-gebaseerd (DEXSeq⁹, diffSplice¹⁰) en gebeurtenisgebaseerd (multivariate analyse van transcriptsplitsing (rMATS)¹¹). De op exon gebaseerde methoden vergelijken de vouwverandering over de omstandigheden van individuele exonen, met een maat voor de totale genplooiverandering om differentieel tot expressie gebracht exongebruik aan te roepen, en berekenen op basis daarvan een maat op genniveau van AS-activiteit. Op gebeurtenissen gebaseerde methoden gebruiken exon-intron-spanning junction reads om specifieke splicinggebeurtenissen zoals exon-overslaan of retentie van introns te detecteren en te classificeren, en onderscheiden deze AS-typen in de uitgang³. Deze methoden bieden dus complementaire weergaven voor een volledige analyse van AS ^12,13. We selecteerden DEXSeq (gebaseerd op het DESeq2¹⁴ DGE-pakket) en diffSplice (gebaseerd op het Limma¹⁰ DGE-pakket) voor de studie omdat ze tot de meest gebruikte pakketten voor differentiële splicinganalyse behoren. rMATS werd gekozen als een populaire methode voor event-based analyse. Een andere populaire op gebeurtenissen gebaseerde methode is MISO (Mixture of Isoforms)¹. Voor APA passen we de exon-gebaseerde aanpak aan.

Figuur 1. Analysepijplijn. Stroomdiagram van de stappen die in de analyse zijn gebruikt. Stappen omvatten: het verkrijgen van de gegevens, het uitvoeren van kwaliteitscontroles en leesuitlijning gevolgd door het tellen van reads met behulp van annotaties voor bekende exonen, introns en pA-sites, filteren om lage tellingen en normalisatie te verwijderen. PolyA-seq-gegevens werden geanalyseerd voor alternatieve pA-locaties met behulp van diffSplice/DEXSeq-methoden, bulk RNA-Seq werd geanalyseerd voor alternatieve splicing op exon-niveau met diffSplice/DEXseq-methoden en AS-gebeurtenissen geanalyseerd met rMATS. Klik hier om een grotere versie van deze figuur te bekijken.

De RNA-seq-gegevens die in dit onderzoek zijn gebruikt, zijn verkregen uit Gene Expression Omnibus (GEO) (GSE138691)¹⁵. We gebruikten muis RNA-seq-gegevens uit deze studie met twee conditiegroepen: wild-type (WT) en Muscleblind-like type 1 knockout (Mbnl1 KO) met elk drie replicaties. Om de analyse van het gebruik van differentiële polyadenylation-sites aan te tonen, verkregen we polyA-seq-gegevens van muizenembryoblasten (MEF's) (GEO Accession GSE60487)¹⁶. De gegevens hebben vier conditiegroepen: Wild-type (WT), Muscleblind-like type1/type 2 double knockout (Mbnl1/2 DKO), Mbnl 1/2 DKO met Mbnl3 knockdown (KD) en Mbnl1/2 DKO met Mbnl3 control (Ctrl). Elke conditiegroep bestaat uit twee replicaties.

	GEO Toetreding	SRA Run nummer	Voorbeeldnaam	Conditie	Nabootsen	Weefsel	Sequencing	Leeslengte
RNA-Seq	GSM4116218	SRR10261601	Mbnl1KO_Thymus_1	Mbnl1 knock-out	Rep 1	Thymus	Gekoppeld uiteinde	100 bp
	GSM4116219	SRR10261602	Mbnl1KO_Thymus_2	Mbnl1 knock-out	Rep 2	Thymus	Gekoppeld uiteinde	100 bp
	GSM4116220	SRR10261603	Mbnl1KO_Thymus_3	Mbnl1 knock-out	Rep 3	Thymus	Gekoppeld uiteinde	100 bp
	GSM4116221	SRR10261604	WT_Thymus_1	Wild type	Rep 1	Thymus	Gekoppeld uiteinde	100 bp
	GSM4116222	SRR10261605	WT_Thymus_2	Wild type	Rep 2	Thymus	Gekoppeld uiteinde	100 bp
	GSM4116223	SRR10261606	WT_Thymus_3	Wild type	Rep 3	Thymus	Gekoppeld uiteinde	100 bp
3P-Seq	GSM1480973	SRR1553129	WT_1	Wild type (WT)	Rep 1	Embryonale fibroblasten bij muizen (MEF's)	Single-end	40 bp
	GSM1480974	SRR1553130	WT_2	Wild type (WT)	Rep 2	Embryonale fibroblasten bij muizen (MEF's)	Single-end	40 bp
	GSM1480975	SRR1553131	DKO_1	Mbnl 1/2 dubbele knock-out (DKO)	Rep 1	Embryonale fibroblasten bij muizen (MEF's)	Single-end	40 bp
	GSM1480976	SRR1553132	DKO_2	Mbnl 1/2 dubbele knock-out (DKO)	Rep 2	Embryonale fibroblasten bij muizen (MEF's)	Single-end	40 bp
	GSM1480977	SRR1553133	DKOsiRNA_1	Mbnl 1/2 dubbele knock-out met Mbnl 3 siRNA (KD)	Rep 1	Embryonale fibroblasten bij muizen (MEF's)	Single-end	40 bp
	GSM1480978	SRR1553134	DKOsiRNA_2	Mbnl 1/2 dubbele knock-out met Mbnl 3 siRNA (KD)	Rep 2	Embryonale fibroblasten bij muizen (MEF's)	Single-end	36 bp
	GSM1480979	SRR1553135	DKONTsiRNA_1	Mbnl 1/2 dubbele knock-out met non-targeting siRNA (Ctrl)	Rep 1	Embryonale fibroblasten bij muizen (MEF's)	Single-end	40 bp
	GSM1480980	SRR1553136	DKONTsiRNA_2	Mbnl 1/2 dubbele knock-out met non-targeting siRNA (Ctrl)	Rep 2	Embryonale fibroblasten bij muizen (MEF's)	Single-end	40 bp

Tabel 1. Samenvatting van RNA-Seq en PolyA-seq datasets die voor de analyse zijn gebruikt.

Protocol

1. Installatie van tools en R-pakketten die bij de analyse worden gebruikt

1. Conda is een populaire en flexibele pakketbeheerder die een gemakkelijke installatie van pakketten met hun afhankelijkheden op alle platforms mogelijk maakt. Gebruik 'Anaconda' (conda package manager) om 'conda' te installeren, wat kan worden gebruikt om de tools/pakketten te installeren die nodig zijn voor de analyse.
2. Download 'Anaconda' volgens de systeemvereisten van https://www.anaconda.com/products/individual#Downloads en installeer het door de aanwijzingen in het grafische installatieprogramma te volgen. Installeer alle vereiste pakketten met behulp van 'conda' door het volgende te typen op de Linux-opdrachtregel.
   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. Als u alle R-pakketten wilt downloaden die in het protocol worden gebruikt, typt u de volgende code in de R-console (gestart op de Linux-opdrachtregel door 'R') of Rstudio-console te typen.
   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)
}
   OPMERKING: In dit computationele protocol worden opdrachten gegeven als R Notebook-bestanden (bestanden met de extensie ". Rmd"), R code bestanden (bestanden met de extensie ". R"), of Linux Bash shell scripts (bestanden met de extensie ".sh"). R Notebook (Rmd)-bestanden moeten in RStudio worden geopend met File| Open Bestand... en afzonderlijke codebrokken (wat R-opdrachten of Bash-shellopdrachten kunnen zijn) en voer vervolgens interactief uit door op de groene pijl rechtsboven te klikken. R-codebestanden kunnen worden uitgevoerd door te openen in RStudio, of op de Linux-opdrachtregel door vooraf te gaan met "Rscript", bijvoorbeeld Rscript example.R. Shell-scripts worden uitgevoerd op de Linux-opdrachtregel door het script vooraf te maken met het commando "sh", bijvoorbeeld .sh example.sh.

2. Alternatieve splicing (AS) analyse met behulp van RNA-seq

1. Gegevens downloaden en voorbewerking
   OPMERKING: De onderstaande codefragmenten zijn beschikbaar in het aanvullende codebestand "AS_analysis_RNASeq.Rmd", om de afzonderlijke stappen interactief te volgen, en worden ook geleverd als een bash-script om in batch op de Linux-opdrachtregel te worden uitgevoerd (SH downloading_data_preprocessing.sh).
   1. Het downloaden van de onbewerkte gegevens.
      1. Download de onbewerkte gegevens uit Sequence Read Archive (SRA) met de opdracht 'prefetch' uit SRA toolkit (v2.10.8)¹⁷. Geef de SRA-toetredings-ID's in volgorde in het volgende commando om ze parallel te downloaden met behulp van GNU parallel utility¹⁸. Als u SRA-bestanden van toetredings-id's van SRR10261601 naar SRR10261606 parallel wilt downloaden, gebruikt u het volgende op de Linux-opdrachtregel.
         ​seq 10261601 10261606 | parallel prefetch SRR{}
      2. Pak de fastq-bestanden uit het archief met behulp van de 'fastq-dump'-functie uit de SRA-toolkit. Gebruik GNU parallel en geef de namen van alle SRA-bestanden samen.
         ​parallel -j 3 fastq-dump --gzip --skip-technical --read-filter pass --dumpbase --split-e --clip --origfmt {} :::
      3. Download het referentiegenoom en de annotaties voor muis (Genome assembly GRCm39) van www.ensembl.org met behulp van het volgende op de Linux-opdrachtregel.
         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. Voorbewerking en mapping leest voor aan de genoomassemblage
      1. Kwaliteitscontrole. Beoordeel de kwaliteit van onbewerkte reads met FASTQC (FASTQ Quality Check v0.11.9)¹⁹. Maak een uitvoermap en voer fastqc parallel uit op meerdere invoer fasta-bestanden. Deze stap genereert een kwaliteitsrapport voor elk monster. Onderzoek de rapporten om ervoor te zorgen dat de kwaliteit van de gelezen teksten acceptabel is voordat u verdere analyse uitvoert. (Raadpleeg de gebruikershandleiding voor meer informatie over de rapporten op https://www.bioinformatics.babraham.ac.uk/projects/fastqc/)
         mkdir fastqc_out
parallel "fastqc {} -o fastqc_out" ::: $RAW_DATA/*.fastq.gz
         OPMERKING: Voer indien nodig adapter trimming uit met '^cutadapt'20 of 'trimmomatic'21 om sequencing in flankerende adapters te verwijderen, wat varieert op basis van de grootte van RNA-fragmenten en de leeslengte. In deze analyse hebben we deze stap overgeslagen omdat de fractie van de getroffen reads minimaal was.
      2. Lees uitlijning. De volgende stap in de voorbewerking omvat het in kaart brengen van de reads aan het referentiegenoom. Bouw eerst de index voor het referentiegenoom met behulp van de 'genomeGenerate'-functie van STAR²² en stem vervolgens de onbewerkte reads af op de referentie (als alternatief zijn vooraf gebouwde indexen beschikbaar op de STAR-website en kunnen rechtstreeks worden gebruikt voor uitlijning). Voer de volgende opdrachten uit op de Linux-opdrachtregel.
         #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
         OPMERKING: De STAR-aligner genereert en sorteert BAM-bestanden (Binary Alignment Map) voor elk voorbeeld na uitlijning. Bam-bestanden moeten worden gesorteerd voordat ze doorgaan met verdere stappen.
2. Exon-annotaties voorbereiden.
   1. Voer het aanvullende codebestand "prepare_mm_exon_annotation. R" met de gedownloade annotatie in GTF-indeling (Gene transfer format) om de annotaties voor te bereiden. Als u wilt uitvoeren, typt u het volgende op de Linux-opdrachtregel.
      Rscript prepare_mm_exon_annotation.R annotation.gtf
      OPMERKING: Het GTF-bestand bevat meerdere exon-vermeldingen voor verschillende isovormen. Dit bestand wordt gebruikt om de meerdere transcript-ID's voor elke exon "samen te vouwen". Het is een belangrijke stap om exon telbakken te definiëren.
3. Tellen Leest. De volgende stap is het tellen van het aantal reads dat is toegewezen aan verschillende transcripten / exonen. Zie Aanvullend bestand: "AS_analysis_RNASeq.Rmd".
   1. Laad vereiste bibliotheken:
      packages<- c("Rsubread","tidyverse", "magrittr", "EnhancedVolcano", "edgeR","openxlsx")
invisible(lapply(packages, library, character.only=TRUE))
   2. Laad het verwerkte annotatiebestand dat is verkregen uit de vorige stap (2.2).
      load("mm_exon_anno.RData")
   3. Lees alle bam-bestanden verkregen in stap 2.2.2 als invoer voor 'featureCounts' om reads te tellen. Lees de map met bam-bestanden door eerst elk bestand uit de map te vermelden dat eindigt op .bam. Gebruik 'featureCounts' uit het Rsubread-pakket dat bam-bestanden en verwerkte GTF-annotatie (referentie) als invoer gebruikt om een matrix van tellingen te genereren die aan elke functie zijn gekoppeld met rijen die exonen (functies) vertegenwoordigen en kolommen die voorbeelden vertegenwoordigen.
      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. Voer vervolgens niet-specifieke filtering uit om laag uitgedrukte exonen te verwijderen ("niet-specifiek" geeft aan dat de informatie over de experimentele toestand niet wordt gebruikt in de filtering, om selectiebias te voorkomen). Transformeer de gegevens van ruwe schaal naar tellingen per miljoen (cpm) met behulp van de cpm-functie uit 'edgeR'-pakket²³ en bewaar exonen met tellingen groter dan een instelbare drempel (voor deze dataset wordt één cpm gebruikt) in ten minste drie steekproeven. Verwijder ook genen met slechts één exon.
      # 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, ]
      OPMERKING: Controleer de vereiste parameters voor featureCounts bij het gebruik van verschillende gegevens, bijvoorbeeld voor single-end reads, stel 'isPairedEnd = FALSE' in. Raadpleeg de RSubread-gebruikershandleiding om de opties voor uw gegevens te kiezen en raadpleeg het gedeelte Discussie hieronder.
4. Differentiële splicing en exon gebruiksanalyse. We beschrijven twee alternatieven voor deze stap: DEXSeq en DiffSplice. Beide kunnen worden gebruikt en vergelijkbare resultaten geven. Voor consistentie selecteert u DEXSeq als u de voorkeur geeft aan een DESeq2-pakket voor DGE en DiffSplice gebruikt voor een op Limma gebaseerde DGE-analyse. Zie Aanvullend dossier: "AS_analysis_RNASeq.Rmd".
   1. Het DEXSeq-pakket gebruiken voor differentiële exonanalyse.
      1. Laad de bibliotheek en maak een voorbeeldtabel om het experimentele ontwerp te definiëren.
         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")))
         OPMERKING: De rijnamen moeten consistent zijn met de bam-bestandsnamen die door featureCounts worden gebruikt om de reads te tellen. sampleTable bestaat uit details van elk voorbeeld, waaronder: bibliotheektype en conditie. Dit is nodig om de contrasten of testgroep te definiëren voor het detecteren van differentieel gebruik.
      2. Bereid het exon-informatiebestand voor. Exon-informatie in de vorm van GRanges-objecten (genomische bereiken) (https://bioconductor.org/packages/release/bioc/html/GenomicRanges.html) is vereist als invoer om het DEXSeq-object in de volgende stap te maken. Vergelijk het gen Ids met de gelezen tellingen om een exoninfo-object te maken.
         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. Maak een DEXSeq-object met de functie DEXSeqDataSet. Het DEXSeq-object verzamelt leestellingen, exon-functiegegevens en voorbeeldinformatie. Gebruik de leestellingen die in stap 3 zijn gegenereerd en de exon-informatie die is verkregen uit de vorige stap om het DEXSeq-object te maken op basis van de tellingsmatrix. Het argument sampleData neemt een gegevensframe-invoer die de voorbeelden definieert (en hun kenmerken: bibliotheektype en -voorwaarde), 'ontwerp' gebruikt sampleData om een ontwerpmatrix voor de differentiële tests te genereren met behulp van modelformulenotatie. Merk op dat een significante interactieterm, condition:exon, aangeeft dat de fractie van reads over een gen dat op een bepaald exon valt, afhankelijk is van de experimentele conditie, d.w.z. er is AS. Zie de DEXSeq-documentatie voor een volledige beschrijving van het instellen van de modelformule voor complexere experimentele ontwerpen. Voor functie-informatie zijn exon-id's, overeenkomstige genen en transcripten vereist.
         ​dxd<-DEXSeqDataSet(countData$counts,sampleData=sampleTable, design=~sample+exon+condition:exon,featureID=exoninfo$ExonID,groupID=exoninfo$GeneID,featureRanges=exoninfo, transcripts=transcripts_l)
      4. Normalisatie en dispersieschatting. Voer vervolgens normalisatie uit tussen monsters en schat de variantie van de gegevens, vanwege zowel Poisson-tellingsruis van de discrete aard van RNA-seq als biologische variabiliteit, met behulp van de volgende commando's.
         ​dxd %<>% estimateSizeFactors %>% estimateDispersions %T>% plotDispEsts
      5. Test op differentieel gebruik. Test na de schatting van de variatie op differentieel exongebruik voor elk gen en genereer de resultaten.
         dxd%<>%testForDEU%>%estimateExonFoldChanges(fitExpToVar=
"condition")#Estimate fold changes
dxr=DEXSeqResults(dxd)
      6. Visualiseer splicinggebeurtenissen voor geselecteerde genen met behulp van de volgende opdracht.
         plotDEXSeq(dxr,"Wnk1", displayTranscripts=TRUE, splicing=TRUE,legend
=TRUE,cex.axis=1.2,cex=1.3,lwd=2)
         Onderzoek het R Notebook-bestand "AS_analysis_RNASeq.Rmd" om extra plots te genereren voor genen van belang en om vulkaanplots op verschillende drempels te genereren.
   2. Met behulp van diffSplice van Limma om differentiële splicing te identificeren. Volg het R Notebook bestand "AS_analysis_RNASeq.Rmd". Zorg ervoor dat de stappen 2.1-2.3 zijn gevolgd om invoerbestanden voor te bereiden voordat u verder gaat.
      1. Bibliotheken laden
         library(limma)
library(edgeR)
      2. Niet-specifieke filtering. Extraheer de matrix van afgelezen tellingen verkregen in 2.3. Maak een lijst met functies met behulp van de functie 'DGEList' van het edgeR-pakket, waarbij rijen genen vertegenwoordigen en kolommen monsters vertegenwoordigen.
         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, "\\,")
         OPMERKING: Als niet-specifieke filterstap worden tellingen gefilterd door cpm < 1 in x uit n monsters, waarbij x het minimale aantal replicaties is in elke omstandigheid. n = 6 en x = 3 voor dit voorbeeld gegevens.
      3. Normaliseer de tellingen in monsters, met de functie 'calcNormFactors' uit het 'edgeR'-pakket met behulp van bijgesneden gemiddelde van M-waarden (TMM-normalisatiemethode)²⁴ Het berekent schaalfactoren om bibliotheekgrootten aan te passen.
         ​dge<-calcNormFactors(dge)
      4. Gebruik sampleTable zoals gegenereerd in stap 2.4.1.1 en maak de ontwerpmatrix. De ontwerpmatrix kenmerkt het ontwerp. Zie de Limma User Guide (https://www.bioconductor.org/packages/devel/bioc/vignettes/limma/inst/doc/usersguide.pdf) hoofdstukken 8 &9 voor meer informatie over ontwerpmatrices voor meer geavanceerde experimentele ontwerpen.
         Treat<- factor(sampleTable$condition)
design<- model.matrix(~0+Treat)
colnames(design) <- levels(Treat)
      5. Pas een lineair model per exon aan. Voer de 'voom'-functie van het 'limma'-pakket uit om RNA-seq-gegevens te verwerken om variantie te schatten en precisiegewichten te genereren om te corrigeren voor Poisson-tellingsruis, en transformeer de exon-niveautellingen naar log2-tellingen per miljoen (logCPM). Voer vervolgens lineaire modellering uit met behulp van de 'lmfit'-functie om lineaire modellen aan te passen aan de expressiegegevens voor elk exon. Bereken empirische Bayes-statistieken voor aangepast model met behulp van de 'eBayes'-functie om differentiële exonexpressie te detecteren. Definieer vervolgens een contrastmatrix voor de experimentele vergelijkingen van belang. Gebruik 'contrasten.fit' om coëfficiënten en standaardfouten voor elk vergelijkingspaar te verkrijgen.
         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. Differentiële splicing analyse. Voer 'diffSplice' uit op het aangepaste model om de verschillen in exongebruik van genen tussen wild-type en knock-out te testen en de best gerangschikte resultaten te verkennen met behulp van de 'topSplice'-functie: test = "t" geeft een rangschikking van AS-exonen, test = "simes" geeft een rangschikking van genen.
         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. Visualisatie. Plot de resultaten met de 'plotSplice'-functie, waardoor het gen van belang is voor het geneid-argument. Sla de topresultaten gesorteerd op log Vouw verandering in een object en genereer een vulkaanplot om de exonen tentoon te stellen.
         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. rMATS gebruiken
      1. Zorg ervoor dat de nieuwste versie van rMATS v4.1.1 (ook bekend als rMATS turbo vanwege de kortere verwerkingstijd en minder vereisten van het geheugen) is geïnstalleerd met behulp van conda of github (https://github.com/Xinglab/rmats-turbo/releases/download/v4.1.1/rmats_turbo_v4_1_1.tar.gz) in de werkmap. Volg paragraaf 4.3 in "AS_analysis_RNASeq.Rmd".
      2. Ga naar de map met bam-bestanden die zijn verkregen na toewijzing en bereid tekstbestanden voor, zoals vereist door rMATS, voor de twee voorwaarden door de naam van bam-bestanden (samen met het pad) gescheiden door ',' te kopiëren. De volgende opdrachten moeten worden uitgevoerd op de Linux-opdrachtregel:
         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. Voer rmats.py uit met de twee invoerbestanden die in de vorige stap zijn gegenereerd, samen met het GTF-bestand dat is verkregen in 2.1.1.3. Dit genereert een uitvoermap 'rmats_out' met tekstbestanden die statistieken (p-waarden en inclusieniveaus) beschrijven voor elke splicinggebeurtenis afzonderlijk.
         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
         OPMERKING: De referentieannotatie in de vorm van een GTF-bestand is ook vereist. Controleer de parameters als de gegevens single-end zijn en wijzig de optie -t dienovereenkomstig.
      4. RMATS-resultaten verkennen. Gebruik Bioconductor pakket '^maser'25 om de rMATS resultaten te verkennen. Laad de JCEC-tekstbestanden (Junction and Exon Counts) in het 'maser'-object en filter het resultaat op basis van dekking door ten minste vijf gemiddelde leeswaarden per splicinggebeurtenis op te nemen.
         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. Visualisatie van rMATS-resultaten. Selecteer de significante splicinggebeurtenissen bij False Discovery Rate (FDR) 10% en minimaal 10% verandering in Percent Spliced In (deltaPSI) met behulp van de functie 'topEvents' uit het 'maser'-pakket. Controleer vervolgens de gengebeurtenissen op individuele genen van belang (monstergen-Wnk1) en plot PSI-waarden voor elke splicinggebeurtenis van dat gen. Genereer een vulkaanplot door het gebeurtenistype op te geven.
         #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. Genereer Sashimi-plots voor het splicinggebeurtenissenresultaat verkregen met rMATS in de vorm van tekstbestanden met behulp van het 'rmats2shahimiplot'-pakket. Voer het python-script uit op de Linux-opdrachtregel.
         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
         OPMERKING: Dit proces kan tijdrovend zijn omdat het de Sashimi-plot genereert voor alle resultaten in het gebeurtenissenbestand. Kies de topresultaten (gennamen en exonen) zoals weergegeven door de topEvents-functie van 'maser' en visualiseer de bijbehorende Sashimi-plot.

3. Alternatieve polyadenylering (APA) analyse met behulp van 3' end sequencing

1. Gegevens downloaden en voorbewerking
   OPMERKING: Raadpleeg het aanvullende R Notebook-bestand "APA_analysis_3PSeq_notebook. Rmd" voor de volledige opdrachten voor het downloaden en voorbewerkingsstappen van gegevens, of voer het aanvullende bash-bestand "APA_data_downloading_preprocessing.sh" uit op de Linux-opdrachtregel.
   1. Download gegevens van SRA met de Toetredings-Ids (1553129 tot 1553136).
   2. Trim adapters en omgekeerde aanvulling om de zinstrengvolgorde te verkrijgen.
      OPMERKING: Deze stap is specifiek voor de gebruikte PolyA-seq-test.
   3. Kaart leest naar muis genoom assemblage met behulp van bowtie aligner²⁶.
2. Het voorbereiden van pA-sitesannotaties.
   OPMERKING: De verwerking van het pA-siteannotatiebestand wordt eerst uitgevoerd met behulp van het aanvullende R Notebook-bestand "APA_analysis_3PSeq_notebook. Rmd" (2.1 - 2.6), en dan met behulp van bash bestand "APA_annotation_preparation.sh".
   1. Download de annotatie van pA-sites uit de PolyASite 2.0-database⁶.
   2. Selecteer pA-siteannotaties om 3'-untranslated region (UTR) pA-sites te behouden, die zijn geannoteerd als Terminal Exon (TE) of 1000 nt downstream van een geannoteerde terminal exon (DS) voor downstream-analyse.
   3. Verkrijg pA-sitepieken. Anker op elke pA-splitsingsplaats en visualiseer de gemiddelde leesdekking met bedgereedschappen en deeptools^27,28. De resultaten toonden aan dat de pieken van de in kaart gebrachte reads voornamelijk verspreid waren binnen ~60 bp stroomopwaarts van de splitsingsplaatsen (figuur 5 en aanvullende figuur 5). Daarom werden de coördinaten van pA-sites uitgebreid van het annotatiebestand tot 60 bp stroomopwaarts van hun splitsingsplaatsen. Afhankelijk van het specifieke 3'-eindsequencingprotocol dat wordt gebruikt, moet deze stap worden geoptimaliseerd voor andere assays dan PolyA-seq.
3. Tellen leest
   1. Bereid het annotatiebestand van pA-sites voor.
      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. Pas 'featureCounts' toe om onbewerkte tellingen te verkrijgen. Sla de teltabel op als het RData-bestand "APA_countData.Rdata" voor APA-analyse met behulp van verschillende hulpprogramma's.
      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")
      OPMERKING: Wees bewust van het wijzigen van een van de parameters die worden vermeld in de functie 'featureCounts'. Wijzig de parameter 'strandSpecific' om ervoor te zorgen dat deze consistent is met de sequencingrichting van de gebruikte 3'-eindsequencingtest (empirisch gezien zal het visualiseren van de gegevens in een genoombrowser over genen op plus- en minstrengen dit verduidelijken).
   3. Pas niet-specifieke filtering van countData toe. Filtering kan de statistische robuustheid in differentiële pA-sitegebruikstests aanzienlijk verbeteren. Ten eerste hebben we die genen verwijderd met slechts één pA-site, waarop differentieel pA-sitegebruik niet kan worden gedefinieerd. Ten tweede passen we niet-specifieke filtering toe op basis van dekking: tellingen worden gefilterd door cpm minder dan 1 in x uit n monsters, waarbij x het minimale aantal replicaties is in elke omstandigheid. N = 8 en x = 2 voor dit voorbeeld gegevens.
      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. Differentiële polyadenylation site gebruiksanalyse met behulp van DEXSeq en diffSplice pipelines.
   1. DeXSeq-pakket gebruiken
      OPMERKING: Aangezien een contrastmatrix niet kan worden gedefinieerd voor de DEXSeq-pijplijn, moet de differentiële APA-analyse van elke twee experimentele omstandigheden afzonderlijk worden uitgevoerd. De differentiële APA-analyse van de aandoening WT en de voorwaarde DKO wordt als voorbeeld uitgevoerd om de procedure uit te leggen. Raadpleeg het aanvullend dossier "APA_analysis_3PSeq_notebook. Rmd" voor de stapsgewijze workflow van deze sectie en de differentiële APA-analyse van andere contrasten.
      1. Laad de bibliotheek en maak een voorbeeldtabel om het experimentele ontwerp te definiëren.
         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. Bereid het informatiebestand van pA-sites voor met behulp van Bioconductor-pakket 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. Gebruik de leestellingen die zijn gegenereerd in stap 3.3 en de pA-sitegegevens die zijn verkregen uit de vorige stap om het DEXSeq-object te maken.
         # 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. Definieer het contrastpaar door de niveaus van voorwaarden in het DEXSeq-object te definiëren.
         dxd1$condition<- factor(dxd1$condition, levels= c("WT", "DKO"))
# The contrast pair will be "DKO - WT"
      5. Normalisatie en dispersieschatting. Net als bij RNA-seq-gegevens, voeren 3'-eindsequencinggegevens normalisatie uit tussen monsters (kolomgewijze mediaan van verhoudingen voor elk monster) met behulp van de functie 'estimateSizeFactors' en schatten de variatie van de gegevens met behulp van de functie 'estimateDispersions', en visualiseer vervolgens het dispersieschattingsresultaat met behulp van de functie 'plotDispEsts'.
         ​dxd1 %<>% estimateSizeFactors %>% estimateDispersions %T>% plotDispEsts
      6. Differentiële pA-sitegebruikstests voor elk gen met behulp van de functie 'testForDEU', en schatten vervolgens de vouwverandering van pA-sitegebruik met behulp van de functie 'estimateExonFoldChanges'. Controleer de resultaten met behulp van de functie 'DEXSeqResults' en stel 'FDR < 10%' in als criterium voor significant differentiële pA-sites.
         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. Visualisatie van de differentiële pA-sitegebruiksresultaten met behulp van differentiële APA-plots gegenereerd door de functie 'plotDEXSeq' en vulkaanplot door de functie '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. Met behulp van diffSplice pakket. Raadpleeg het aanvullende R Notebook-bestand "APA_analysis_3PSeq_notebook. Rmd" voor de stapsgewijze workflow van deze sectie.
      1. Definieer interessante contrasten voor differentiële pA-gebruiksanalyse.
         OPMERKING: Deze stap moet worden uitgevoerd na de constructie en verwerking van het DGEList-object, dat is opgenomen in het R Notebook-bestand "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. Visualiseer het resultaat van interessante contrasten (hier "DKO - WT") met behulp van differentiële APA-plots door de functie 'plotSplice' en vulkaanplots met de functie 'EnhancedVolcano'. Raadpleeg het R Notebook-bestand "APA_analysis_3PSeq_notebook. Rmd" 4.2.7 - 4.2.9 voor de visualisatie van andere contrastparen.
         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

Na het uitvoeren van de bovenstaande stapsgewijze workflow, worden de AS- en APA-analyse-uitvoer en representatieve resultaten in de vorm van tabellen en gegevensplots weergegeven, die als volgt worden gegenereerd.

ALS:
De belangrijkste output van de AS-analyse (aanvullende tabel 1 voor diffSplice; Tabel 2 voor DEXSeq) is een lijst van exonen die differentieel gebruik tussen omstandigheden laten zien, en een lijst van genen die significante algehele splicingactiviteit van een of meer van de samenstellende exonen vertonen, gerangschikt naar statistische significantie. Aanvullende tabel 1, tabblad 2 toont significante exonen, met kolommen met differentiële FC van exon versus rust, per exon niet-gecorrigeerde p-waarde en aangepaste p-waarde (Benjamini-Hockberg-correctie). Drempels op de aangepaste p-waarden geven een set exonen met gedefinieerde FDR. Aanvullende tabel 1, tabblad 3 toont een gerangschikte lijst van genen die significantie van de totale splicingactiviteit vertonen, met een kolom met de voor genniveau gecorrigeerde p-waarde berekend met behulp van de Simes-methode. Vergelijkbare gegevens zijn weergegeven in tabel 2 voor DEXSeq. Aanvullende figuur 1 en aanvullende figuur 2 tonen differentiële splicingpatronen in Mbnl1-, Tcf7- en Lef1-genen die experimenteel zijn gevalideerd in het gepubliceerde artikel dat met de gegevens is gepresenteerd¹⁵. De auteurs hebben experimentele validatie van vijf genen aangetoond: Mbnl1, Mbnl2, Lef1, Tcf7 en Ncor2. Onze aanpak detecteerde differentiële splicing pattens in al deze genen. Hier presenteren we de FDR-niveaus voor elk gen met respectievelijk DEXSeq, diffSplice en rMATS zoals verkregen in aanvullende tabellen 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) en Ncor2 (9.2E-11, 1.6E-06, 0).

Figuur 2 toont een vergelijking tussen outputs verkregen uit drie verschillende tools en illustreert alternatieve splicingpatronen in het Wnk1-gen . De vulkaanplots zijn weergegeven in figuur 2A (diffSplice) en figuur 2B (DEXSeq). Nog eens drie hoog gerangschikte genen worden getoond in aanvullende figuur 1 (diffSplice) en aanvullende figuur 2 (DEXSeq). De y-as toont statistische significantie (-log10 P-waarden) en de x-as toont de effectgrootte (vouwverandering). Genen in het kwadrant linksboven of rechtsboven wijzen op substantiële FC en sterk statistisch bewijs van echte verschillen.

Figuur 2. Vergelijking van alternatieve splicingresultaten verkregen uit diffSplice, DEXSeq en rMATS. |
(A) Vulkaan plot (links) van RNA-Seq van Limma diffSplice analyse: De x-as toont log exon vouw verandering; de y-as toont -log₁₀ p-waarde. Elk punt komt overeen met een exon. Horizontale stippellijn bij p-waarde=1E-5; verticale stippellijnen bij tweevoudige verandering (FC). Rode exonen vertonen een aanzienlijke FC en statistische significantie. Differentiële splicing plot (rechts): Splicing patronen worden getoond voor een voorbeeld gen Wnk1 waarbij de x-as exon id per transcript toont; de y-as toont exon relatieve log fold verandering (het verschil tussen de logFC van het exon en de totale logFC voor alle andere exonen). Exonen die in rood zijn gemarkeerd, tonen statistisch significante differentiële expressie (FDR < 0,1).
(B) Vulkaanplot (links) en Differentieel exongebruik (rechts) van RNA-Seq verkregen uit DEXSeq-analyse. Wnk1-gen vertoont significant differentieel gebruik van exonen tussen WT en Mbnl1 knock-out gemarkeerd in roze, die overeenkomen met dezelfde differentiële exonen in (A).
(C) Vulkaan plot (links) en Sashimi plot (rechts) voor Wnk1 verkregen uit rMATS analyse. Vulkaanplot met significante overgeslagen (cassette) exon (SE) gebeurtenis in wild-type in vergelijking met knock-out bij 10% FDR met verandering in procent gesplitst in (PSI of ΔΨ) waarden > 0,1. De x-as toont verandering in PSI-waarden tussen omstandigheden en de y-as toont log P-waarde. De Sashimi-plot toont een overgeslagen exongebeurtenis in het Wnk1-gen , overeenkomend met een significant differentieel exon in (A) en (B). Elke rij vertegenwoordigt een RNA-Seq-monster: drie replicaties van wild-type en Mbnl1 knock-out. De hoogte toont de leesdekking in RPKM en de verbindingsbogen geven junction reads over exonen weer. Geannoteerde genmodel alternatieve isovormen worden onderaan de plot weergegeven. Het onderste paneel van C illustreert junction reads die worden gebruikt om de PSI-statistiek te berekenen.
(D) Venn-diagram waarin het aantal significante differentiële exonen wordt vergeleken dat met de verschillende methoden wordt verkregen. Klik hier om een grotere versie van deze figuur te bekijken.

Figuur 2 Een (rechterpaneel) toont een schematische weergave van exonverschillen van een van de best gerangschikte genen, met logFC op de y-as en exon-nummer op de x-as. Dit voorbeeld toont een cassette exon variërend tussen voorwaarden voor gen Wnk1. De differentiële exon-gebruiksplot van DEXSeq toont bewijs van differentiële splicing op vijf exon-locaties in de buurt van Wnk1.6.45. De gemarkeerde exonen in roze worden waarschijnlijk uitgesplitst in Mbnl1 KO-monsters in vergelijking met WT. Deze exonen zijn complementair aan de resultaten verkregen door diffSplice die een vergelijkbaar patroon vertoont op de specifieke genomische positie. Meer voorbeelden zijn weergegeven in aanvullende figuur 1 en aanvullende figuur 2. Een meer gedetailleerd beeld om interessante resultaten te bevestigen kan worden gegeven door dekking (wiggle) tracks in RPM (Reads per million) eenheden in de UCSC (University of Santa Cruz) of IGV (Integrative Genomics Viewer) genoombrowsers (niet getoond) te vergelijken, samen met visuele correlatie met andere sporen van belang, zoals bekende genmodellen, conservering en andere genoombrede testen.

De rMATS-uitvoertabel bevat significante alternatieve splicinggebeurtenissen gecategoriseerd op type (aanvullende tabel 3). Figuur 2C toont een vulkaanplot van genen die alternatief gesplitst zijn, met de effectgrootte gemeten door de differentiële "procent gesplitst in" (PSI of ΔΨ) statistiek van¹¹.

PSI verwijst naar het percentage reads dat consistent is met de opname van een cassette exon (d.w.z. leest mapping naar de cassette exon zelf of junction reads overlappen de exon) vergeleken met reads consistent met exon uitsluiting, d.w.z. junction reads over aangrenzende upstream en downstream exonen (het onderste paneel van figuur 2C). Het rechterpaneel van figuur 2C toont de sashimi-plot van het Wnk1-gen met differentiële splicinggebeurtenis bedekt met dekkingssporen voor het gen, met een overgeslagen exon in Mbnl1 KO. Bogen die exonen verbinden, tonen het aantal junctielezingen (leest het kruisen van een uitgesplitst intron). Verschillende tabbladen van aanvullende tabel 3 tonen significante lezingen van elk type gebeurtenis dat knooppunten met exongrenzen overspant (knooppunttellingen en exontellingen (JCEC)). Figuur 2D vergelijkt de significante differentieel gesplitste exonen die door de drie gereedschappen worden gedetecteerd.

Figuur 3. Alternatieve splicinggebeurtenissen verkregen door rMATS-analyse. A) Soorten AS-evenementen. Deze figuur is gebaseerd op rMATS-documentatie¹¹ waarin de splicinggebeurtenis wordt uitgelegd met constitutieve en alternatief gesplitste exonen. Gelabeld met het nummer van elke gebeurtenis op FDR 10%. B) Sashimi-plot van Add3-gen vertoont overgeslagen exon (SE). C) Sashimi plot van Baz2b gen vertoont alternatieve 5 'splice site (A5SS). D) Sashimi-plot van het Lsm14b-gen met alternatieve 3 'splice-site (A3SS). E) Sashimi-plot van Mta1-gen dat elkaar uitsluitende exonen (MXE) vertoont. F) Sashimi-plot van Arpp21-gen met behoud van intron (RI). Rode rijen vertegenwoordigen drie replica's van wild-type en oranje rijen vertegenwoordigen Mbnl1 knock-out replica's. De x-as komt overeen met genomische coördinaten en strenginformatie, y-as toont dekking in RPKM. Klik hier om een grotere versie van deze figuur te bekijken.

Figuur 3 illustreert soorten splicinggebeurtenissen SE, A5SS, A3SS, MXE en RI met behulp van Sashimi-plots van de belangrijkste significante genen van die gebeurtenissen. Bij het vergelijken van de drie replicaties van zowel WT- als Mbnl1_KO, werden in totaal 1272 SE-gebeurtenissen, 130 A5SS, 116 A3SS,215 MXE- en 313 RI-gebeurtenissen gedetecteerd bij FDR 10%. Sashimi plot illustreert het type gebeurtenis met behulp van topgenen als voorbeeld.

APA:
De output van de APA-analyse is vergelijkbaar met de AS-analyse op exon-niveau. Een tabel met topgenen gerangschikt naar differentiële APA-activiteit in de 3'UTR wordt verstrekt (aanvullende tabel 4 en aanvullende tabel 5). Figuur 4A toont de vulkaanplots van genen door differentiële APA-activiteit in 3'UTR's gegenereerd met behulp van diffSplice en DEXSeq afzonderlijk. Figuur 4B toont de Venn-plot waarin de significant verschillende pA-sitegebruiksresultaten van verschillende pijpleidingen worden vergeleken. Figuur 4C en 4D tonen de diagrammatische weergave van differentieel pA-sitegebruik in de 3'UTR van gen Fosl2 (figuur 4C) en Papola (figuur 4D) gegenereerd met behulp van zowel diffSplice als DEXSeq, die experimenteel zijn gevalideerd om significante distale tot proximale verschuiving (Fosl2) en significante proximale tot distale verschuiving (Papola) van pA-sitegebruik in DKO te laten zien, respectievelijk. Meer voorbeelden zijn weergegeven in aanvullende figuur 3 en aanvullende figuur 4.

Figuur 4. Alternatieve polyadenylation plots door diffSplice en DEXSeq. A) Vulkaanplots van PolyA-seq-gegevens gegenereerd met behulp van diffSplice en DEXSeq. X-as toont log pA site fold verandering; y-as toont -log₁₀ p-waarde. Elk punt komt overeen met een pA-site. Horizontale stippellijn bij p-waarde=1E-5; verticale stippellijnen bij 2-voudige FC. Rode exonen vertonen een aanzienlijke FC en statistische significantie. B) Venn-plot waarin de significante differentiële pA-gebruiksresultaten van verschillende pijpleidingen worden vergeleken. C-D) Differentiële APA-plots gegenereerd met behulp van diffSplice en DEXSeq die de proximale, interne en distale pA-sites voor het Fosl2 - en Papola-gen weergeven. Cijfers worden gegenereerd door dezelfde functie als figuur 2 (B), maar met pA-locaties die exonen vervangen. Klik hier om een grotere versie van deze figuur te bekijken.

Figuur 5 is een diagnostische plot om de verwachte leesverdeling rond geannoteerde pA-splitsingsplaatsen voor de gebruikte PolyA-seq-test te bevestigen. Het toont de gemiddelde dekking in flankerende gebieden verankerd op bekende pA-splitsingsplaatsen op genoombreed niveau. In dit geval wordt de verwachte stapelingen stroomopwaarts van de sites gevisualiseerd. De leesverdelingen verankerd op pA-locaties voor alle PolyA-seq-monsters zijn weergegeven in aanvullende figuur 5.

Figuur 5. Mean coverage plot rond pA cleavage sites. De splitsingsplaats voor PolyA-seq-gegevens wordt weergegeven door een verticale stippellijn. X-as toont de basispositie ten opzichte van pA-splitsingsplaatsen, tot 100 nucleotiden stroomopwaarts en stroomafwaarts; y-as toont de gemiddelde leesdekking over alle pA-splitsingssites, genormaliseerd door bibliotheekgrootte in CPM. Klik hier om een grotere versie van deze figuur te bekijken.

De differentiële APA-resultaten van verschillende contrasten die door dezelfde pijplijn worden gegenereerd, kunnen worden vergeleken en geverifieerd door de leesdekking van representatieve significant differentiële pA-sites in de genoombrowser te visualiseren. Figuur 6A is de Venn-plot waarin het significant differentiële pA-sitegebruik van verschillende contrasten van diffSplice wordt vergeleken. Figuur 6B-D zijn de IGV-snapshots van de leesdekking op pA-locaties voor verschillende genen, die de patronen laten zien die consistent zijn met die ontdekt in de APA-analyse met behulp van diffSplice. Figuur 6B valideert de significante proximale tot distale verschuiving van pA-sitegebruik voor gen Paip2, dat uniek wordt gedetecteerd in het contrast DKO versus WT, maar niet in andere twee contrasten KD versus WT en Ctr versus WT. Figuur 6C valideert de significante distale naar proximale verschuiving van pA-sitegebruik voor gen Ccl2 uniek gedetecteerd in het contrast KD versus WT, terwijl figuur 6D het significante differentiële pA-gebruik van alle contrasten voor gen Cacna2d1 valideert.

Figuur 6. Contrastvergelijking en verificatie van diffSplice-resultaten. A) Venn-diagram waarin significante differentiële pA-sitegebruiksresultaten van verschillende contrasten worden vergeleken die zijn verkregen uit diffSplice. B-D) IGV snapshot visualiseert pA pieken dekking van genen Paip2, Ccl2 en Cacna2d1 over omstandigheden. Elke track vertegenwoordigt de leesdekking in een specifieke toestand. Klik hier om een grotere versie van deze figuur te bekijken.

Aanvullende figuur 1. RNA-Seq analyse van differentiële splicing met Limma diffSplice. (A) Vulkaan plot van RNA-Seq van Limma diffSplice analyse: De x-as toont log exon vouw verandering; de y-as toont -log₁₀ p-waarde. Elk punt komt overeen met een exon. Horizontale stippellijn bij p-waarde=1E-5; verticale stippellijnen bij tweevoudige verandering (FC). Rode exonen vertonen een aanzienlijke FC en statistische significantie. (B-D) Differentiële splicingplots: Splicingpatronen worden bijvoorbeeld weergegeven met respectievelijk genen Mbnl1, Tcf7 en Lef1, waarbij de x-as exon-id per transcript toont; de y-as toont exon relatieve log fold verandering (het verschil tussen de logFC van het exon en de totale logFC voor alle andere exonen). Exonen die in rood zijn gemarkeerd, tonen statistisch significante differentiële expressie (FDR < 0,1). Klik hier om dit bestand te downloaden.

Aanvullende figuur 2. RNA-Seq analyse van differentieel exon gebruik met DEXSeq. (A) Vulkaan plot. (B-D) Differentieel exongebruik van RNA-Seq verkregen uit DEXSeq-analyse. Genen Mbnl1, Tcf7 en Lef1 vertonen respectievelijk significant differentieel gebruik van exonen tussen WT en Mbnl1 knock-out gemarkeerd in roze, die overeenkomen met dezelfde differentiële exonen in aanvullende figuur 1. Klik hier om dit bestand te downloaden.  

Aanvullende figuur 3. Alternatieve polyadenylation plots door diffSplice. A) Vulkaanplots van PolyA-seq-gegevens gegenereerd met behulp van diffSplice in drie contrastparen verkregen uit de muis PolyA-seq-gegevens, waaronder dubbele knock-out (DKO) versus wild-type (WT), knock-down (KD) versus WT, en controle (Ctrl) versus WT. X-as toont log pA site fold verandering; y-as toont -log₁₀ p-waarde. Elk punt komt overeen met een pA-site. Horizontale stippellijn bij p-waarde=1E-5; verticale stippellijnen bij 2-voudige FC. Rode pA-sites vertonen een aanzienlijke FC- en statistische significantie. B) Differentiële APA-plots gegenereerd met behulp van diffSplice die de proximale, interne en distale pA-sites toont voor de hoog gerangschikte genen S100a7a, Tpm1 en Smc6. Klik hier om dit bestand te downloaden.  

Aanvullende figuur 4. Differentiële pA-gebruiksanalyse door DEXSeq-pijplijn. A) Vulkaanplots van PolyA-seq-gegevens gegenereerd met behulp van DEXSeq in drie paren verkregen uit de muis PolyA-seq-gegevens, waaronder dubbele knock-out (DKO) versus wild-type (WT), knock-down (KD) versus WT, en controle (Ctrl) versus WT. X-as toont log pA site fold change; y-as toont -log₁₀ p-waarde. Elk punt komt overeen met een pA-site. Horizontale stippellijn bij p-waarde=1E-5; verticale stippellijnen bij 2-voudige FC. Rode pA-sites vertonen een aanzienlijke FC- en statistische significantie. B) Differentiële APA-plots gegenereerd met behulp van DEXSeq die de proximale, interne en distale pA-sites tonen voor de hoog gerangschikte genen S100a7a, Tpm1 en Smc6. Klik hier om dit bestand te downloaden.  

Aanvullende figuur 5. Gemiddelde dekking plot en heatmaps rond pA splitsing sites.  Dekking getoond voor vier aandoeningen, met genen op voorwaartse en achterwaartse strengen afzonderlijk weergegeven. X-as toont de basispositie ten opzichte van pA-splitsingsplaatsen, tot 100 nucleotiden stroomopwaarts en stroomafwaarts; y-as verwijst naar de gemiddelde dekking op overeenkomstige relatieve basisposities over alle pA-splitsingsplaatsen. Heatmaps bieden een alternatieve weergave, waarbij elke pA-splitsingssite wordt weergegeven als een afzonderlijke rij op de x-as, geordend op dekking. Intensiteit toont leesdekking (zie legenda). Klik hier om dit bestand te downloaden.

Aanvullende tabel 1. diffSplice output van de AS-analyse. Het eerste tabblad definieert de kolomnamen voor de oorspronkelijke uitvoer die wordt gepresenteerd in de tabbladen tweede (exon-niveau) en derde (genniveau). Klik hier om deze tabel te downloaden.

Aanvullende tabel 2. DEXSeq output van de AS analyse. Het eerste tabblad definieert de kolomnamen voor de oorspronkelijke uitvoer die wordt gepresenteerd in de tweede (exon-niveau) en derde (gen-niveau) tabbladen. Klik hier om deze tabel te downloaden.

Aanvullende tabel 3. rMATS-output van de AS-analyse. Het eerste tabblad definieert de kolomnamen voor het samenvattingsbestand (tabblad 2) en de JCEC-bestanden voor elke gebeurtenis (tabblad 3-7). Klik hier om deze tabel te downloaden.

Aanvullende tabel 4. diffSplice output van de APA-analyse. Het eerste tabblad definieert de kolomnamen voor de oorspronkelijke uitvoer die wordt gepresenteerd in de tweede tabbladen (pA-siteniveau) en derde (genniveau). Klik hier om deze tabel te downloaden.

Aanvullende tabel 5. DEXSeq output van de APA analyse. Het eerste tabblad definieert de kolomnamen voor de oorspronkelijke uitvoer die wordt gepresenteerd in de tweede tabbladen (pA-siteniveau) en derde (genniveau). Klik hier om deze tabel te downloaden.

Aanvullende tabel 6. Een samenvatting van het aantal significant veranderde exonen voor AS- en pA-sites voor APA. Klik hier om deze tabel te downloaden.

Aanvullende tabel 7. Een overzicht van tools en pakketten die worden gebruikt in de AS/APA-analyse. Klik hier om deze tabel te downloaden.

Discussion

In deze studie evalueerden we exon-gebaseerde en event-gebaseerde benaderingen om AS en APA te detecteren in bulk RNA-Seq en 3'-end sequencinggegevens. De exon-gebaseerde AS-benaderingen produceren zowel een lijst van differentieel tot expressie gebrachte exonen als een rangschikking op genniveau, geordend op basis van de statistische significantie van de totale differentiële splicingactiviteit op genniveau (tabellen 1-2, 4-5). Voor het diffSplice-pakket wordt het differentiële gebruik bepaald door gewogen lineaire modellen op exon-niveau aan te brengen om de differentiële log fold-change van een exon te schatten ten opzichte van de gemiddelde log fold-change van de andere exonen binnen hetzelfde gen (per exon FC genoemd). De statistische significantie op genniveau wordt berekend door individuele exon-level significantietests te combineren tot een gengewijze test volgens de Simes-methode¹⁰. Deze rangschikking op genniveau differentiële splicingactiviteit kan vervolgens worden gebruikt om een gensetverrijkingsanalyse (GSEA) uit te voeren van de belangrijkste betrokken routes¹⁰. DEXSeq gebruikt een vergelijkbare strategie, door een gegeneraliseerd lineair model te plaatsen om differentieel exongebruik te meten, hoewel verschillend in bepaalde stappen zoals filtering, normalisatie en dispersieschatting. Bij het vergelijken van de top 500 gerangschikte exonen met AS-activiteit en APA met behulp van DEXSeq en DiffSplice, vonden we een overlap van respectievelijk 310 exonen en 300 pA-sites, wat de concordantie van de twee exon-gebaseerde benaderingen aantoont, wat ook werd aangetoond in een eerdere studie ²⁹. Het wordt aanbevolen om een combinatie te gebruiken van zowel een exon-gebaseerde (ofwel DEXSeq of diffSplice) als een op gebeurtenissen gebaseerde benadering voor uitgebreide detectie en classificatie van AS. Voor APA kunnen gebruikers kiezen voor DEXSeq of diffSplice: van beide methoden is aangetoond dat ze goed presteren in een breed scala aan transcriptomics-experimenten²⁹.

Bij het voorbereiden van de RNA-seq-bibliotheek voor een AS-analyse is het belangrijk om een strengspecifiek bulk-RNA-seq-protocol⁸ te gebruiken, omdat een groot deel van de genen in gewervelde genomen elkaar overlappen op verschillende strengen en een niet-strengspecifiek protocol niet in staat is om deze overlappende regio's te onderscheiden, waardoor de uiteindelijke exondetectie wordt verstoord. Een andere overweging is leesdiepte, waarbij splicinganalyses diepere sequencing vereisen dan DGE, bijvoorbeeld 30-60 miljoen reads per sample, versus 5-25 miljoen reads per sample voor DGE (https://sapac.support.illumina.com/bulletins/2017/04/considerations-for-rna-seq-read-length-and-coverage-.html). Alle tools die in het protocol worden gedemonstreerd, ondersteunen zowel single-end als paired-end sequencing-gegevens. Als alleen bekende genannotaties worden gebruikt om junction reads te detecteren, kunnen single-ended kortere reads (≥ 50 bp) worden gebruikt, hoewel de novo detectie van nieuwe splice junctions profiteert van paired-end en langere (≥ 100bp) leest^30,31. De keuze van het RNA-extractieprotocol - polyA-selectie of rRNA-depletie - hangt af van de kwaliteit van RNA en de experimentele vraag - zie³¹ voor een discussie. Afhankelijk van de details van de bibliotheekconstructie, zijn wijzigingen nodig in de voorbeeldscripts die hier worden gegeven voor de parameters van leesuitlijning, functietelling en rMATS. Bij het berekenen van de initiële exon level read counts met behulp van featureCounts, of vergelijkbare methoden, moet ervoor worden gezorgd dat de functie-opties correct worden geconfigureerd voor tellingen en strandedness: in featureCounts stellen we het argument "strandSpecific" op de juiste manier in voor het strand-specifieke RNA-seq-protocol dat wordt gebruikt; en voor kwantificering op exon-niveau wordt verwacht dat een read zal worden toegewezen over aangrenzende exonen, en dus stellen we de parameter allowMultiOverlap in op TRUE. Voor APA zijn er verschillende 3'-eindvolgprotocollen⁶ die variëren in de precieze locatie van pieken ten opzichte van de pA-site. Voor onze voorbeeldgegevens bepalen we dat de piek 60 bp stroomopwaarts van de pA-site is, zoals weergegeven in figuur 5, en deze analyse zal moeten worden aangepast voor andere 3'-eindsequencingprotocollen.

In dit protocol beperken we de scope tot de bespreking van differentiële analyses op het niveau van individuele exonen, en splicing events bestaande uit aangrenzende exon-intron combinaties. We bespreken niet de klasse van analyses op basis van isoform de novo reconstructie zoals Cufflinks, Cuffdiff³², RSEM³³, Kallisto³⁴ die tot doel hebben de absolute en relatieve expressie van volledige alternatieve isovormen te detecteren en te kwantificeren. De exon- en event-based methoden zijn gevoeliger voor het detecteren van individuele splicing events³⁰ en bieden in veel gevallen alle informatie die nodig is voor verdere analyse, zonder dat er isoform-level kwantificering nodig is.

De nieuwste versie van de bronbestanden in dit protocol is beschikbaar op https://github.com/jiayuwen/AS_APA_JoVE

Materials

Name	Company	Catalog Number	Comments
Not relevent for computational study