Waiting
Login processing...

Trial ends in Request Full Access Tell Your Colleague About Jove
JoVE Journal Biology
Click here for the English version

Biology

Identifikation af alternativ splejsning og polyadenylering i RNA-seq-data

Published: June 24, 2021 doi: 10.3791/62636
Automatic Translation
1, 1, 2, 1
1Department of Genome Sciences, The John Curtin School of Medical Research, The Australian National University, 2Department of Biology, New York University

Summary

Alternativ splejsning (AS) og alternativ polyadenylering (APA) udvider mangfoldigheden af transkriptionsisoformer og deres produkter. Her beskriver vi bioinformatiske protokoller til analyse af bulk RNA-seq og 3 'slutsekventeringsassays for at detektere og visualisere AS og APA, der varierer på tværs af eksperimentelle forhold.

Abstract

Ud over den typiske analyse af RNA-Seq til måling af differentiel genekspression (DGE) på tværs af eksperimentelle / biologiske tilstande kan RNA-seq-data også bruges til at udforske andre komplekse reguleringsmekanismer på exon-niveau. Alternativ splejsning og polyadenylering spiller en afgørende rolle i et gens funktionelle mangfoldighed ved at generere forskellige isoformer til regulering af genekspression på post-transkriptionelt niveau, og begrænsning af analyser til hele genniveauet kan gå glip af dette vigtige regulatoriske lag. Her demonstrerer vi detaljerede trin for trin analyser til identifikation og visualisering af differentiel exon og polyadenylering site brug på tværs af betingelser ved hjælp af Bioconductor og andre pakker og funktioner, herunder DEXSeq, diffSplice fra Limma-pakken og rMATS.

Introduction

RNA-seq har været meget udbredt gennem årene typisk til estimering af differentiel genekspression og genopdagelse1. Derudover kan det også bruges til at estimere varierende exon-niveau brug på grund af gen, der udtrykker forskellige isoformer, hvilket bidrager til en bedre forståelse af genregulering på post-transkriptionelt niveau. Størstedelen af eukaryote gener genererer forskellige isoformer ved alternativ splejsning (AS) for at øge mangfoldigheden af mRNA-ekspression. AS-begivenheder kan opdeles i forskellige mønstre: spring over komplette exons (SE), hvor en ("kassette") exon fjernes helt ud af udskriften sammen med dens flankerende introns; alternativ (donor) 5' splejsningssted (A5SS) og alternativ 3' (acceptor) splejsningssted (A3SS), når to eller flere splejsningssteder er til stede i hver ende af en exon retention af introns (RI), når en intron bevares inden for det modne mRNA-transkript og gensidig udelukkelse af exonbrug (MXE), hvor kun en af de to tilgængelige exoner kan bevares ad gangen 2,3. Alternativ polyadenylering (APA) spiller også en vigtig rolle i reguleringen af genekspression ved hjælp af alternative poly (A) steder til at generere flere mRNA-isoformer fra et enkelt transkript4. De fleste polyadenyleringssteder (pA'er) er placeret i den 3' uoversatte region (3' UTR'er), der genererer mRNA-isoformer med forskellige 3' UTR-længder. Da 3' UTR er det centrale knudepunkt for genkendelse af regulatoriske elementer, kan forskellige 3' UTR-længder påvirke mRNA-lokalisering, stabilitet og translation5. Der er en klasse af 3 'ende sekventeringsassays optimeret til at detektere APA, der adskiller sig i detaljerne i protokollen6. Rørledningen beskrevet her er designet til PolyA-seq, men kan tilpasses til andre protokoller som beskrevet.

I denne undersøgelse præsenterer vi en pipeline af differentielle exonanalysemetoder7,8 (figur 1), som kan opdeles i to brede kategorier: exonbaseret (DEXSeq9, diffSplice 10) og begivenhedsbaseret (replikeret multivariat analyse af transkriptionsplejsning (rMATS)11). De exon-baserede metoder sammenligner foldændringen på tværs af betingelser for individuelle exoner mod et mål for den samlede genfoldændring for at kalde differentielt udtrykt exon-brug og ud fra det beregne et mål på genniveau for AS-aktivitet. Hændelsesbaserede metoder bruger exon-intron-spanning junction-læsninger til at registrere og klassificere specifikke splejsningshændelser såsom exon-spring eller tilbageholdelse af introns og skelne mellem disse AS-typer i output3. Disse metoder giver således supplerende synspunkter til en komplet analyse af AS12,13. Vi valgte DEXSeq (baseret på DESeq214 DGE-pakken) og diffSplice (baseret på Limma10 DGE-pakken) til undersøgelsen, da de er blandt de mest anvendte pakker til differentiel splejsningsanalyse. rMATS blev valgt som en populær metode til begivenhedsbaseret analyse. En anden populær begivenhedsbaseret metode er MISO (Blanding af isoformer)1. For APA tilpasser vi den exon-baserede tilgang.

Figure 1
Figur 1. Analyse pipeline. Rutediagram over de trin, der bruges i analysen. Trin inkluderer: indhentning af data, udførelse af kvalitetskontrol og læsejustering efterfulgt af optælling af læsninger ved hjælp af kommentarer til kendte exons, introns og pA-websteder, filtrering for at fjerne lave tællinger og normalisering. PolyA-seq-data blev analyseret for alternative pA-steder ved hjælp af diffSplice/DEXSeq-metoder, bulk RNA-Seq blev analyseret for alternativ splejsning på exonniveau med diffSplice/DEXseq-metoder og AS-hændelser analyseret med rMATS. Klik her for at se en større version af denne figur.

RNA-seq-dataene, der blev brugt i denne undersøgelse, blev erhvervet fra Gene Expression Omnibus (GEO) (GSE138691)15. Vi brugte mus RNA-seq data fra denne undersøgelse med to tilstandsgrupper: wild-type (WT) og muskelblind-lignende type 1 knockout (Mbnl1 KO) med tre replikater hver. For at demonstrere differentiel analyse af brugen af polyadenyleringsstedet opnåede vi museembryofibroblaster (MEF'er) PolyA-seq-data (GEO-tiltrædelse GSE60487)16. Dataene har fire tilstandsgrupper: Wild-type (WT), Muscleblind-lignende type1/type 2 dobbelt knockout (Mbnl1/2 DKO), Mbnl 1/2 DKO med Mbnl3 knockdown (KD) og Mbnl1/2 DKO med Mbnl3 kontrol (Ctrl). Hver tilstandsgruppe består af to replikater.

GEO tiltrædelse SRA Kør nummer Eksempel på navn Betingelse Kopiere Væv Sekventering Læs længde
RNA-Seq GSM4116218 SRR10261601 Mbnl1KO_Thymus_1 Mbnl1 knockout Rep 1 Thymus Parret ende 100 bp
GSM4116219 SRR10261602 Mbnl1KO_Thymus_2 Mbnl1 knockout Repræsentant 2 Thymus Parret ende 100 bp
GSM4116220 SRR10261603 Mbnl1KO_Thymus_3 Mbnl1 knockout Rep 3 Thymus Parret ende 100 bp
GSM4116221 SRR10261604 WT_Thymus_1 Vild type Rep 1 Thymus Parret ende 100 bp
GSM4116222 SRR10261605 WT_Thymus_2 Vild type Repræsentant 2 Thymus Parret ende 100 bp
GSM4116223 SRR10261606 WT_Thymus_3 Vild type Rep 3 Thymus Parret ende 100 bp
3P-Seq GSM1480973 SRR1553129 WT_1 Vildtype (WT) Rep 1 Museembryonale fibroblaster (MEF'er) Single-end 40 bp
GSM1480974 SRR1553130 WT_2 Vildtype (WT) Repræsentant 2 Museembryonale fibroblaster (MEF'er) Single-end 40 bp
GSM1480975 SRR1553131 DKO_1 Mbnl 1/2 dobbelt knockout (DKO) Rep 1 Museembryonale fibroblaster (MEF'er) Single-end 40 bp
GSM1480976 SRR1553132 DKO_2 Mbnl 1/2 dobbelt knockout (DKO) Repræsentant 2 Museembryonale fibroblaster (MEF'er) Single-end 40 bp
GSM1480977 SRR1553133 DKOsiRNA_1 Mbnl 1/2 dobbelt knockout med Mbnl 3 siRNA (KD) Rep 1 Museembryonale fibroblaster (MEF'er) Single-end 40 bp
GSM1480978 SRR1553134 DKOsiRNA_2 Mbnl 1/2 dobbelt knockout med Mbnl 3 siRNA (KD) Repræsentant 2 Museembryonale fibroblaster (MEF'er) Single-end 36 bp
GSM1480979 SRR1553135 DKONTsiRNA_1 Mbnl 1/2 dobbelt knockout med ikke-målrettet siRNA (Ctrl) Rep 1 Museembryonale fibroblaster (MEF'er) Single-end 40 bp
GSM1480980 SRR1553136 DKONTsiRNA_2 Mbnl 1/2 dobbelt knockout med ikke-målrettet siRNA (Ctrl) Repræsentant 2 Museembryonale fibroblaster (MEF'er) Single-end 40 bp

Tabel 1. Resumé af RNA-Seq og PolyA-seq datasæt, der anvendes til analysen.

Subscription Required. Please recommend JoVE to your librarian.

Protocol

1. Installation af værktøjer og R-pakker, der anvendes i analysen

  1. Conda er en populær og fleksibel pakkehåndtering, der muliggør praktisk installation af pakker med deres afhængigheder på tværs af alle platforme. Brug 'Anaconda' (conda-pakkehåndtering) til at installere 'conda', som kan bruges til at installere de værktøjer/pakker, der kræves til analysen.
  2. Download 'Anaconda' i henhold til systemkravene fra https://www.anaconda.com/products/individual#Downloads og installer det ved at følge vejledningen i grafisk installationsprogram. Installer alle nødvendige pakker ved hjælp af 'conda' ved at 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 at downloade alle de R-pakker, der bruges i protokollen, skal du skrive følgende kode i R-konsollen (startet på Linux-kommandolinjen ved at 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)
    }
    BEMÆRK: I denne beregningsprotokol vil kommandoer blive givet som enten R Notebook-filer (filer med udvidelse ". Rmd"), R-kodefiler (filer med udvidelse ". R"), eller Linux Bash shell-scripts (filer med udvidelse ".sh"). R Notebook (Rmd) filer skal åbnes i RStudio ved hjælp af File| Åbn Filer..., og individuelle kodestykker (som kan være R-kommandoer eller Bash shell-kommandoer) kør derefter interaktivt ved at klikke på den grønne pil øverst til højre. R-kodefiler kan køres ved at åbne i RStudio eller på Linux-kommandolinjen ved at indlede med "Rscript", f.eks. Rscript example.R. Shell-scripts køres på Linux-kommandolinjen ved at indlede scriptet med kommandoen "sh", f.eks.sh example.sh.

2. Alternativ splejsning (AS) analyse ved hjælp af RNA-seq

  1. Download og forbehandling af data
    BEMÆRK: De kodestykker, der er kommenteret nedenfor, er tilgængelige i den supplerende kodefil "AS_analysis_RNASeq.Rmd", for at følge de enkelte trin interaktivt og leveres også som et bash-script, der skal køres i batch på Linux-kommandolinjen (Sh downloading_data_preprocessing.sh).
    1. Download af rådata.
      1. Download rådata fra Sequence Read Archive (SRA) ved hjælp af kommandoen 'prefetch' fra SRA-værktøjssæt (v2.10.8)17. Giv SRA-tiltrædelses-id'erne i rækkefølge i følgende kommando for at downloade dem parallelt ved hjælp af GNU parallelværktøj18. For at downloade SRA-filer med tiltrædelses-id'er fra SRR10261601 til SRR10261606 parallelt skal du bruge følgende på Linux-kommandolinjen.
        ​seq 10261601 10261606 | parallel prefetch SRR{}
      2. Uddrag fastq-filerne fra arkivet ved hjælp af 'fastq-dump' -funktionen fra SRA-værktøjssæt. Brug GNU parallel og giv navnene på alle SRA-filer sammen.
        ​parallel -j 3 fastq-dump --gzip --skip-technical --read-filter pass --dumpbase --split-e --clip --origfmt {} :::
      3. Download referencegenomet og anmærkningerne til mus (Genome assembly GRCm39) fra www.ensembl.org ved hjælp af 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 kortlægning af aflæsninger til genomsamlingen
      1. Kvalitetskontrol. Vurder kvaliteten af rå læsninger med FASTQC (FASTQ Quality Check v0.11.9)19. Opret en outputmappe, og kør fastqc med parallel på flere fasta-inputfiler. Dette trin genererer en kvalitetsrapport for hver prøve. Undersøg rapporterne for at sikre, at kvaliteten af læsningerne er acceptabel, før du foretager yderligere analyse. (Se brugervejledningen for at forstå rapporterne på https://www.bioinformatics.babraham.ac.uk/projects/fastqc/)
        mkdir fastqc_out
        parallel "fastqc {} -o fastqc_out" ::: $RAW_DATA/*.fastq.gz
        BEMÆRK: Udfør om nødvendigt adaptertrimning med 'cutadapt'20 eller 'trimmomatic'21 for at fjerne sekventering i flankerende adaptere, som varierer afhængigt af RNA-fragmentstørrelse og læselængde. I denne analyse sprang vi dette trin over, da brøkdelen af berørte læsninger var minimal.
      2. Læs justering. Det næste trin i forbehandlingen omfatter kortlægning af læsningerne til referencegenomet. Først skal du opbygge indekset for referencegenomet ved hjælp af funktionen 'genomeGenerate' i STAR22 og derefter justere de rå læsninger til referencen (alternativt er forudbyggede indekser tilgængelige fra STAR-webstedet og kan bruges direkte til justering). Kø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
        BEMÆRK: STAR-justeringsprogrammet genererer og sorterer BAM-filer (Binary Alignment Map) for hver prøve efter læsejustering. Bam-filer skal sorteres, før du fortsætter til yderligere trin.
  2. Forberedelse af Exon-annoteringer.
    1. Kør den supplerende kodefil "prepare_mm_exon_annotation. R" med den downloadede anmærkning i GTF-format (Gene transfer format) for at forberede anmærkningerne. For at køre skal du skrive følgende på Linux-kommandolinjen.
      Rscript prepare_mm_exon_annotation.R annotation.gtf
      BEMÆRK: Den GTF fil indeholder flere exon poster for forskellige isoforms. Denne fil bruges til at "kollapse" de flere transkriptions-id'er for hver exon. Det er et vigtigt skridt at definere exon-tællekasser.
  3. Optælling læser. Det næste trin er at tælle antallet af læsninger, der er kortlagt til forskellige udskrifter / exons. Se supplerende sagsfremstilling: "AS_analysis_RNASeq.Rmd".
    1. Indlæs påkrævede biblioteker:
      packages<- c("Rsubread","tidyverse", "magrittr", "EnhancedVolcano", "edgeR","openxlsx")
      invisible(lapply(packages, library, character.only=TRUE))
    2. Indlæs den behandlede annoteringsfil hentet fra forrige trin (2.2).
      load("mm_exon_anno.RData")
    3. Læs alle bam-filer, der er hentet i trin 2.2.2, som input til 'featureCounts' for at tælle læsninger. Læs mappen, der indeholder bam-filer, ved først at liste hver fil fra den mappe, der slutter med .bam. Brug 'featureCounts' fra Rsubread-pakken, der tager bam-filer og behandlet GTF-annotation (reference) som input til at generere en matrix af tællinger, der er knyttet til hver funktion med rækker, der repræsenterer exons (funktioner) og kolonner, der repræsenterer 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. Udfør derefter ikke-specifik filtrering for at fjerne lavt udtrykte exoner ("ikke-specifik" angiver, at de eksperimentelle tilstandsoplysninger ikke bruges i filtreringen for at undgå selektionsforstyrrelser). Transformér dataene fra rå skala til tællinger pr. million (cpm) ved hjælp af cpm-funktionen fra 'edgeR'-pakke23 , og hold exoner med tællinger, der er større end en indstillelig tærskel (for dette datasæt bruges en cpm) i mindst tre prøver. Fjern også gener med kun en 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, ]
      BEMÆRK: Kontroller de nødvendige parametre for featureCounts, når du bruger forskellige data, for eksempel til single-end-læsninger, indstil 'isPairedEnd = FALSE'. Se RSubread-brugervejledningen for at vælge indstillingerne for dine data, og se afsnittet Diskussion nedenfor.
  4. Analyse af differentiel splejsning og exon-brug. Vi beskriver to alternativer til dette trin: DEXSeq og DiffSplice. Begge kan bruges og give lignende resultater. For konsistens skal du vælge DEXSeq, hvis du foretrækker en DESeq2-pakke til DGE og bruge DiffSplice til en Limma-baseret DGE-analyse. Se supplerende sagsfremstilling: "AS_analysis_RNASeq.Rmd".
    1. Brug af DEXSeq-pakke til differentiel exonanalyse.
      1. Indlæs bibliotek, og opret en eksempeltabel for at definere det eksperimentelle design.
        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")))
        BEMÆRK: Rækkenavnene skal være i overensstemmelse med de bam-filnavne, der bruges af featureCounts til at tælle læsningerne. sampleTable består af detaljer om hver prøve, som omfatter: bibliotekstype og tilstand. Dette er nødvendigt for at definere kontrasterne eller testgruppen til påvisning af differentieret brug.
      2. Forbered exon-informationsfilen. Exon-information i form af GRanges (genomiske intervaller) objekter (https://bioconductor.org/packages/release/bioc/html/GenomicRanges.html) kræves som input for at oprette DEXSeq-objektet i næste trin. Match gen-id'erne med læsetællingerne for at oprette 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. Opret DEXSeq-objekt ved hjælp af DEXSeqDataSet-funktionen. DEXSeq-objektet samler læsetællinger, exon-funktionsoplysninger og eksempeloplysninger. Brug de læsetællinger, der blev genereret i trin 3, og exonoplysningerne fra det forrige trin til at oprette DEXSeq-objektet ud fra tællematrixen. Argumentet sampleData tager et datarammeinput, der definerer prøverne (og deres attributter: bibliotekstype og tilstand), 'design' bruger sampleData til at generere en designmatrix til differentialtesten ved hjælp af modelformelnotation. Bemærk, at et signifikant interaktionsudtryk, condition:exon, indikerer, at brøkdelen af læsninger over et gen, der falder på en bestemt exon, afhænger af den eksperimentelle tilstand, dvs. der er AS. Se DEXSeq-dokumentationen for en fuldstændig beskrivelse af indstilling af modelformlen for mere komplekse eksperimentelle designs. For funktionsinformation kræves exon id'er, tilsvarende gen og transkripter.
        ​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 spredning estimering. Udfør derefter normalisering mellem prøver og estimer variansen af dataene på grund af både Poisson-tællestøj fra den diskrete karakter af RNA-seq og biologisk variabilitet ved hjælp af følgende kommandoer.
        ​dxd %<>% estimateSizeFactors %>% estimateDispersions %T>% plotDispEsts
      5. Test for differentieret brug. Efter estimering af variation testes for differentiel exonbrug for hvert gen og genererer resultaterne.
        dxd%<>%testForDEU%>%estimateExonFoldChanges(fitExpToVar=
        "condition")#Estimate fold changes
        dxr=DEXSeqResults(dxd)
      6. Visualiser splejsningshændelser for udvalgte gener ved hjælp af følgende kommando.
        plotDEXSeq(dxr,"Wnk1", displayTranscripts=TRUE, splicing=TRUE,legend
        =TRUE,cex.axis=1.2,cex=1.3,lwd=2)
        Undersøg R Notebook-filen "AS_analysis_RNASeq.Rmd" for at generere yderligere plots for gener af interesse og for at generere vulkanplotter ved forskellige tærskler.
    2. Brug af diffSplice fra Limma til at identificere differentiel splejsning. Følg filen R Notebook "AS_analysis_RNASeq.Rmd". Sørg for, at trin 2.1-2.3 er blevet fulgt for at forberede inputfiler, før du fortsætter yderligere.
      1. Indlæs biblioteker
        library(limma)
        library(edgeR)
      2. Ikke-specifik filtrering. Matrixen af læsetællinger i punkt 2.3 uddrages. Opret en liste over funktioner ved hjælp af funktionen 'DGEList' fra edgeR-pakken, hvor rækker repræsenterer gener, og kolonner repræsenterer 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, "\\,")
        BEMÆRK: Som et ikke-specifikt filtreringstrin filtreres tællinger efter cpm < 1 in x ud af n prøver, hvor x er det mindste antal replikater under enhver tilstand. n = 6 og x = 3 for dette eksempel data.
      3. Normaliser tællingerne på tværs af prøver med funktionen 'calcNormFactors' fra 'edgeR'-pakken ved hjælp af Trimmet middelværdi af M-værdier (TMM-normaliseringsmetode)24 Den beregner skaleringsfaktorer for at justere biblioteksstørrelser.
        ​dge<-calcNormFactors(dge)
      4. Brug sampleTable som genereret i trin 2.4.1.1, og opret designmatrixen. Designmatrixen karakteriserer designet. Se Limma User Guide (https://www.bioconductor.org/packages/devel/bioc/vignettes/limma/inst/doc/usersguide.pdf) kapitel 8 & 9 for detaljer om designmatricer til mere avancerede eksperimentelle designs.
        Treat<- factor(sampleTable$condition)
        design<- model.matrix(~0+Treat)
        colnames(design) <- levels(Treat)
      5. Tilpas en lineær model pr. exon. Kør 'voom' -funktionen af 'limma' -pakke for at behandle RNA-seq-data for at estimere varians og generere præcisionsvægte for at korrigere for Poisson-tællestøj og omdanne exon-niveautællingerne til log2-tællinger pr. Million (logCPM). Kør derefter lineær modellering ved hjælp af 'lmfit' -funktion for at tilpasse lineære modeller til udtryksdataene for hver exon. Beregn empirisk Bayes-statistik for tilpasset model ved hjælp af 'eBayes' -funktion til at detektere differentiel exon-ekspression. Definer derefter en kontrastmatrix til de eksperimentelle sammenligninger af interesse. Brug 'contrasts.fit' til at opnå koefficienter og standardfejl 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. Differentiel splejsningsanalyse. Kør 'diffSplice' på den tilpassede model for at teste forskellene i exon-brug af gener mellem vildtype og knockout og udforske de bedst rangerede resultater ved hjælp af 'topSplice' -funktion: test = "t" giver en rangordning af AS-exoner, test = "simes" giver en rangordning af 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. Plot resultaterne med funktionen 'plotSplice', hvilket giver gen af interesse for geneidargumentet. Gem de øverste resultater sorteret efter log Fold skift til et objekt og generer et vulkanplot for at udstille exonerne.
        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. Brug af rMATS
      1. Sørg for, at den nyeste version af rMATS v4.1.1 (også kendt som rMATS turbo på grund af den reducerede behandlingstid og mindre hukommelseskrav) er installeret enten ved hjælp af conda eller github (https://github.com/Xinglab/rmats-turbo/releases/download/v4.1.1/rmats_turbo_v4_1_1.tar.gz) i arbejdsmappen. Følg afsnit 4.3 i "AS_analysis_RNASeq.Rmd".
      2. Gå til mappen, der indeholder bam-filer, der er opnået efter kortlægning, og forbered tekstfiler, som krævet af rMATS, for de to betingelser ved at kopiere navnet på bam-filer (sammen med stien) adskilt af ','. Følgende kommandoer skal kø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. Kør rmats.py med de to inputfiler, der blev genereret i det foregående trin, sammen med GTF-filen opnået i 2.1.1.3. Dette genererer en outputmappe 'rmats_out', der indeholder tekstfiler, der beskriver statistik (p-værdier og inklusionsniveauer) for hver splejsningshændelse 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
        BEMÆRK: Referenceanmærkningen i form af en GTF-fil er også påkrævet. Kontroller parametrene, hvis dataene er single-end, og skift -t-indstillingen i overensstemmelse hermed.
      4. Udforskning af rMATS-resultater. Brug biolederpakken 'maser'25 til at udforske rMATS-resultaterne. Indlæs JCEC-tekstfilerne (Junction and exon counts) i 'maser'-objektet, og filtrer resultatet baseret på dækning ved at inkludere mindst fem gennemsnitlige læsninger pr. splejsningshændelse.
        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 af rMATS-resultater. Vælg de signifikante splejsningshændelser ved False Discovery Rate (FDR) 10% og minimum 10% ændring i Procent splejset ind (deltaPSI) ved hjælp af 'topEvents' -funktionen fra 'maser' -pakken. Kontroller derefter genhændelserne for individuelle gener af interesse (prøvegen-Wnk1) og plot PSI-værdier for hver splejsningshændelse af det gen. Generer et vulkanplot ved at angive hændelsestypen.
        #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-plots for det splejsningshændelsesresultat, der opnås med rMATS i form af tekstfiler ved hjælp af pakken 'rmats2shahimiplot'. Kør python-scriptet 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
        BEMÆRK: Denne proces kan være tidskrævende, da den genererer Sashimi-plottet for alle resultaterne i begivenhedsfilen. Vælg de øverste resultater (gennavne og exons) som vist af topEvents-funktionen fra 'maser', og visualiser det tilsvarende Sashimi-plot.

3. Alternativ polyadenylering (APA) analyse ved hjælp af 3 'slutsekventering

  1. Download og forbehandling af data
    BEMÆRK: Se den supplerende R Notebook-fil "APA_analysis_3PSeq_notebook. Rmd" for de komplette kommandoer til dataoverførsel og forbehandlingstrin, eller kør den supplerende bash-fil "APA_data_downloading_preprocessing.sh" på Linux-kommandolinjen.
    1. Download data fra SRA med tiltrædelses-id'erne (1553129 til 1553136).
    2. Trimadaptere og omvendt komplement for at opnå sansestrengsekvensen.
      BEMÆRK: Dette trin er specifikt for det anvendte PolyA-seq-assay.
    3. Kort læser til musegenomsamling ved hjælp af bowtie aligner26.
  2. Forberedelse af annoteringer af pA-websteder.
    BEMÆRK: Behandlingen af pA-webstedsannoteringsfilen udføres først ved hjælp af supplerende R Notebook-fil "APA_analysis_3PSeq_notebook. Rmd" (2.1 - 2.6), og derefter bruge bash-fil "APA_annotation_preparation.sh".
    1. Download pA-webstedsannotering fra PolyASite 2.0-databasen6.
    2. Vælg pA-webstedsannoteringer for at bevare 3'-uoversatte region (UTR) pA-websteder, der er kommenteret som Terminal Exon (TE) eller 1000 nt nedstrøms for en kommenteret terminal exon (DS) til downstream-analyse.
    3. Få pA-site toppe. Anker ved hvert pA-spaltningssted, og visualiser den gennemsnitlige læsedækning ved hjælp af sengeredskaber og deeptools27,28. Resultaterne viste, at toppene af de kortlagte læsninger hovedsageligt var spredt inden for ~ 60 bp opstrøms for spaltningsstederne (figur 5 og supplerende figur 5). Derfor blev koordinaterne for pA-websteder udvidet fra annotationsfilen til 60 bp opstrøms for deres spaltningssteder. Afhængigt af den specifikke 3'-slutsekventeringsprotokol, der anvendes, skal dette trin optimeres til andre assays end PolyA-seq.
  3. Optælling af læsninger
    1. Forbered annoteringsfilen for pA-websteder.
      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. Anvend 'featureCounts' for at erhverve rå tællinger. Gem tælletabellen som RData-filen "APA_countData.Rdata" til APA-analyse ved hjælp af forskellige værktøjer.
      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")
      BEMÆRK: Vær opmærksom på at ændre nogen af de parametre, der er angivet i funktionen 'featureCounts'. Rediger parameteren 'strengspecifik' for at sikre, at den er i overensstemmelse med sekventeringsretningen for det anvendte 3'-endesekventeringsassay (empirisk vil visualisering af dataene i en genombrowser over gener på plus- og minusstrenge afklare dette).
    3. Anvend ikke-specifik filtrering af countData. Filtrering kan forbedre den statistiske robusthed betydeligt i differentierede test af brugen af pA-websteder. For det første fjernede vi de gener med kun et pA-sted, hvor differentieret pA-webstedsbrug ikke kan defineres. For det andet anvender vi ikke-specifik filtrering baseret på dækning: tællinger filtreres med cpm mindre end 1 in x ud af n prøver, hvor x er det mindste antal replikater under enhver tilstand. N = 8 og x = 2 for dette eksempel data.
      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. Analyse af brugen af differentieret polyadenyleringssted ved hjælp af DEXSeq- og diffSplice-rørledninger.
    1. Brug af DEXSeq-pakken
      BEMÆRK: Da en kontrastmatrix ikke kan defineres for DEXSeq-rørledningen, skal den differentielle APA-analyse af hver to eksperimentelle betingelser udføres separat. Den differentielle APA-analyse af tilstanden WT og tilstanden DKO udføres som et eksempel for at forklare proceduren. Der henvises til tillægssagen "APA_analysis_3PSeq_notebook. Rmd" til den trinvise arbejdsgang i dette afsnit og den differentielle APA-analyse af andre kontraster.
      1. Indlæs bibliotek, og opret en eksempeltabel for at definere det eksperimentelle design.
        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. Forbered pA-webstedets informationsfil ved hjælp af 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. Brug de læsetællinger, der blev genereret i trin 3.3, og pA-webstedsoplysningerne fra det forrige trin til at oprette 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 kontrastparret ved at definere niveauerne af betingelser i DEXSeq-objektet.
        dxd1$condition<- factor(dxd1$condition, levels= c("WT", "DKO"))
        # The contrast pair will be "DKO - WT"
      5. Normalisering og spredning estimering. I lighed med RNA-seq-data udfører 3'-slutsekventeringsdata normalisering mellem prøver (kolonnevis median af forhold for hver prøve) ved hjælp af funktionen 'estimateSizeFactors' og estimerer variationen af dataene ved hjælp af funktionen 'estimateDispersions' og visualiser derefter dispersionsestimeringsresultatet ved hjælp af funktionen 'plotDispEsts'.
        ​dxd1 %<>% estimateSizeFactors %>% estimateDispersions %T>% plotDispEsts
      6. Test af differentiel pA-stedbrug for hvert gen ved hjælp af funktionen 'testForDEU', og estimer derefter foldændringen af pA-stedbrug ved hjælp af funktionen 'estimateExonFoldChanges'. Kontroller resultaterne ved hjælp af funktionen 'DEXSeqResults' og indstil 'FDR < 10%' som kriterium for signifikant differentierede pA-steder.
        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 af de differentielle pA-stedbrugsresultater ved hjælp af differentielle APA-plots genereret af funktionen 'plotDEXSeq' og vulkanplot af funktionen '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. Brug af diffSplice-pakke. Se den supplerende R Notebook-fil "APA_analysis_3PSeq_notebook. Rmd" til den trinvise arbejdsgang i dette afsnit.
      1. Definer kontraster af interesse for differentiel pA-brugsanalyse.
        BEMÆRK: Dette trin skal udføres efter konstruktion og behandling af DGEList-objektet, som er inkluderet 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 af kontraster af interesse (her "DKO - WT") ved hjælp af differentielle APA-plots efter funktionen 'plotSplice' og vulkanplots med funktionen 'EnhancedVolcano'. Se filen "APA_analysis_3PSeq_notebook. Rmd" 4.2.7 - 4.2.9 til visualisering af 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")

Subscription Required. Please recommend JoVE to your librarian.

Representative Results

Efter at have kørt ovenstående trinvise arbejdsgang er AS- og APA-analyseoutput og repræsentative resultater i form af tabeller og dataplot, der genereres som følger.

SOM:
Det vigtigste output af AS-analysen (supplerende tabel 1 for diffSplice; Tabel 2 for DEXSeq) er en liste over exoner, der viser differentieret brug på tværs af betingelser, og en liste over gener, der viser signifikant samlet splejsningsaktivitet af en eller flere af dens bestanddele, rangeret efter statistisk signifikans. Supplerende tabel 1, fane 2 viser signifikante exoner med kolonner, der viser differentieret FC af exon versus hvile, pr. exon ujusteret p-værdi og justeret p-værdi (Benjamini-Hockberg-korrektion). Tærskelværdi på de justerede p-værdier vil give et sæt exoner med defineret FDR. Supplerende tabel 1, fane 3 viser en rangordnet liste over gener, der viser betydningen af den samlede splejsningsaktivitet, med en kolonne, der viser justeret p-værdi på genniveau beregnet ved hjælp af Simes-metoden. Lignende data er vist i tabel 2 for DEXSeq. Supplerende figur 1 og supplerende figur 2 viser differentielt splejsningsmønster i Mbnl1-, Tcf7- og Lef1-gener, som er blevet eksperimentelt valideret i den offentliggjorte artikel, der præsenteres med dataene15. Forfatterne har vist eksperimentel validering af fem gener - Mbnl1, Mbnl2, Lef1, Tcf7 og Ncor2. Vores tilgang opdagede differentielle splejsningspattens i alle disse gener. Her præsenterer vi FDR-niveauerne for hvert gen ved hjælp af henholdsvis DEXSeq, diffSplice og rMATS som opnået i supplerende tabel 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 mellem output opnået fra tre forskellige værktøjer og illustrerer alternative splejsningsmønstre i Wnk1-genet . Vulkanplottene er vist i figur 2A (diffSplice) og figur 2B (DEXSeq). Yderligere tre højt rangerede gener er vist i supplerende figur 1 (diffSplice) og supplerende figur 2 (DEXSeq). Y-aksen viser statistisk signifikans (-log10 P-værdier), og x-aksen viser effektstørrelse (foldændring). Gener placeret øverst til venstre eller højre kvadrant indikerer betydelige FC og stærke statistiske beviser for ægte forskelle.

Figure 2
Figur 2. Sammenligning af alternative splejsningsresultater opnået fra diffSplice, DEXSeq og rMATS. |
(A) Vulkanplot (til venstre) af RNA-Seq fra Limma diffSplice analyse: X-aksen viser log exon fold ændring; Y-aksen viser -log10 p-værdi. Hvert punkt svarer til en exon. Vandret stiplet linje ved p-værdi = 1E-5; lodrette stiplede linjer ved to gange ændring (FC). Røde exons viser betydelig FC og statistisk signifikans. Differentiel splejsningsplot (højre): Splejsningsmønstre udstilles for et eksempel på gen Wnk1, hvor x-aksen viser exon id pr. transkription; y-aksen viser exon relativ log fold ændring (forskellen mellem exon's logFC og den samlede logFC for alle andre exons). Exons fremhævet med rødt viser statistisk signifikant differentiel ekspression (FDR < 0,1).
(B) Vulkanplot (venstre) og differentiel exonbrug (højre) af RNA-Seq opnået ved DEXSeq-analyse. Wnk1-genet viser signifikant differentiel brug af exoner mellem WT og Mbnl1 knock-out fremhævet med pink, hvilket svarer til de samme differentielle exoner i (A).
(C) Vulkanplot (venstre) og Sashimi-plot (højre) for Wnk1 opnået ved rMATS-analyse. Vulkanplot, der viser signifikant sprunget (kassette) exon (SE) begivenhed i vildtype sammenlignet med knockout ved 10% FDR med ændring i procent splejset i (PSI eller ΔΨ) værdier > 0,1. X-aksen viser ændring i PSI-værdier på tværs af betingelser, og y-aksen viser log P-værdi. Sashimi-plottet viser en sprunget exon-begivenhed i Wnk1-genet, svarende til en signifikant differentiel exon i (A) og (B). Hver række repræsenterer en RNA-Seq-prøve: tre replikater af vildtype og Mbnl1 knock-out. Højden viser læsedækning i RPKM, og forbindelsesbuerne viser krydslæsninger på tværs af exons. Annoteret genmodel alternative isoformer er vist i bunden af plottet. Det nederste panel af C illustrerer krydslæsninger, der bruges til at beregne PSI-statistikken.
(D) Venn-diagram, der sammenligner antallet af signifikante differentielle exoner opnået ved de forskellige metoder. Klik her for at se en større version af denne figur.

Figur 2 Et (højre panel) viser en diagrammatisk visning af exonforskelle i et af de højest rangerede gener, der viser logFC på y-aksen og exon-tallet på x-aksen. Dette eksempel viser en kassette exon varierende mellem betingelser for gen Wnk1. Det differentielle exonbrugsplot fra DEXSeq viser tegn på differentiel splejsning på fem exon-steder nær Wnk1.6.45. De fremhævede exoner i pink vil sandsynligvis blive splejset ud i Mbnl1 KO-prøver sammenlignet med WT. Disse exoner supplerer de resultater, der opnås ved diffSplice, som viser et lignende mønster på den specifikke genomiske position. Flere eksempler er vist i supplerende figur 1 og supplerende figur 2. En mere detaljeret visning for at bekræfte interessante resultater kan gives ved at sammenligne dækningsspor (wiggle) i RPM (Læser pr. Million) enheder i UCSC (University of Santa Cruz) eller IGV (Integrative Genomics Viewer) genombrowsere (ikke vist) sammen med visuel korrelation med andre spor af interesse, såsom kendte genmodeller, bevarelse og andre genom-dækkende assays.

rMATS-outputtabel viser vigtige alternative splejsningshændelser kategoriseret efter type (supplerende tabel 3). Figur 2C viser et vulkanplot af gener, der er alternative splejset, med effektstørrelsen målt ved differential "procent splejset i" (PSI eller ΔΨ) statistik på11.

PSI henviser til procentdelen af læsninger, der er i overensstemmelse med inkluderingen af en kassette exon (dvs. læser kortlægning til kassetten exon selv eller kryds læser overlappende exon) sammenlignet med læsninger i overensstemmelse med exon udelukkelse, dvs. kryds læser på tværs af tilstødende opstrøms og nedstrøms exons (Det nederste panel i figur 2C). Det højre panel i figur 2C viser sashimi-plot af Wnk1-genet med differentiel splejsningshændelse overlejret på dækningsspor for genet med en sprunget exon i Mbnl1 KO. Buer, der forbinder exoner, viser antallet af krydslæsninger (læser krydser en splejset intron). Forskellige faner i supplerende tabel 3 viser signifikante læsninger af hver type hændelse, der spænder over kryds med exongrænser (krydstællinger og exontællinger (JCEC)). Figur 2D sammenligner de signifikante differentielt splejsede exoner, der er registreret af de tre værktøjer.

Figure 3
Figur 3. Alternative splejsningshændelser erhvervet ved rMATS-analyse. A) Typer af AS-begivenheder. Dette tal er tilpasset fra rMATS dokumentation11, der forklarer splejsningsbegivenheden med konstitutive og alternativt splejsede exoner. Mærket med antallet af hver begivenhed på FDR 10%. B) Sashimi-plot af Add3-genet, der udviser sprunget exon (SE). C) Sashimi-plot af Baz2b-genet, der udviser alternativt 5'-splejsningssted (A5SS). D) Sashimi-plot af Lsm14b-genet, der udviser alternativt 3'-splejsningssted (A3SS). E) Sashimi-plot af Mta1-genet, der udviser gensidigt udelukkende exoner (MXE). F) Sashimi-plot af Arpp21-genet, der udviser bevaret intron (RI). Røde rækker repræsenterer tre replikater af vildtype, og orange rækker repræsenterer Mbnl1 knock-out replikater. X-aksen svarer til genomiske koordinater og strenginformation, y-aksen viser dækning i RPKM. Klik her for at se en større version af denne figur.

Figur 3 illustrerer typer af splejsningshændelser SE, A5SS, A3SS, MXE og RI ved hjælp af Sashimi-plots af de vigtigste signifikante gener af disse begivenheder. Ved sammenligning af de tre replikater af både WT og Mbnl1_KO blev der påvist i alt 1272 SE-hændelser, 130 A5SS, 116 A3SS, 215 MXE og 313 RI-hændelser ved FDR 10%. Sashimi plot illustrerer typen af begivenhed ved hjælp af top gener som et eksempel.

APA:
Outputtet fra APA-analysen svarer til AS-analysen på exon-niveau. Der findes en tabel over topgener rangeret efter differentiel APA-aktivitet i 3'UTR (supplerende tabel 4 og supplerende tabel 5). Figur 4A viser vulkanplotterne af gener ved differentiel APA-aktivitet i 3'UTR'er genereret ved hjælp af diffSplice og DEXSeq separat. Figur 4B viser Venn-plottet, der sammenligner de signifikant differentierede pA-brugsresultater erhvervet fra forskellige rørledninger. Figur 4C og 4D viser den diagrammatiske repræsentation af differentiel pA-stedbrug i 3'UTR af genet Fosl2 (figur 4C) og Papola (figur 4D) genereret ved hjælp af både diffSplice og DEXSeq, som er eksperimentelt valideret til at vise signifikant distal til proksimalt skift (Fosl2) og signifikant proksimalt til distalt skift (Papola) af pA-webstedsbrug i DKO, henholdsvis. Flere eksempler er vist i supplerende figur 3 og supplerende figur 4.

Figure 4
Figur 4. Alternative polyadenyleringsdiagrammer ved diffSplice og DEXSeq. A) Vulkandiagrammer af PolyA-seq-data genereret ved hjælp af diffSplice og DEXSeq. X-akse viser log pA site fold ændring; y-aksen viser -log10 p-værdi. Hvert punkt svarer til et pA-sted. Vandret stiplet linje ved p-værdi = 1E-5; lodrette stiplede linjer ved 2-fold FC. Røde exons viser betydelig FC og statistisk signifikans. B) Venn-plot , der sammenligner de betydelige differentierede pA-webstedsbrugsresultater, der er erhvervet fra forskellige rørledninger. C-D) Differentielle APA-plots genereret ved hjælp af diffSplice og DEXSeq, der viser de proksimale, interne og distale pA-steder for Fosl2 - og Papola-genet . Tallene genereres ved samme funktion som figur 2 (B), men med pA-steder, der erstatter exons. Klik her for at se en større version af denne figur.

Figur 5 er et diagnostisk plot for at bekræfte den forventede læsefordeling omkring kommenterede pA-spaltningssteder for det anvendte PolyA-seq-assay. Det viser den gennemsnitlige dækning i flankerende regioner forankret på kendte pA-spaltningssteder på genomniveau. I dette tilfælde visualiseres den forventede bunke af læsninger opstrøms for webstederne. Læsefordelingerne forankret på pA-steder for alle PolyA-seq-prøver er vist i supplerende figur 5.

Figure 5
Figur 5. Gennemsnitlig dækningsplot omkring pA-spaltningssteder. Spaltningsstedet for PolyA-seq-data vises med lodret stiplet linje. X-akse viser basisposition i forhold til pA-spaltningssteder, op til 100 nukleotider opstrøms og nedstrøms; y-aksen viser den gennemsnitlige læsedækning over alle pA-spaltningssteder, normaliseret efter biblioteksstørrelse i CPM. Klik her for at se en større version af denne figur.

De differentielle APA-resultater af forskellige kontraster genereret af den samme rørledning kan sammenlignes og verificeres ved at visualisere læsedækningen af repræsentative signifikant differentierede pA-steder i genombrowseren. Figur 6A er Venn-plottet, der sammenligner den signifikant differentierede pA-webstedsbrug af forskellige kontraster erhvervet fra diffSplice. Figur 6B-D er IGV-snapshots af læsedækningen på pA-steder for forskellige gener, som viser mønstrene i overensstemmelse med dem, der blev opdaget i APA-analysen ved hjælp af diffSplice. Figur 6B validerer det signifikante proksimale til distale skift af pA-stedbrug for gen Paip2, som entydigt detekteres i kontrasten DKO vs WT, men ikke i andre to kontraster KD vs WT og Ctr vs WT. Figur 6C validerer det signifikante distale til proksimale skift af pA-stedbrug for gen Ccl2 detekteret entydigt i kontrasten KD vs WT, mens figur 6D validerer den signifikante differentielle pA-brug af alle kontraster for genet Cacna2d1.

Figure 6
Figur 6. Kontrast sammenligning og verifikation af diffSplice resultater. A) Venn-diagram, der sammenligner signifikante differentielle pA-webstedsbrugsresultater af forskellige kontraster erhvervet fra diffSplice. B-D) IGV snapshot visualisering pA toppe dækning af gener Paip2, Ccl2 og Cacna2d1 på tværs af betingelser. Hvert spor repræsenterer læsedækningen i en bestemt tilstand. Klik her for at se en større version af denne figur.

Supplerende figur 1. RNA-Seq analyse af differentiel splejsning med Limma diffSplice. (A) Vulkanplot af RNA-Seq fra Limma diffSplice analyse: X-aksen viser log exon fold ændring; Y-aksen viser -log10 p-værdi. Hvert punkt svarer til en exon. Vandret stiplet linje ved p-værdi = 1E-5; lodrette stiplede linjer ved dobbelt ændring (FC). Røde exons viser betydelig FC og statistisk signifikans. (B-D) Differentierede splejsningsområder: Splejsningsmønstre udstilles for eksempel generne Mbnl1, Tcf7 og Lef1, hvor x-aksen viser exon id pr. transkription; y-aksen viser exon relativ log fold ændring (forskellen mellem exon's logFC og den samlede logFC for alle andre exons). Exons fremhævet med rødt viser statistisk signifikant differentiel ekspression (FDR < 0,1). Klik her for at downloade denne fil.

Supplerende figur 2. RNA-Seq analyse af differentiel exon brug med DEXSeq. (A) Vulkan plot. (B-D) Differentiel exon-brug af RNA-Seq opnået ved DEXSeq-analyse. Generne Mbnl1, Tcf7 og Lef1 viser henholdsvis signifikant differentieret brug af exoner mellem WT og Mbnl1 knock-out fremhævet med pink, hvilket svarer til de samme differentielle exoner i supplerende figur 1. Klik her for at downloade denne fil.  

Supplerende figur 3. Alternative polyadenyleringsplots ved diffSplice. A) Vulkanplots af PolyA-seq-data genereret ved hjælp af diffSplice i tre kontrastpar opnået fra musen PolyA-seq-data, herunder dobbelt knockout (DKO) vs vildtype (WT), knock-down (KD) vs WT og kontrol (Ctrl) vs WT. X-akse viser log pA site fold ændring; y-aksen viser -log10 p-værdi. Hvert punkt svarer til et pA-sted. Vandret stiplet linje ved p-værdi = 1E-5; lodrette stiplede linjer ved 2-fold FC. Røde pA-websteder viser betydelig FC og statistisk signifikans. B) Differentielle APA-plots genereret ved hjælp af diffSplice, der viser de proksimale, interne og distale pA-steder for de højt rangerede gener S100a7a, Tpm1 og Smc6. Klik her for at downloade denne fil.  

Supplerende figur 4. Differentiel pA-brugsanalyse af DEXSeq-pipeline. A) Vulkanplots af PolyA-seq-data genereret ved hjælp af DEXSeq i tre par opnået fra musen PolyA-seq-data, herunder dobbelt knockout (DKO) vs vildtype (WT), knock-down (KD) vs WT og kontrol (Ctrl) vs WT. X-akse viser log pA site fold ændring; y-aksen viser -log10 p-værdi. Hvert punkt svarer til et pA-sted. Vandret stiplet linje ved p-værdi = 1E-5; lodrette stiplede linjer ved 2-fold FC. Røde pA-websteder viser betydelig FC og statistisk signifikans. B) Differentielle APA-plots genereret ved hjælp af DEXSeq, der viser de proksimale, interne og distale pA-steder for de højt rangerede gener S100a7a, Tpm1 og Smc6. Klik her for at downloade denne fil.  

Supplerende figur 5. Gennemsnitlig dækningsplot og varmekort omkring pA-spaltningssteder.  Dækning vist for fire tilstande, med gener på fremadgående og omvendte tråde vist separat. X-akse viser basisposition i forhold til pA-spaltningssteder, op til 100 nukleotider opstrøms og nedstrøms; y-akse henviser til den gennemsnitlige dækning ved tilsvarende relative basispositioner på tværs af alle pA-spaltningssteder. Heatmaps giver en alternativ visning, hvor hvert pA-spaltningssted vises som en separat række på x-aksen, sorteret efter dækning. Intensitet viser læsedækning (se forklaring). Klik her for at downloade denne fil.

Supplerende tabel 1. diffSplice output af AS-analysen. Første fane definerer kolonnenavnene for de oprindelige output, der præsenteres i anden fane (exon-niveau) og tredje (genniveau). Klik her for at downloade denne tabel.

Supplerende tabel 2. DEXSeq-output af AS-analysen. Første fane definerer kolonnenavnene for det oprindelige output, der præsenteres i anden fane (exon-niveau) og tredje (genniveau). Klik her for at downloade denne tabel.

Supplerende tabel 3. rMATS-output af AS-analysen. Første fane definerer kolonnenavnene for oversigtsfilen (fane 2) og JCEC-filerne for hver hændelse (fane 3-7). Klik her for at downloade denne tabel.

Supplerende tabel 4. diffSplice output af APA-analysen. Første fane definerer kolonnenavnene for det oprindelige output, der præsenteres i anden (pa-webstedsniveau) og tredje (genniveau) faner. Klik her for at downloade denne tabel.

Supplerende tabel 5. DEXSeq-output af APA-analysen. Første fane definerer kolonnenavnene for det oprindelige output, der præsenteres i anden (pa-webstedsniveau) og tredje (genniveau) faner. Klik her for at downloade denne tabel.

Supplerende tabel 6. En oversigt over antallet af væsentligt ændrede exoner for AS- og pA-websteder for APA. Klik her for at downloade denne tabel.

Supplerende tabel 7. En oversigt over værktøjer og pakker, der blev brugt i AS/APA-analysen. Klik her for at downloade denne tabel.

Subscription Required. Please recommend JoVE to your librarian.

Discussion

I denne undersøgelse evaluerede vi exon-baserede og begivenhedsbaserede tilgange til at detektere AS og APA i bulk RNA-Seq og 3 'slutsekventeringsdata. De exonbaserede AS-tilgange producerer både en liste over differentielt udtrykte exoner og en rangordning på genniveau ordnet efter den statistiske signifikans af den samlede differentielle splejsningsaktivitet på genniveau (tabel 1-2, 4-5). For diffSplice-pakken bestemmes differentiel brug ved at tilpasse vægtede lineære modeller på exon-niveau for at estimere differential log fold-change af en exon mod den gennemsnitlige log fold-ændring af de andre exoner inden for samme gen (kaldet per exon FC). Den statistiske signifikans på genniveau beregnes ved at kombinere individuelle signifikanstest på exon-niveau til en gen-klog test ved hjælp af Simes-metoden10. Denne rangordning efter differentiel splejsningsaktivitet på genniveau kan efterfølgende bruges til at udføre en gensætberigelsesanalyse (GSEA) af de involverede nøgleveje10. DEXSeq bruger en lignende strategi ved at tilpasse en generaliseret lineær model til måling af differentiel exon-brug, selvom den adskiller sig i visse trin såsom filtrering, normalisering og dispersionsestimering. Ved sammenligning af de top 500 rangerede exoner, der viser AS-aktivitet og APA ved hjælp af DEXSeq og DiffSplice, fandt vi en overlapning på henholdsvis 310 exoner og 300 pA-steder, hvilket demonstrerede overensstemmelsen mellem de to exon-baserede tilgange, hvilket også blev demonstreret i en tidligere undersøgelse 29. Det anbefales at anvende en kombination af både en exonbaseret (enten DEXSeq eller diffSplice) og en hændelsesbaseret tilgang til omfattende detektion og klassificering af AS. For APA kan brugerne vælge enten DEXSeq eller diffSplice: begge metoder har vist sig at fungere godt på tværs af en lang række transkriptomiske eksperimenter29.

Ved forberedelsen af RNA-seq-biblioteket til en AS-analyse er det vigtigt at bruge en strengspecifik bulk-RNA-seq-protokol8, da en stor del af generne i hvirveldyrgenomer overlapper hinanden på forskellige strenge, og en ikke-strengspecifik protokol ikke er i stand til at skelne mellem disse overlappende regioner, hvilket forvirrer den endelige exondetektion. En anden overvejelse er læsedybde med splejsningsanalyser, der kræver dybere sekventering end DGE, f.eks. 30-60 millioner læsninger pr. prøve, mod 5-25 millioner læsninger pr. prøve for DGE (https://sapac.support.illumina.com/bulletins/2017/04/considerations-for-rna-seq-read-length-and-coverage-.html). Alle de værktøjer, der demonstreres i protokollen, understøtter både single-end og paired-end sekventeringsdata. Hvis kun kendte genannoteringer bruges til at detektere krydslæsninger, kan enkeltendede kortere læsninger (≥ 50 bp) bruges, selvom de novo-detektion af nye splejsningskryds drager fordel af parret ende og længere (≥ 100bp) læser30,31. Valget af RNA-ekstraktionsprotokol - enten polyA-selektion eller rRNA-udtømning - afhænger af kvaliteten af RNA og det eksperimentelle spørgsmål - se31 for en diskussion. Afhængigt af detaljerne i bibliotekskonstruktionen kræves der ændringer af de eksempelscripts, der er angivet her for parametrene for læsejustering, funktionstælling og rMATS. Ved beregning af de indledende læsetællinger på exonniveau ved hjælp af featureCounts eller lignende metoder skal man sørge for at konfigurere funktionsindstillingerne korrekt for tællinger og stranding: i featureCounts indstiller vi argumentet "strengspecifik" korrekt for den anvendte strengspecifikke RNA-seq-protokol; og for kvantificering på exon-niveau forventes det, at en læsning vil kortlægge over tilstødende exons, og derfor indstiller vi parameteren allowMultiOverlap til TRUE. For APA er der forskellige 3'-endesekventeringsprotokoller6, som varierer i den nøjagtige placering af toppe i forhold til pA-stedet. For vores eksempeldata bestemmer vi, at toppen er 60 bp opstrøms for pA-stedet som vist i figur 5, og denne analyse skal tilpasses til andre 3'-slutsekventeringsprotokoller.

I denne protokol begrænser vi omfanget til diskussion af differentielle analyser på niveau med individuelle exoner og splejsningsbegivenheder bestående af tilstødende exon-intron-kombinationer. Vi diskuterer ikke klassen af analyser baseret på isoform de novo rekonstruktion såsom Cufflinks, Cuffdiff32, RSEM 33, Kallisto34, der sigter mod at detektere og kvantificere den absolutte og relative ekspression af hele alternative isoformer. De exon- og hændelsesbaserede metoder er mere følsomme til påvisning af individuelle splejsningshændelser30 og giver i mange tilfælde alle de oplysninger, der er nødvendige for yderligere analyse, uden behov for kvantificering på isoformniveau.

Den seneste version af kildefilerne i denne protokol er tilgængelige på https://github.com/jiayuwen/AS_APA_JoVE

Subscription Required. Please recommend JoVE to your librarian.

Disclosures

Forfatterne har intet at afsløre.

Acknowledgments

Denne undersøgelse blev støttet af et australsk forskningsråd (ARC) Future Fellowship (FT16010043) og ANU Futures Scheme.

Materials

Name Company Catalog Number Comments
Not relevent for computational study

DOWNLOAD MATERIALS LIST

References

  1. Katz, Y., Wang, E. T., Airoldi, E. M., Burge, C. B. Analysis and design of RNA sequencing experiments for identifying isoform regulation. Nature Methods. 7 (12), 1009-1015 (2010).
  2. Wang, Y., et al. Mechanism of alternative splicing and its regulation. Biomedical Reports. 3 (2), 152-158 (2015).
  3. Mehmood, A., et al. Systematic evaluation of differential splicing tools for RNA-seq studies. Briefings in Bioinformatics. 21 (6), 2052-2065 (2020).
  4. Movassat, M., et al. Coupling between alternative polyadenylation and alternative splicing is limited to terminal introns. RNA Biology. 13 (7), 646-655 (2016).
  5. Tian, B., Manley, J. L. Alternative polyadenylation of mRNA precursors. Nature Reviews Molecular Cell Biology. 18 (1), 18-30 (2017).
  6. Herrmann, C. J., et al. PolyASite 2.0: a consolidated atlas of polyadenylation sites from 3' end sequencing. Nucleic Acids Research. 48 (1), 174-179 (2020).
  7. Liu, R., Loraine, A. E., Dickerson, J. A. Comparisons of computational methods for differential alternative splicing detection using RNA-seq in plant systems. BMC Bioinformatics. 15 (1), 364 (2014).
  8. Conesa, A., et al. A survey of best practices for RNA-seq data analysis. Genome Biology. 17 (1), 13 (2016).
  9. Anders, S., Reyes, A., Huber, W. Detecting differential usage of exons from RNA-seq data. Genome Research. 22 (10), 2008-2017 (2012).
  10. Ritchie, M. E., et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Research. 43 (7), 47 (2014).
  11. Shen, S., et al. rMATS: Robust and flexible detection of differential alternative splicing from replicate RNA-Seq data. Proceedings of the National Academy of Sciences. 111 (51), 5593-5601 (2014).
  12. Mehmood, A., et al. Systematic evaluation of differential splicing tools for RNA-seq studies. Briefings in bioinformatics. 21 (6), 2052-2065 (2020).
  13. Kanitz, A., et al. Comparative assessment of methods for the computational inference of transcript isoform abundance from RNA-seq data. Genome biology. 16 (1), 1-26 (2015).
  14. Love, M. I., Huber, W., Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology. 15 (12), 550 (2014).
  15. Sznajder, L. J., et al. Loss of MBNL1 induces RNA misprocessing in the thymus and peripheral blood. Nature Communications. 11, 1-11 (2020).
  16. Batra, R., et al. Loss of MBNL leads to disruption of developmentally regulated alternative polyadenylation in RNA-mediated disease. Molecular Cell. 56 (2), 311-322 (2014).
  17. Leinonen, R., Sugawara, H., Shumway, M., et al. The sequence read archive. Nucleic acids research. 39, suppl_1 19-21 (2010).
  18. Tange, O. GNU parallel-the command-line power tool. 36, 42-47 (2011).
  19. Andrews, S. FastQC: a quality control tool for high throughput sequence data. Bioinformatics. , Available from: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/2010 (2011).
  20. Martin, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet journal. 17 (1), 10-12 (2011).
  21. Bolger, A. M., Lohse, M., Usadel, B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 30 (15), 2114-2120 (2014).
  22. Dobin, A., et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 29 (1), Oxford, England. 15-21 (2013).
  23. Robinson, M. D., McCarthy, D. J., Smyth, G. K. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 26 (1), 139-140 (2010).
  24. Robinson, M. D., Oshlack, A. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biology. 11 (3), 25 (2010).
  25. Veiga, D. F. T. maser: Mapping Alternative Splicing Events to pRoteins. R package version 1.4.0. , (2019).
  26. Langmead, B., Trapnell, C., Pop, M., Salzberg, S. L. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biology. 10 (13), 25 (2009).
  27. Quinlan, A. R., Hall, I. M. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 26 (6), 841-842 (2010).
  28. Ramírez, F., Dündar, F., Diehl, S., Grüning, B. A., Manke, T. deepTools: a flexible platform for exploring deep-sequencing data. Nucleic acids research. 42 (1), 187-191 (2014).
  29. Merino, G. A., Conesa, A., Fernández, E. A. A benchmarking of workflows for detecting differential splicing and differential expression at isoform level in human RNA-seq studies. Briefings in bioinformatics. 20 (2), 471-481 (2019).
  30. Chhangawala, S., Rudy, G., Mason, C. E., Rosenfeld, J. A. The impact of read length on quantification of differentially expressed genes and splice junction detection. Genome biology. 16 (1), 1-10 (2015).
  31. Conesa, A., et al. A survey of best practices for RNA-seq data analysis. Genome Biol. 17, 13 (2016).
  32. Trapnell, C., et al. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc. 7 (3), 562-578 (2012).
  33. Li, B., Dewey, C. N. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 12, 323 (2011).
  34. Bray, N. L., Pimentel, H., Melsted, P., Pachter, L. Near-optimal probabilistic RNA-seq quantification. Nat Biotechnol. 34 (5), 525-527 (2016).

Tags

Biologi udgave 172
Identifikation af alternativ splejsning og polyadenylering i RNA-seq-data
Play Video
PDF DOI DOWNLOAD MATERIALS LIST

Cite this Article

Dixit, G., Zheng, Y., Parker, B.,… More

Dixit, G., Zheng, Y., Parker, B., Wen, J. Identification of Alternative Splicing and Polyadenylation in RNA-seq Data. J. Vis. Exp. (172), e62636, doi:10.3791/62636 (2021).

Less
Copy Citation Download Citation Reprints and Permissions
View Video

Get cutting-edge science videos from JoVE sent straight to your inbox every month.

Waiting X
Simple Hit Counter