Identifisering av alternativ spleising og polyadenylering i RNA-seq-data

Summary

Alternativ spleising (AS) og alternativ polyadenylering (APA) utvider mangfoldet av transkripsjonsisoformer og deres produkter. Her beskriver vi bioinformatiske protokoller for å analysere bulk RNA-seq og 3 'end sekvenseringsanalyser for å oppdage og visualisere AS og APA som varierer på tvers av eksperimentelle forhold.

Abstract

I tillegg til den typiske analysen av RNA-Seq for å måle differensialgenuttrykk (DGE) på tvers av eksperimentelle / biologiske forhold, kan RNA-seq-data også brukes til å utforske andre komplekse reguleringsmekanismer på eksonnivå. Alternativ spleising og polyadenylering spiller en avgjørende rolle i det funksjonelle mangfoldet av et gen ved å generere forskjellige isoformer for å regulere genuttrykk på posttranskripsjonsnivå, og begrensende analyser til hele gennivået kan gå glipp av dette viktige reguleringslaget. Her demonstrerer vi detaljerte trinnvise analyser for identifisering og visualisering av differensial ekson- og polyadenyleringsstedsbruk på tvers av forhold, ved bruk av Bioconductor og andre pakker og funksjoner, inkludert DEXSeq, diffSplice fra Limma-pakken og rMATS.

Introduction

RNA-seq har blitt mye brukt gjennom årene, typisk for å estimere differensial genuttrykk og genoppdagelse¹. I tillegg kan det også brukes til å estimere varierende eksonnivåbruk på grunn av genuttrykk for forskjellige isoformer, og dermed bidra til en bedre forståelse av genregulering på posttranskripsjonsnivå. Flertallet av eukaryote gener genererer forskjellige isoformer ved alternativ spleising (AS) for å øke mangfoldet av mRNA-uttrykk. AS-hendelser kan deles inn i forskjellige mønstre: hopping av komplette eksoner (SE) der en ("kassett") ekson fjernes helt ut av transkripsjonen sammen med dens flankerende introner; alternativ (donor) 5' valg av spleisested (A5SS) og alternativ 3' (akseptor) spleisestedsvalg (A3SS) når to eller flere skjøtesteder er til stede i hver ende av en ekson; oppbevaring av introner (RI) når et intron beholdes i det modne mRNA-transkriptet og gjensidig utelukkelse av eksonbruk (MXE) der bare ett av de to tilgjengelige eksonene kan beholdes om gangen ^2,3. Alternativ polyadenylering (APA) spiller også en viktig rolle i regulering av genuttrykk ved bruk av alternative poly (A) steder for å generere flere mRNA-isoformer fra et enkelt transkripsjon⁴. De fleste polyadenyleringssteder (pAs) er lokalisert i 3 'uoversatte regionen (3' UTRs), og genererer mRNA-isoformer med forskjellige 3 'UTR-lengder. Siden 3' UTR er det sentrale knutepunktet for å gjenkjenne regulatoriske elementer, kan forskjellige 3' UTR-lengder påvirke mRNA-lokalisering, stabilitet og oversettelse⁵. Det er en klasse med 3 'endesekvenseringsanalyser optimalisert for å oppdage APA som er forskjellige i detaljene i protokollen⁶. Rørledningen beskrevet her er designet for PolyA-seq, men kan tilpasses for andre protokoller som beskrevet.

I denne studien presenterer vi en pipeline av differensielle eksonanalysemetoder^7,8 (figur 1), som kan deles inn i to brede kategorier: eksonbasert (DEXSeq^9, diffSplice 10) og hendelsesbasert (replikat multivariat analyse av transkripsjonsspleising (rMATS)¹¹). De eksonbaserte metodene sammenligner foldeendringen på tvers av forholdene til individuelle eksoner, mot et mål på total genfoldendring for å kalle differensielt uttrykt eksonbruk, og fra det beregnes et gennivåmål for AS-aktivitet. Hendelsesbaserte metoder bruker ekson-intron-spennende kryssavlesninger for å oppdage og klassifisere spesifikke spleisehendelser som eksonhopping eller oppbevaring av introner, og skille disse AS-typene i utgangen³. Dermed gir disse metodene komplementære visninger for en fullstendig analyse av AS^12,13. Vi valgte DEXSeq (basert på DESeq2¹⁴ DGE-pakken) og diffSplice (basert på Limma¹⁰ DGE-pakken) for studien, da de er blant de mest brukte pakkene for differensiell spleisingsanalyse. rMATS ble valgt som en populær metode for hendelsesbasert analyse. En annen populær hendelsesbasert metode er MISO (Blanding av isoformer)¹. For TFO tilpasser vi den eksonbaserte tilnærmingen.

Figur 1. Analyse rørledning. Flytskjema for trinnene som brukes i analysen. Trinn inkluderer: innhenting av data, utførelse av kvalitetskontroller og lesejustering etterfulgt av telling av avlesninger ved hjelp av merknader for kjente eksoner, introner og pA-nettsteder, filtrering for å fjerne lave tellinger og normalisering. PolyA-seq-data ble analysert for alternative pA-lokaliteter med diffSplice/DEXSeq-metoder, bulk RNA-Seq ble analysert for alternativ spleising på eksonnivå med diffSplice/DEXseq-metoder, og AS-hendelser analysert med rMATS. Vennligst klikk her for å se en større versjon av denne figuren.

RNA-seq-dataene som ble brukt i denne undersøkelsen ble hentet fra Gene Expression Omnibus (GEO) (GSE138691)¹⁵. Vi brukte mus RNA-seq data fra denne studien med to tilstandsgrupper: villtype (WT) og Muskelblind-lignende type 1 knockout (Mbnl1 KO) med tre replikasjoner hver. For å demonstrere differensiell analyse av bruk av polyadenyleringsstedet, oppnådde vi museembryofibroblaster (MEFs) PolyA-seq-data (GEO Accession GSE60487)¹⁶. Dataene har fire tilstandsgrupper: Wild-type (WT), Muscleblind-lignende type1/type 2 dobbel knockout (Mbnl1/2 DKO), Mbnl 1/2 DKO med Mbnl3 knockdown (KD) og Mbnl1/2 DKO med Mbnl3 kontroll (Ctrl). Hver tilstandsgruppe består av to replikaer.

	GEO-tiltredelse	SRA Run nummer	Eksempel på navn	Betingelse	Replikere	Vev	Sekvensering	Lese lengde
RNA-Seq	GSM4116218	SRR10261601	Mbnl1KO_Thymus_1	Mbnl1 knockout	Rep 1	Thymus	Par-ende	100 bp
	GSM4116219	SRR10261602	Mbnl1KO_Thymus_2	Mbnl1 knockout	Rep 2	Thymus	Par-ende	100 bp
	GSM4116220	SRR10261603	Mbnl1KO_Thymus_3	Mbnl1 knockout	Rep 3	Thymus	Par-ende	100 bp
	GSM4116221	SRR10261604	WT_Thymus_1	Vill type	Rep 1	Thymus	Par-ende	100 bp
	GSM4116222	SRR10261605	WT_Thymus_2	Vill type	Rep 2	Thymus	Par-ende	100 bp
	GSM4116223	SRR10261606	WT_Thymus_3	Vill type	Rep 3	Thymus	Par-ende	100 bp
3P-Seq	GSM1480973	SRR1553129	WT_1	Vill type (WT)	Rep 1	Mus embryonale fibroblaster (MEFs)	Single-end	40 bp
	GSM1480974	SRR1553130	WT_2	Vill type (WT)	Rep 2	Mus embryonale fibroblaster (MEFs)	Single-end	40 bp
	GSM1480975	SRR1553131	DKO_1	Mbnl 1/2 dobbel knockout (DKO)	Rep 1	Mus embryonale fibroblaster (MEFs)	Single-end	40 bp
	GSM1480976	SRR1553132	DKO_2	Mbnl 1/2 dobbel knockout (DKO)	Rep 2	Mus embryonale fibroblaster (MEFs)	Single-end	40 bp
	GSM1480977	SRR1553133	DKOsiRNA_1	Mbnl 1/2 dobbel knockout med Mbnl 3 siRNA (KD)	Rep 1	Mus embryonale fibroblaster (MEFs)	Single-end	40 bp
	GSM1480978	SRR1553134	DKOsiRNA_2	Mbnl 1/2 dobbel knockout med Mbnl 3 siRNA (KD)	Rep 2	Mus embryonale fibroblaster (MEFs)	Single-end	36 bp
	GSM1480979	SRR1553135	DKONTsiRNA_1	Mbnl 1/2 dobbel knockout med ikke-målretting siRNA (Ctrl)	Rep 1	Mus embryonale fibroblaster (MEFs)	Single-end	40 bp
	GSM1480980	SRR1553136	DKONTsiRNA_2	Mbnl 1/2 dobbel knockout med ikke-målretting siRNA (Ctrl)	Rep 2	Mus embryonale fibroblaster (MEFs)	Single-end	40 bp

Tabell 1. Sammendrag av RNA-Seq og PolyA-seq datasett som brukes til analysen.

Protocol

1. Installasjon av verktøy og R-pakker som brukes i analysen

1. Conda er en populær og fleksibel pakkebehandling som muliggjør praktisk installasjon av pakker med deres avhengigheter på tvers av alle plattformer. Bruk 'Anaconda' (conda package manager) for å installere 'conda' som kan brukes til å installere verktøyene / pakkene som kreves for analysen.
2. Last ned 'Anaconda' i henhold til systemkravene fra https://www.anaconda.com/products/individual#Downloads og installer den ved å følge instruksjonene i grafisk installasjonsprogram. Installer alle nødvendige pakker ved hjelp av 'conda' ved å skrive følgende på Linux-kommandolinjen.
   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. For å laste ned alle R-pakkene som brukes i protokollen, skriv inn følgende kode i R-konsollen (startet på Linux-kommandolinjen ved å skrive 'R') eller Rstudio-konsollen.
   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)
}
   MERK: I denne beregningsprotokollen vil kommandoer bli gitt som enten R Notebook-filer (filer med utvidelse ". Rmd"), R-kodefiler (filer med utvidelse ". R"), eller Linux Bash shell scripts (filer med utvidelsen ".sh"). R Notebook (Rmd) filer skal åpnes i RStudio bruker File | Åpne Fil..., og individuelle kodebiter (som kan være R-kommandoer eller Bash-skallkommandoer) kjører deretter interaktivt ved å klikke på den grønne pilen øverst til høyre. R-kodefiler kan kjøres ved å åpne i RStudio, eller på Linux-kommandolinjen ved å prefacing med "Rscript", f.eks Rscript example.R. Shell-skript kjøres på Linux-kommandolinjen ved å prefacingere skriptet med "sh" -kommandoen f.eks.sh example.sh.

2. Alternativ spleising (AS) analyse ved bruk av RNA-seq

1. Nedlasting og forhåndsbehandling av data
   MERK: Kodebitene som er kommentert nedenfor, er tilgjengelige i tilleggskodefilen "AS_analysis_RNASeq.Rmd", for å følge de enkelte trinnene interaktivt, og er også gitt som en bash script som skal kjøres i batch på Linux kommandolinjen (sh downloading_data_preprocessing.sh).
   1. Nedlasting av rådata.
      1. Last ned rådata fra Sequence Read Archive (SRA) ved hjelp av kommandoen 'prefetch' fra SRA toolkit (v2.10.8)¹⁷. Gi SRA Accession ID-ene i rekkefølge i følgende kommando for å laste dem ned parallelt ved hjelp av GNU parallellverktøy¹⁸. For å laste ned SRA-filer med tiltredelses-ID-er fra SRR10261601 til SRR10261606 parallelt, bruk følgende på Linux-kommandolinjen.
         ​seq 10261601 10261606 | parallel prefetch SRR{}
      2. Pakk ut fastq-filene fra arkivet ved hjelp av 'fastq-dump' -funksjonen fra SRA toolkit. Bruk GNU parallelt og gi navnene på alle SRA-filer sammen.
         ​parallel -j 3 fastq-dump --gzip --skip-technical --read-filter pass --dumpbase --split-e --clip --origfmt {} :::
      3. Last ned referansegenomet og merknadene for Mouse (Genome assembly GRCm39) fra www.ensembl.org ved å bruke følgende på Linux-kommandolinjen.
         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. Forbehandling og kartlegging leser til genomsamlingen
      1. Kvalitetskontroll. Vurder kvaliteten på råavlesninger med FASTQC (FASTQ Quality Check v0.11.9)¹⁹. Opprett en utdatamappe og kjør fastqc parallelt på flere input fasta-filer. Dette trinnet genererer en kvalitetsrapport for hvert utvalg. Undersøk rapportene for å sikre at kvaliteten på lesingene er akseptabel før du gjør videre analyse. (Se brukerhåndboken for å forstå rapportene på https://www.bioinformatics.babraham.ac.uk/projects/fastqc/)
         mkdir fastqc_out
parallel "fastqc {} -o fastqc_out" ::: $RAW_DATA/*.fastq.gz
         MERK: Utfør om nødvendig adaptertrimming med 'cutadapt'20 eller 'trimmomatic'21 for å fjerne sekvensering i flankerende adaptere, som varierer basert på RNA-fragmentstørrelse og leselengde. I denne analysen hoppet vi over dette trinnet da brøkdelen av lesingene som ble berørt var minimal.
      2. Les justering. Det neste trinnet i forbehandlingen inkluderer kartlegging av avlesningene til referansegenomet. For det første, bygg indeksen for referansegenomet ved hjelp av 'genomeGenerate'-funksjonen til STAR²² , og juster deretter råavlesningene til referansen (alternativt er forhåndsbygde indekser tilgjengelige fra STAR-nettstedet og kan brukes direkte til justering). Kjør følgende kommandoer på Linux-kommandolinjen.
         #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
         MERK: STAR-justeringen vil generere og sortere BAM-filer (Binary Alignment Map) for hver prøve etter lesejustering. BAM-filer må sorteres før du går videre til ytterligere trinn.
2. Forbereder Exon-merknader.
   1. Kjør tilleggskodefilen "prepare_mm_exon_annotation. R" med den nedlastede merknaden i GTF (Gene transfer format) format for å forberede merknadene. For å kjøre, skriv inn følgende på Linux-kommandolinjen.
      Rscript prepare_mm_exon_annotation.R annotation.gtf
      MERK: GTF-filen inneholder flere eksonoppføringer for forskjellige isoformer. Denne filen brukes til å "kollapse" flere transkripsjons-IDer for hvert ekson. Det er et viktig skritt å definere eksontellingskasser.
3. Telling leser. Det neste trinnet er å telle antall leser kartlagt til forskjellige transkripsjoner / eksoner. Se Tilleggsfil: "AS_analysis_RNASeq.Rmd".
   1. Last inn nødvendige biblioteker:
      packages<- c("Rsubread","tidyverse", "magrittr", "EnhancedVolcano", "edgeR","openxlsx")
invisible(lapply(packages, library, character.only=TRUE))
   2. Last inn den behandlede merknadsfilen hentet fra forrige trinn (2.2).
      load("mm_exon_anno.RData")
   3. Les alle bam-filer hentet i trinn 2.2.2 som input for 'featureCounts' for å telle avlesninger. Les mappen som inneholder bam-filer ved først å liste opp hver fil fra katalogen som slutter med .bam. Bruk 'featureCounts' fra Rsubread-pakken som tar bam-filer og behandlet GTF-merknad (referanse) som input for å generere en matrise av tellinger knyttet til hver funksjon med rader som representerer eksoner (funksjoner) og kolonner som representerer prøver.
      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. Deretter utfører du ikke-spesifikk filtrering for å fjerne lavt uttrykte eksoner ("ikke-spesifikk" indikerer at den eksperimentelle tilstandsinformasjonen ikke brukes i filtreringen, for å unngå utvalgsskjevheter). Transformer dataene fra råskala til antall per million (cpm) ved hjelp av cpm-funksjonen fra 'edgeR'-pakke²³ , og hold eksoner med tellinger som er større enn en settbar terskel (for dette datasettet brukes én cpm) i minst tre utvalg. Fjern også gener med bare ett ekson.
      # 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, ]
      MERK: Kontroller de nødvendige parameterne for featureCounts når du bruker forskjellige data, for eksempel for single-end reads, set 'isPairedEnd = FALSE'. Se brukerveiledningen for RSubread for å velge alternativene for dataene dine, og se diskusjonsdelen nedenfor.
4. Differensiell spleising og eksonbruksanalyse. Vi beskriver to alternativer for dette trinnet: DEXSeq og DiffSplice. Begge kan brukes og gi lignende resultater. For konsistens, velg DEXSeq hvis du foretrekker en DESeq2-pakke for DGE og bruk DiffSplice for en Limma-basert DGE-analyse. Se Tilleggsfil: "AS_analysis_RNASeq.Rmd".
   1. Bruk av DEXSeq-pakke for differensiell eksonanalyse.
      1. Last inn bibliotek og opprett en eksempeltabell for å definere den eksperimentelle utformingen.
         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")))
         MERK: Radnavnene skal være i samsvar med bam-filnavnene som brukes av featureCounts for å telle lesingene. sampleTable består av detaljer om hver prøve som inkluderer: bibliotek-type og tilstand. Dette er nødvendig for å definere kontrastene eller testgruppen for å oppdage differensiell bruk.
      2. Klargjør eksoninformasjonsfilen. Eksoninformasjon i form av GRanges (genomiske områder) objekter (https://bioconductor.org/packages/release/bioc/html/GenomicRanges.html) kreves som en inngang for å lage DEXSeq-objektet i neste trinn. Match gen-ID-ene med lesetallene for å opprette exoninfo-objekt.
         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. Opprett DEXSeq-objekt ved hjelp av DEXSeqDataSet-funksjonen. DEXSeq-objektet samler sammen lesetellinger, informasjon om eksonfunksjoner og eksempelinformasjon. Bruk lesetallene som ble generert i trinn 3, og eksoninformasjonen fra forrige trinn til å opprette DEXSeq-objektet fra tellematrisen. SampleData-argumentet tar en datarammeinngang som definerer prøvene (og deres attributter: bibliotektype og tilstand), 'design' bruker sampleData til å generere en designmatrise for differensialtestingen ved hjelp av modellformelnotasjon. Merk at et signifikant interaksjonsbegrep, condition:exon, indikerer at brøkdelen av lesninger over et gen som faller på et bestemt ekson, avhenger av den eksperimentelle tilstanden, dvs. Se DEXSeq-dokumentasjon for en fullstendig beskrivelse av å sette modellformelen for mer komplekse eksperimentelle design. For funksjonsinformasjon kreves ekson-ID-er, tilsvarende gen og transkripsjoner.
         ​dxd<-DEXSeqDataSet(countData$counts,sampleData=sampleTable, design=~sample+exon+condition:exon,featureID=exoninfo$ExonID,groupID=exoninfo$GeneID,featureRanges=exoninfo, transcripts=transcripts_l)
      4. Normalisering og dispersjonsestimering. Deretter utfører du normalisering mellom prøver og estimerer variansen av dataene, på grunn av både Poisson-tellestøy fra den diskrete naturen til RNA-seq og biologisk variabilitet, ved hjelp av følgende kommandoer.
         ​dxd %<>% estimateSizeFactors %>% estimateDispersions %T>% plotDispEsts
      5. Test for differensiell bruk. Etter estimering av variasjon, test for differensiell eksonbruk for hvert gen og generer resultatene.
         dxd%<>%testForDEU%>%estimateExonFoldChanges(fitExpToVar=
"condition")#Estimate fold changes
dxr=DEXSeqResults(dxd)
      6. Visualiser skjøtehendelser for utvalgte gener ved hjelp av følgende kommando.
         plotDEXSeq(dxr,"Wnk1", displayTranscripts=TRUE, splicing=TRUE,legend
=TRUE,cex.axis=1.2,cex=1.3,lwd=2)
         Undersøk R Notebook-filen "AS_analysis_RNASeq.Rmd" for å generere flere tomter for gener av interesse og for å generere vulkanplott ved forskjellige terskler.
   2. Bruke diffSplice fra Limma for å identifisere differensiell spleising. Følg R Notebook-filen "AS_analysis_RNASeq.Rmd". Sørg for at trinn 2.1-2.3 er fulgt for å forberede inndatafiler før du fortsetter videre.
      1. Last inn biblioteker
         library(limma)
library(edgeR)
      2. Ikke-spesifikk filtrering. Trekk ut matrisen av lesetall oppnådd i 2,3. Opprett en liste over funksjoner ved hjelp av DGEList-funksjonen fra edgeR-pakken, der rader representerer gener og kolonner representerer prøver.
         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, "\\,")
         MERK: Som et ikke-spesifikt filtreringstrinn filtreres antall med cpm < 1 tommer x av n prøver, der x er det minste antallet replikasjoner i en hvilken som helst tilstand. n = 6 og x = 3 for disse eksempeldataene.
      3. Normaliser tellingene på tvers av prøver, med funksjonen 'calcNormFactors' fra 'edgeR'-pakken ved hjelp av Trimmet middelverdi av M-verdier (TMM-normaliseringsmetode)²⁴ Den vil beregne skaleringsfaktorer for å justere bibliotekstørrelser.
         ​dge<-calcNormFactors(dge)
      4. Bruk sampleTable som generert i trinn 2.4.1.1 og opprett utformingsmatrisen. Designmatrisen karakteriserer designet. Se Limma brukerhåndbok (https://www.bioconductor.org/packages/devel/bioc/vignettes/limma/inst/doc/usersguide.pdf) kapittel 8 og 9 for detaljer om designmatriser for mer avanserte eksperimentelle design.
         Treat<- factor(sampleTable$condition)
design<- model.matrix(~0+Treat)
colnames(design) <- levels(Treat)
      5. Tilpass en lineær modell per ekson. Kjør 'voom'-funksjonen til 'limma'-pakken for å behandle RNA-seq-data for å estimere varians og generere presisjonsvekter for å korrigere for Poisson-tellestøy, og transformere eksonnivåtellingene til log2-tellinger per million (logCPM). Kjør deretter lineær modellering ved hjelp av 'lmfit' -funksjonen for å tilpasse lineære modeller til uttrykksdataene for hvert ekson. Beregne empirisk Bayes-statistikk for tilpasset modell ved hjelp av 'eBayes'-funksjonen for å oppdage differensialeksonuttrykk. Deretter definerer du en kontrastmatrise for eksperimentelle sammenligninger av interesse. Bruk 'contrasts.fit' for å oppnå koeffisienter og standardfeil for hvert sammenligningspar.
         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. Differensiell spleisingsanalyse. Kjør 'diffSplice' på den monterte modellen for å teste forskjellene i eksonbruk av gener mellom wild-type og knockout og utforske de topprangerte resultatene ved hjelp av 'topSplice' -funksjon: test = "t" gir en rangering av AS exons, test = "simes" gir en rangering av gener.
         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. Visualisering. Plott resultatene med 'plotSplice' -funksjonen, noe som gir gen av interesse i geneidargumentet. Lagre de øverste resultatene sortert etter logg Fold endre til et objekt og generere en vulkan plott for å vise eksoner.
         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. Bruke rMATS
      1. Kontroller at den nyeste versjonen av rMATS v4.1.1 (også kjent som rMATS turbo på grunn av redusert behandlingstid og mindre krav til minne) er installert enten ved hjelp av conda eller github (https://github.com/Xinglab/rmats-turbo/releases/download/v4.1.1/rmats_turbo_v4_1_1.tar.gz) i arbeidskatalogen. Følg avsnitt 4.3 i "AS_analysis_RNASeq.Rmd".
      2. Gå til mappen som inneholder bam-filer oppnådd etter kartlegging og klargjør tekstfiler, som kreves av rMATS, for de to forholdene ved å kopiere navnet på bam-filer (sammen med banen) atskilt med ','. Følgende kommandoer skal kjøres på Linux-kommandolinjen:
         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. Kjør rmats.py med de to inngangsfilene som ble generert i forrige trinn, sammen med GTF-filen oppnådd i 2.1.1.3. Dette vil generere en utdatamappe 'rmats_out' som inneholder tekstfiler som beskriver statistikk (p-verdier og inkluderingsnivåer) for hver skjøtehendelse separat.
         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
         MERK: Referansemerknaden i form av en GTF-fil er også nødvendig. Kontroller parametrene hvis dataene er single-end, og endre alternativet -t tilsvarende.
      4. Utforske rMATS-resultater. Bruk Bioconductor-pakken '^maser'25 for å utforske rMATS-resultatene. Last tekstfilene Junction and exon counts (JCEC) inn i 'maser'-objektet og filtrer resultatet basert på dekning ved å inkludere minst fem gjennomsnittlige avlesninger per skjøtehendelse.
         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. Visualisering av rMATS-resultater. Velg de signifikante skjøtehendelsene ved False Discovery Rate (FDR) 10 % og minimum 10 % endring i Percent Spliced In (deltaPSI) ved hjelp av 'topEvents'-funksjonen fra 'maser'-pakken. Deretter sjekker du genhendelsene for individuelle gener av interesse (prøvegen-Wnk1) og plotter PSI-verdier for hver spleisingshendelse av det genet. Generer et vulkanplott ved å angi hendelsestypen.
         #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. Generer Sashimi-plott for spleisehendelsesresultatet oppnådd med rMATS i form av tekstfiler ved hjelp av 'rmats2shahimiplot' -pakken. Kjør python-skriptet på Linux-kommandolinjen.
         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
         MERK: Denne prosessen kan være tidkrevende, da den vil generere Sashimi-plottet for alle resultatene i hendelsesfilen. Velg toppresultatene (gennavn og eksoner) som vist av topEvents-funksjonen fra 'maser' og visualiser det tilsvarende Sashimi-plottet.

3. Alternativ polyadenylering (APA) analyse ved bruk av 3 'end sekvensering

1. Nedlasting og forhåndsbehandling av data
   MERK: Se den supplerende R Notebook-filen "APA_analysis_3PSeq_notebook. Rmd" for de komplette kommandoene for nedlasting av data og forhåndsbehandlingstrinn, eller kjør den supplerende bash-filen "APA_data_downloading_preprocessing.sh" på Linux-kommandolinjen.
   1. Last ned data fra SRA med tiltredelses-ID-ene (1553129 til 1553136).
   2. Trim adaptere og omvendt komplement for å oppnå sansestrengsekvensen.
      MERK: Dette trinnet er spesifikt for PolyA-seq-analysen som brukes.
   3. Kartet leser til mus genom montering ved hjelp av bowtie aligner²⁶.
2. Klargjøre merknader for pA-nettsteder.
   MERK: Behandlingen av pA-nettstedets merknadsfil utføres først ved hjelp av supplerende R Notebook-fil "APA_analysis_3PSeq_notebook. Rmd" (2.1 - 2.6), og deretter bruke bash-filen "APA_annotation_preparation.sh".
   1. Last ned merknader fra pA-nettsteder fra PolyASite 2.0-databasen⁶.
   2. Velg pA-nettstedsmerknader for å beholde 3'-uoversatte region (UTR) pA-nettsteder, som er kommentert som Terminal Exon (TE) eller 1000 nt nedstrøms for en kommentert terminalekson (DS) for nedstrømsanalyse.
   3. Få pA-områdetopper. Forankre på hvert pA-spaltningssted, og visualiser gjennomsnittlig lesedekning ved hjelp av bedtools og deeptools^27,28. Resultatene viste at toppene i de kartlagte målingene hovedsakelig var spredt innenfor ~60 bp oppstrøms spaltningsstedene (figur 5 og supplerende figur 5). Derfor ble koordinatene til pA-nettsteder utvidet fra merknadsfilen til 60 bp oppstrøms for spaltningsstedene. Avhengig av den spesifikke 3'-endesekvenseringsprotokollen som brukes, må dette trinnet optimaliseres for andre analyser enn PolyA-seq.
3. Telling leser
   1. Klargjør merknadsfilen for pA-nettsteder.
      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. Bruk 'featureCounts' for å skaffe rå teller. Lagre telletabellen som RData-filen "APA_countData.Rdata" for APA-analyse ved hjelp av forskjellige verktøy.
      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")
      MERK: Vær bevisst på å endre noen av parameterne som er oppført i funksjonen 'featureCounts'. Endre parameteren 'strandSpecific' for å sikre at den er i samsvar med sekvenseringsretningen til 3'-endesekvenseringsanalysen som brukes (empirisk vil visualisering av dataene i en genomleser over gener på pluss- og minustråder avklare dette).
   3. Bruk ikke-spesifikk filtrering av countData. Filtrering kan forbedre den statistiske robustheten betydelig i differensielle pA-brukstester. Først fjernet vi disse genene med bare ett pA-nettsted, der differensielt pA-bruk av nettstedet ikke kan defineres. For det andre bruker vi ikke-spesifikk filtrering basert på dekning: antall filtreres av cpm mindre enn 1 i x av n prøver, hvor x er det minste antall replikasjoner i alle forhold. N = 8 og x = 2 for disse eksempeldataene.
      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. Differensiell analyse av bruk av polyadenyleringssted ved bruk av DEXSeq og diffSplice-rørledninger.
   1. Bruke DEXSeq-pakken
      MERK: Siden en kontrastmatrise ikke kan defineres for DEXSeq-rørledningen, må den differensielle APA-analysen av hver to eksperimentelle betingelser utføres separat. Den differensielle APA-analysen av tilstanden WT og tilstanden DKO utføres som et eksempel for å forklare prosedyren. Se tilleggsfilen "APA_analysis_3PSeq_notebook. Rmd" for trinnvis arbeidsflyt i denne delen og differensiell APA-analyse av andre kontraster.
      1. Last inn bibliotek og opprett en eksempeltabell for å definere den eksperimentelle utformingen.
         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. Klargjør informasjonsfilen for pA-nettsteder ved hjelp av Bioconductor-pakken 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. Bruk lesetallene som ble generert i trinn 3.3, og informasjonen om pA-området fra forrige trinn til å opprette DEXSeq-objektet.
         # 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. Definer kontrastparet ved å definere nivåene av betingelser i DEXSeq-objektet.
         dxd1$condition<- factor(dxd1$condition, levels= c("WT", "DKO"))
# The contrast pair will be "DKO - WT"
      5. Normalisering og dispersjonsestimering. I likhet med RNA-seq-data, utfører 3' endesekvenseringsdata normalisering mellom prøver (kolonnevis median av forhold for hver prøve) ved hjelp av 'estimateSizeFactors' -funksjonen, og estimerer variasjonen av dataene ved hjelp av 'estimateDispersions' -funksjonen, og visualiser deretter spredningsestimeringsresultatet ved hjelp av 'plotDispEsts' -funksjonen.
         ​dxd1 %<>% estimateSizeFactors %>% estimateDispersions %T>% plotDispEsts
      6. Differensiell pA-brukstesting for hvert gen ved hjelp av funksjonen 'testForDEU', estimerer deretter foldeendringen av pA-nettstedbruk ved hjelp av funksjonen 'estimateExonFoldChanges'. Kontroller resultatene med funksjonen 'DEXSeqResults' og sett 'FDR < 10%' som kriterium for signifikant differensielle pA-nettsteder.
         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. Visualisering av differensialet pA-bruksresultater ved bruk av differensielle APA-plott generert av funksjonen 'plotDEXSeq' og vulkanplott av funksjonen '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. Bruker diffSplice-pakke. Se den supplerende R Notebook-filen "APA_analysis_3PSeq_notebook. Rmd" for den trinnvise arbeidsflyten i denne delen.
      1. Definer interessekontraster for differensiell pA-bruksanalyse.
         MERK: Dette trinnet skal utføres etter konstruksjon og behandling av DGEList-objektet, som er inkludert i R Notebook-filen "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. Visualiser resultatet av interessekontraster (her "DKO - WT") ved hjelp av differensielle APA-plott med funksjonen 'plotSplice' og vulkanplott med funksjonen 'EnhancedVolcano'. Se R Notebook-filen "APA_analysis_3PSeq_notebook. Rmd" 4.2.7 - 4.2.9 for visualisering av andre kontrastpar.
         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

Etter å ha kjørt den trinnvise arbeidsflyten ovenfor, er AS- og TFO-analyseutgangene og representative resultater i form av tabeller og dataplott, generert som følger.

SOM:
Hovedutgangen av AS-analysen (supplerende tabell 1 for diffSplice; Tabell 2 for DEXSeq) er en liste over eksoner som viser differensiell bruk på tvers av forhold, og en liste over gener som viser signifikant samlet skjøteaktivitet av en eller flere av dens bestanddeler, rangert etter statistisk signifikans. Supplerende tabell 1, tab 2 viser signifikante eksoner, med kolonner som viser differensial FC for ekson versus hvile, per ekson ujustert p-verdi og justert p-verdi (Benjamini-Hockberg-korreksjon). Terskling på de justerte p-verdiene vil gi et sett med eksoner med definert FDR. Supplerende tabell 1, tab 3 viser en rangert liste over gener som viser betydningen av samlet skjøteaktivitet, med en kolonne som viser gennivåjustert p-verdi beregnet ved hjelp av Simes-metoden. Tilsvarende data er vist i tabell 2 for DEXSeq. Supplerende figur 1 og supplerende figur 2 viser differensiell skjøtemønster i Mbnl1-, Tcf7- og Lef1-gener som er eksperimentelt validert i den publiserte artikkelen presentert med dataene¹⁵. Forfatterne har vist eksperimentell validering av fem gener - Mbnl1, Mbnl2, Lef1, Tcf7 og Ncor2. Vår tilnærming oppdaget differensielle spleisingspattens i alle disse genene. Her presenterer vi FDR-nivåene for hvert gen ved bruk av henholdsvis DEXSeq, diffSplice og rMATS som oppnådd i tilleggstabellene 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) og Ncor2 (9.2E-11, 1.6E-06, 0).

Figur 2 viser en sammenligning mellom utganger hentet fra tre ulike verktøy og illustrerer alternative skjøtemønstre i Wnk1-genet. Vulkanplottene er vist i figur 2A (diffSplice) og figur 2B (DEXSeq). Ytterligere tre høyt rangerte gener er vist i supplerende figur 1 (diffSplice) og supplerende figur 2 (DEXSeq). Y-aksen viser statistisk signifikans (-log10 P-verdier) og x-aksen viser effektstørrelse (foldeendring). Gener plassert øverst til venstre eller høyre kvadranter indikerer betydelig FC og sterke statistiske bevis på ekte forskjeller.

Figur 2. Sammenligning av alternative skjøteresultater oppnådd fra diffSplice, DEXSeq og rMATS. |
(A) Vulkanplott (venstre) av RNA-Seq fra Limma diffSplice-analyse: X-aksen viser log exon fold endring; y-aksen viser -log₁₀ p-verdi. Hvert punkt tilsvarer en ekson. Horisontal stiplet linje ved p-verdi = 1E-5; vertikale stiplede linjer ved todelt endring (FC). Røde eksoner viser betydelig FC og statistisk signifikans. Differensiell spleisingsplott (til høyre): Spleisingsmønstre vises for et eksempelgen Wnk1 der x-aksen viser ekson-id per transkripsjon; y-aksen viser ekson relativ loggfoldendring (forskjellen mellom eksons logFC og den samlede logFC for alle andre eksoner). Eksoner uthevet i rødt viser statistisk signifikant differensialuttrykk (FDR < 0,1).
(B) Vulkanplott (venstre) og differensiell eksonbruk (høyre) av RNA-Seq hentet fra DEXSeq-analyse. Wnk1-genet viser signifikant differensiell bruk av eksoner mellom WT og Mbnl1 knock-out uthevet i rosa, som tilsvarer de samme differensialeksonene i (A).
(C) Vulkanplott (venstre) og Sashimi-plott (høyre) for Wnk1 hentet fra rMATS-analyse. Vulkanplott som viser signifikant hoppet over (kassett) ekson (SE) hendelse i villtype sammenlignet med knockout ved 10% FDR med endring i prosent spleiset i (PSI eller ΔΨ) verdier > 0,1. X-aksen viser endring i PSI-verdier på tvers av betingelser, og y-aksen viser logg-P-verdi. Sashimi-plottet viser en hoppet eksonhendelse i Wnk1-genet , tilsvarende en signifikant differensialekson i (A) og (B). Hver rad representerer en RNA-Seq-prøve: tre replikasjoner av villtype og Mbnl1 knock-out. Høyden viser lesedekning i RPKM, og forbindelsesbuene viser kryssavlesninger over eksoner. Annoterte genmodellalternative isoformer er vist nederst i plottet. Det nederste panelet på C illustrerer kryssavlesninger som brukes til å beregne PSI-statistikken.
(D) Venn-diagram som sammenligner antall signifikante differensialeksoner oppnådd ved de forskjellige metodene. Vennligst klikk her for å se en større versjon av denne figuren.

Figur 2 A (høyre panel) viser en diagrammatisk visning av eksonforskjeller for et av de høyest rangerte genene, og viser logFC på y-aksen og eksontallet på x-aksen. Dette eksemplet viser en kassettekson som varierer mellom tilstandene for genet Wnk1. Det differensielle eksonbruksplottet fra DEXSeq viser tegn på differensiell spleising på fem eksonsteder nær Wnk1.6.45. De uthevede eksonene i rosa vil sannsynligvis bli spleiset ut i Mbnl1 KO-prøver sammenlignet med WT. Disse eksonene er komplementære til resultatene oppnådd ved diffSplice som viser et lignende mønster ved den spesifikke genomiske posisjonen. Flere eksempler er vist i supplerende figur 1 og supplerende figur 2. En mer detaljert visning for å bekrefte interessante resultater kan gis ved å sammenligne dekning (wiggle) spor i RPM (Leser per million) enheter i UCSC (University of Santa Cruz) eller IGV (Integrative Genomics Viewer) genom nettlesere (ikke vist), sammen med visuell korrelasjon med andre spor av interesse, for eksempel kjente genmodeller, bevaring og andre genom-wide analyser.

rMATS-utgangstabellen viser signifikante alternative skjøtehendelser kategorisert etter type (tilleggstabell 3). Figur 2C viser et vulkanplott av gener som er alternative skjøtet, med effektstørrelsen målt ved differensiell "prosent spleiset inn" (PSI eller ΔΨ) statistikk på¹¹.

PSI refererer til prosentandelen av avlesninger som er i samsvar med inkluderingen av en kassettekson (dvs. leser kartlegging til selve kassetteksonen eller krysslesingen overlapper eksonet) sammenlignet med avlesninger i samsvar med eksonekskludering, dvs. krysslesninger over tilstøtende oppstrøms og nedstrøms eksoner (Det nederste panelet i figur 2C). Det høyre panelet i figur 2C viser sashimi-plottet av Wnk1-genet med differensiell spleisingshendelse overlagt på dekningsspor for genet, med en hoppet ekson i Mbnl1 KO. Buer som slutter seg til eksoner viser antall kryssavlesninger (leser krysser en spleiset ut intron). Ulike faner i tilleggstabell 3 viser signifikante avlesninger av hver type hendelse som spenner over kryss med eksongrenser (krysstellinger og eksontellinger (JCEC)). Figur 2D sammenligner de signifikante differensielt skjøtede eksonene som oppdages av de tre verktøyene.

Figur 3. Alternative skjøtehendelser oppnådd ved rMATS-analyse. A) Typer AS-hendelser. Dette tallet er tilpasset fra rMATS-dokumentasjon¹¹ som forklarer skjøtehendelsen med konstitutive og alternativt skjøtede eksoner. Merket med nummeret på hver hendelse på FDR 10%. B) Sashimi-plottet av Add3-genet som viser hoppet over ekson (SE). C) Sashimi-plottet av Baz2b-genet som viser alternativ 5' spleisested (A5SS). D) Sashimi-plottet av Lsm14b-genet som viser alternativ 3' spleisested (A3SS). E) Sashimi-plottet av Mta1-genet som viser gjensidig utelukkende eksoner (MXE). F) Sashimi-plottet av Arpp21-genet som viser beholdt intron (RI). Røde rader representerer tre replikaer av villtype og oransje rader representerer Mbnl1 knock-out-replikaer. X-aksen tilsvarer genomiske koordinater og strenginformasjon, y-aksen viser dekning i RPKM. Vennligst klikk her for å se en større versjon av denne figuren.

Figur 3 illustrerer typer spleisingshendelser SE, A5SS, A3SS, MXE og RI ved hjelp av Sashimi-plott av de viktigste genene til disse hendelsene. Ved sammenligning av de tre replikasjonene av både WT og Mbnl1_KO, ble totalt 1272 SE-hendelser, 130 A5SS, 116 A3SS, 215 MXE og 313 RI-hendelser detektert ved FDR 10%. Sashimi-plottet illustrerer typen hendelse ved å bruke toppgener som eksempel.

TFO:
Resultatet fra TFO-analysen ligner på eksonnivå AS-analysen. En tabell over toppgener rangert etter differensiell APA-aktivitet i 3'UTR er gitt (supplerende tabell 4 og supplerende tabell 5). Figur 4A viser vulkanplottene av gener ved differensiell APA-aktivitet i 3'UTRs generert ved hjelp av diffSplice og DEXSeq separat. Figur 4B viser Venn-plottet som sammenligner de signifikant differensielle resultatene fra pA-bruken av nettstedet som er oppnådd fra forskjellige rørledninger. Figur 4C og 4D viser den diagrammatiske representasjonen av differensiell pA-stedsbruk i 3'UTR av genet Fosl2 (figur 4C) og Papola (figur 4D) generert ved hjelp av både diffSplice og DEXSeq, som er eksperimentelt validert for å vise signifikant distal til proksimal skift (Fosl2) og signifikant proksimalt til distalt skift (Papola) av pA-bruk i DKO, henholdsvis. Flere eksempler er vist i supplerende figur 3 og supplerende figur 4.

Figur 4. Alternative polyadenyleringsplott med diffSplice og DEXSeq. A) Vulkanplott av PolyA-seq-data generert ved hjelp av diffSplice og DEXSeq. X-aksen viser logg pA nettsted fold endring; y-aksen viser -log₁₀ p-verdi. Hvert punkt tilsvarer et pA-nettsted. Horisontal stiplet linje ved p-verdi = 1E-5; vertikale stiplede linjer ved 2 ganger FC. Røde eksoner viser betydelig FC og statistisk signifikans. B) Venn-plott som sammenligner de signifikante differensielle pA-bruksresultatene som er oppnådd fra forskjellige rørledninger. C-D) Differensielle APA-plott generert ved hjelp av diffSplice og DEXSeq som viser proksimale, interne og distale pA-steder for Fosl2 - og Papola-genet . Tallene genereres av samme funksjon som figur 2 (B), men med pA-steder som erstatter eksoner. Vennligst klikk her for å se en større versjon av denne figuren.

Figur 5 er et diagnostisk plott for å bekrefte forventet lesefordeling rundt kommenterte pA-spaltningssteder for PolyA-seq-analysen som brukes. Den viser gjennomsnittlig dekning i flankerende regioner forankret på kjente pA-spaltningssteder på genomnivå. I dette tilfellet visualiseres den forventede bunken av avlesninger oppstrøms nettstedene. Lesefordelingene forankret på pA-lokaliteter for alle PolyA-seq-prøver er vist i supplerende figur 5.

Figur 5. Gjennomsnittlig dekningsplott rundt pA-spaltningssteder. Spaltningsstedet for PolyA-seq-data vises med vertikal stiplet linje. X-aksen viser baseposisjon i forhold til pA-spaltningssteder, opptil 100 nukleotider oppstrøms og nedstrøms; y-aksen viser gjennomsnittlig lesedekning over alle pA-spaltningssteder, normalisert etter bibliotekstørrelse i CPM. Vennligst klikk her for å se en større versjon av denne figuren.

De differensielle APA-resultatene av forskjellige kontraster generert av samme rørledning kan sammenlignes og verifiseres ved å visualisere lesedekningen til representative signifikant differensielle pA-steder i genomleseren. Figur 6A er Venn-plottet som sammenligner den signifikant differensielle pA-bruken av forskjellige kontraster ervervet fra diffSplice. Figur 6B-D er IGV-øyeblikksbildene av lesedekningen på pA-steder for forskjellige gener, som viser mønstrene som stemmer overens med de som ble oppdaget i APA-analysen ved hjelp av diffSplice. Figur 6B validerer den signifikante proksimale til distale forskyvningen av pA-stedsbruk for gen Paip2, som er unikt påvist i kontrasten DKO vs WT, men ikke i de to andre kontrastene KD vs WT, og Ctr vs WT. Figur 6C validerer den signifikante distale til proksimale forskyvningen av pA-stedsbruk for gen Ccl2 oppdaget unikt i kontrasten KD vs WT, mens figur 6D validerer den signifikante differensielle pA-bruken av alle kontrastene for genet Cacna2d1.

Figur 6. Kontrastsammenligning og verifisering av diffSplice-resultater. A) Venn-diagram som sammenligner signifikante differensielle pA-bruksresultater av forskjellige kontraster oppnådd fra diffSplice. B-D) IGV øyeblikksbilde visualisere pA topper dekning av gener Paip2, Ccl2 og Cacna2d1 på tvers av forhold. Hvert spor representerer lesedekningen i en bestemt tilstand. Vennligst klikk her for å se en større versjon av denne figuren.

Supplerende figur 1. RNA-Seq analyse av differensial spleising med Limma diffSplice. (A) Vulkanplott av RNA-Seq fra Limma diffSplice analyse: X-aksen viser log exon fold endring; y-aksen viser -log₁₀ p-verdi. Hvert punkt tilsvarer en ekson. Horisontal stiplet linje ved p-verdi = 1E-5; vertikale stiplede linjer ved todelt endring (FC). Røde eksoner viser betydelig FC og statistisk signifikans. (B-D) Differensielle skjøteplott: Spleisingsmønstre vises for eksempel gener Mbnl1, Tcf7 og Lef1, henholdsvis hvor x-aksen viser ekson-id per transkripsjon; y-aksen viser ekson relativ loggfoldendring (forskjellen mellom eksons logFC og den samlede logFC for alle andre eksoner). Eksoner uthevet i rødt viser statistisk signifikant differensialuttrykk (FDR < 0,1). Klikk her for å laste ned denne filen.

Supplerende figur 2. RNA-Seq-analyse av differensiell eksonbruk med DEXSeq. (A) Vulkan tomt. (B-D) Differensiell eksonbruk % av RNA-Seq oppnådd fra DEXSeq-analyse. Genene Mbnl1, Tcf7 og Lef1 viser signifikant differensiell bruk av eksoner mellom WT og Mbnl1 knock-out uthevet i rosa, som tilsvarer de samme differensialeksonene i tilleggsfigur 1. Klikk her for å laste ned denne filen.  

Supplerende figur 3. Alternative polyadenyleringsplott ved diffSplice. A) Vulkanplott av PolyA-seq-data generert ved hjelp av diffSplice i tre kontrastpar hentet fra musen PolyA-seq-data, inkludert dobbel knockout (DKO) vs wild-type (WT), knock-down (KD) vs WT, og kontroll (Ctrl) vs WT. X-akse viser logg pA nettsted fold endring; y-aksen viser -log₁₀ p-verdi. Hvert punkt tilsvarer et pA-nettsted. Horisontal stiplet linje ved p-verdi = 1E-5; vertikale stiplede linjer ved 2 ganger FC. Røde pA-nettsteder viser betydelig FC og statistisk signifikans. B) Differensielle APA-plott generert ved hjelp av diffSplice som viser de proksimale, interne og distale pA-stedene for de høyt rangerte genene S100a7a, Tpm1 og Smc6. Klikk her for å laste ned denne filen.  

Supplerende figur 4. Differensiell pA-bruksanalyse ved hjelp av DEXSeq-pipeline. A) Vulkanplott av PolyA-seq-data generert ved hjelp av DEXSeq i tre par hentet fra musen PolyA-seq-data, inkludert dobbel knockout (DKO) vs wild-type (WT), knock-down (KD) vs WT, og kontroll (Ctrl) vs WT. X-akse viser logg pA nettsted fold endring; y-aksen viser -log₁₀ p-verdi. Hvert punkt tilsvarer et pA-nettsted. Horisontal stiplet linje ved p-verdi = 1E-5; vertikale stiplede linjer ved 2 ganger FC. Røde pA-nettsteder viser betydelig FC og statistisk signifikans. B) Differensielle APA-plott generert ved hjelp av DEXSeq som viser de proksimale, interne og distale pA-stedene for de høyt rangerte genene S100a7a, Tpm1 og Smc6. Klikk her for å laste ned denne filen.  

Supplerende figur 5. Gjennomsnittlig dekningsplott og varmekart rundt pA-spaltningssteder.  Dekning vist for fire forhold, med gener på fremover og bakover tråder vist separat. X-aksen viser baseposisjon i forhold til pA-spaltningssteder, opptil 100 nukleotider oppstrøms og nedstrøms; y-akse refererer til gjennomsnittlig dekning ved tilsvarende relative basisposisjoner på tvers av alle pA-spaltningssteder. Varmekart gir en alternativ visning, med hvert pA-spaltningssted vist som en egen rad på x-aksen, sortert etter dekning. Intensitet viser lesedekning (se forklaring). Klikk her for å laste ned denne filen.

Supplerende tabell 1. diffSplice utgang av AS analysen. Første fane definerer kolonnenavnene for de opprinnelige utdataene som presenteres i andre (ekson-nivå) og tredje (gennivå) faner. Klikk her for å laste ned denne tabellen.

Supplerende tabell 2. DEXSeq-utdata fra AS-analysen. Første fane definerer kolonnenavnene for de opprinnelige utdataene som presenteres i andre (eksonnivå) og tredje (gennivå) faner. Klikk her for å laste ned denne tabellen.

Supplerende tabell 3. rMATS-utdata fra AS-analysen. Den første kategorien definerer kolonnenavnene for sammendragsfilen (tab 2) og JCEC-filene for hver hendelse (tab 3–7). Klikk her for å laste ned denne tabellen.

Supplerende tabell 4. diffSplice utgang av TFO-analysen. Første fane definerer kolonnenavnene for den opprinnelige utdataene presentert i andre (pA site-level) og tredje (gen-nivå) faner. Klikk her for å laste ned denne tabellen.

Supplerende tabell 5. DEXSeq-utgangen fra TFO-analysen. Første fane definerer kolonnenavnene for den opprinnelige utdataene presentert i andre (pA site-level) og tredje (gen-nivå) faner. Klikk her for å laste ned denne tabellen.

Supplerende tabell 6. Et sammendrag av antall signifikant endrede eksoner for AS- og pA-nettsteder for TFO. Klikk her for å laste ned denne tabellen.

Supplerende tabell 7. Et sammendrag av verktøy og pakker som brukes i AS/APA-analysen. Klikk her for å laste ned denne tabellen.

Discussion

I denne studien evaluerte vi eksonbaserte og hendelsesbaserte tilnærminger for å oppdage AS og APA i bulk RNA-Seq og 3 'end sekvenseringsdata. De eksonbaserte AS-tilnærmingene produserer både en liste over differensielt uttrykte eksoner og en gennivårangering sortert etter statistisk signifikans av samlet gennivådifferensiell spleisingsaktivitet (tabell 1-2, 4-5). For diffSplice-pakken bestemmes differensiell bruk ved å tilpasse vektede lineære modeller på eksonnivå for å estimere differensialloggfoldingsendringen av et ekson mot gjennomsnittlig loggfoldeendring av de andre eksonene i samme gen (kalt per exon FC). Den statistiske signifikansen på gennivå beregnes ved å kombinere individuelle eksonnivå-signifikanstester til en genvis test ved Simes-metoden¹⁰. Denne rangeringen etter differensiell spleisingsaktivitet på gennivå kan deretter brukes til å utføre en gensettberikelsesanalyse (GSEA) av viktige veier involvert¹⁰. DEXSeq bruker en lignende strategi ved å tilpasse en generalisert lineær modell for å måle differensiell eksonbruk, men forskjellig i visse trinn som filtrering, normalisering og dispersjonsestimering. Ved sammenligning av de 500 beste rangerte eksonene som viser AS-aktivitet og APA ved bruk av DEXSeq og DiffSplice, fant vi en overlapping på henholdsvis 310 eksoner og 300 pA-steder, noe som demonstrerte samsvaret mellom de to eksonbaserte tilnærmingene, som også ble demonstrert i en tidligere studie ²⁹. Det anbefales å bruke en kombinasjon av både en eksonbasert (enten DEXSeq eller diffSplice) og en hendelsesbasert tilnærming for omfattende deteksjon og klassifisering av AS. For APA kan brukerne velge enten DEXSeq eller diffSplice: begge metodene har vist seg å fungere godt over et bredt spekter av transkriptomikkeksperimenter²⁹.

Ved utarbeidelse av RNA-seq-biblioteket for en AS-analyse er det viktig å bruke en strengspesifikk bulk RNA-seq-protokoll⁸, da en stor del av gener i vertebrate genomer overlapper på forskjellige tråder, og en ikke-strengspesifikk protokoll ikke klarer å skille disse overlappende regionene, forvirrende endelig eksondeteksjon. En annen vurdering er lesedybde, med spleisingsanalyser som krever dypere sekvensering enn DGE, for eksempel 30-60 millioner avlesninger per prøve, mot 5-25 millioner avlesninger per prøve for DGE (https://sapac.support.illumina.com/bulletins/2017/04/considerations-for-rna-seq-read-length-and-coverage-.html). Alle verktøyene som er demonstrert i protokollen, støtter både single-end og paired-end sekvenseringsdata. Hvis bare kjente genmerknader brukes til å oppdage kryssavlesninger, kan single-ended kortere avlesninger (≥ 50 bp) brukes, selv om de novo-deteksjon av nye spleisekryss drar nytte av parrede ender og lengre (≥ 100bp) leser^30,31. Valget av RNA-ekstraksjonsprotokoll - enten polyA-seleksjon eller rRNA-uttømming - avhenger av kvaliteten på RNA og det eksperimentelle spørsmålet - se³¹ for en diskusjon. Avhengig av detaljene i bibliotekkonstruksjonen, vil det være nødvendig med endringer i eksempelskriptene som er gitt her for parametrene for lesejustering, funksjonstelling og rMATS. Ved beregning av de første lesetellingene på eksonnivå ved hjelp av featureCounts, eller lignende metoder, må det tas hensyn til å konfigurere funksjonsalternativene riktig for tellinger og strandedhet: i featureCounts setter vi "strandSpecific" -argumentet passende for den strengspesifikke RNA-seq-protokollen som brukes; og for kvantifisering på eksonnivå forventes det at en lesning vil kartlegge over tilstøtende eksoner, og derfor setter vi allowMultiOverlap-parameteren til TRUE. For APA er det forskjellige 3 'endesekvenseringsprotokoller⁶ som varierer i den nøyaktige plasseringen av topper i forhold til pA-området. For våre eksempeldata bestemmer vi at toppen er 60 bp oppstrøms for pA-nettstedet som vist i figur 5, og denne analysen må tilpasses for andre 3 'endesekvenseringsprotokoller.

I denne protokollen begrenser vi omfanget til diskusjon av differensialanalyser på nivå med individuelle eksoner, og spleisingshendelser bestående av tilstøtende ekson-intron-kombinasjoner. Vi diskuterer ikke klassen av analyser basert på isoform de novo rekonstruksjon som mansjettknapper, Cuffdiff³², RSEM 33, Kallisto³⁴ som tar sikte på å oppdage og kvantifisere det absolutte og relative uttrykket for hele alternative isoformer. Ekson- og hendelsesbaserte metoder er mer følsomme for å oppdage individuelle skjøtehendelser³⁰ og gir i mange tilfeller all informasjon som trengs for videre analyse, uten behov for kvantifisering på isoformnivå.

Den nyeste versjonen av kildefilene i denne protokollen er tilgjengelig på https://github.com/jiayuwen/AS_APA_JoVE

Materials

Name	Company	Catalog Number	Comments
Not relevent for computational study