Waiting
Login processing...

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

Biology

Identifiering av alternativ skarvning och 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 skarvning (AS) och alternativ polyadenylation (APA) utökar mångfalden av transkriptisoformer och deras produkter. Här beskriver vi bioinformatiska protokoll för att analysera bulk RNA-seq och 3 'end sequencing assays för att detektera och visualisera AS och APA som varierar över experimentella förhållanden.

Abstract

Förutom den typiska analysen av RNA-Seq för att mäta differentiellt genuttryck (DGE) över experimentella / biologiska förhållanden, kan RNA-seq-data också användas för att utforska andra komplexa regleringsmekanismer på exonnivå. Alternativ skarvning och polyadenylation spelar en avgörande roll för en gens funktionella mångfald genom att generera olika isoformer för att reglera genuttryck på post-transkriptionsnivå, och att begränsa analyser till hela gennivån kan missa detta viktiga regleringsskikt. Här demonstrerar vi detaljerade stegvisa analyser för identifiering och visualisering av differentiell användning av exon- och polyadenyleringsplatser över förhållanden, med hjälp av Bioconductor och andra paket och funktioner, inklusive DEXSeq, diffSplice från Limma-paketet och rMATS.

Introduction

RNA-seq har använts i stor utsträckning genom åren vanligtvis för att uppskatta differentiellt genuttryck och genupptäckt1. Dessutom kan den också användas för att uppskatta varierande användning på exonnivå på grund av genuttryck för olika isoformer, vilket bidrar till en bättre förståelse av genreglering på post-transkriptionsnivå. Majoriteten av eukaryota gener genererar olika isoformer genom alternativ skarvning (AS) för att öka mångfalden av mRNA-uttryck. AS-händelser kan delas in i olika mönster: hoppning av fullständiga exoner (SE) där en ("kassett") exon helt tas bort ur transkriptet tillsammans med dess flankerande introner; alternativt (givare) 5' val av skarvplats (A5SS) och alternativt 3' (acceptor) skarvplatsval (A3SS) när två eller flera skarvställen finns i vardera änden av en exon; retention av introner (RI) när en intron behålls inom det mogna mRNA-transkriptet och ömsesidig uteslutning av exonanvändning (MXE) där endast en av de två tillgängliga exonerna kan behållas vid en tidpunkt 2,3. Alternativ polyadenylation (APA) spelar också en viktig roll för att reglera genuttryck med hjälp av alternativa poly (A) -ställen för att generera flera mRNA -isoformer från ett enda transkript4. De flesta polyadenyleringsställen (pAs) är belägna i 3'-oöversatta regionen (3'UTR), vilket genererar mRNA-isoformer med olika 3'UTR-längder. Eftersom 3'UTR är det centrala navet för att känna igen regulatoriska element, kan olika 3' UTR-längder påverka mRNA-lokalisering, stabilitet och översättning5. Det finns en klass av 3 'slutsekvenseringsanalyser optimerade för att upptäcka APA som skiljer sig åt i detaljerna i protokollet6. Rörledningen som beskrivs här är konstruerad för PolyA-seq, men kan anpassas för andra protokoll enligt beskrivningen.

I denna studie presenterar vi en pipeline av differentiella exonanalysmetoder7,8 (figur 1), som kan delas in i två breda kategorier: exonbaserad (DEXSeq9, diffSplice 10) och händelsebaserad (replikera multivariat analys av transkriptskarvning (rMATS)11). De exonbaserade metoderna jämför vikförändringen över förhållanden för enskilda exoner, mot ett mått på den totala genvecksförändringen för att kalla differentiellt uttryckt exonanvändning, och utifrån det beräkna ett gennivåmått på AS-aktivitet. Händelsebaserade metoder använder exon-intron-spanning junction reads för att upptäcka och klassificera specifika skarvningshändelser, till exempel exon-hoppning eller retention av introner, och särskilja dessa AS-typer i utdata3. Således ger dessa metoder kompletterande vyer för en fullständig analys av AS12,13. Vi valde DEXSeq (baserat på DESeq214 DGE-paketet) och diffSplice (baserat på Limma10 DGE-paketet) för studien eftersom de är bland de mest använda paketen för differentialsplitsningsanalys. rMATS valdes som en populär metod för händelsebaserad analys. En annan populär händelsebaserad metod är MISO (Blandning av isoformer)1. För APA anpassar vi det exobaserade tillvägagångssättet.

Figure 1
Figur 1. Analys pipeline. Flödesschema över stegen som används i analysen. Stegen inkluderar: att hämta data, utföra kvalitetskontroller och läsa justering följt av att räkna läsningar med hjälp av anteckningar för kända exoner, introner och pA-platser, filtrering för att ta bort låga räkningar och normalisering. PolyA-seq-data analyserades för alternativa pA-ställen med diffSplice/DEXSeq-metoder, bulk-RNA-Seq analyserades för alternativ skarvning på exonnivå med diffSplice/DEXseq-metoder och AS-händelser analyserades med rMATS. Klicka här för att se en större version av denna siffra.

RNA-seq-data som används i denna undersökning förvärvades från Gene Expression Omnibus (GEO) (GSE138691)15. Vi använde mus-RNA-seq-data från denna studie med två tillståndsgrupper: vildtyp (WT) och muskelblindliknande typ 1-knockout (Mbnl1 KO) med tre replikat vardera. För att demonstrera differentiell polyadenylationsplatsanvändningsanalys erhöll vi musembryofibroblaster (MEF) PolyA-seq-data (GEO Accession GSE60487)16. Data har fyra villkorsgrupper: Wild-type (WT), Muscleblind-like type1/type 2 double knockout (Mbnl1/2 DKO), Mbnl 1/2 DKO med Mbnl3 knockdown (KD) och Mbnl1/2 DKO med Mbnl3-kontroll (Ctrl). Varje villkorsgrupp består av två replikat.

GEO-anslutning SRA Run-nummer Exempel på namn Tillstånd Replikera Vävnad Sekvensering Läslängd
RNA-Seq GSM4116218 SRR10261601 Mbnl1KO_Thymus_1 Mbnl1 knockout Rep 1 Bräss Parat slut 100 bp
GSM4116219 SRR10261602 Mbnl1KO_Thymus_2 Mbnl1 knockout Rep 2 Bräss Parat slut 100 bp
GSM4116220 SRR10261603 Mbnl1KO_Thymus_3 Mbnl1 knockout Rep 3 Bräss Parat slut 100 bp
GSM4116221 SRR10261604 WT_Thymus_1 Vild typ Rep 1 Bräss Parat slut 100 bp
GSM4116222 SRR10261605 WT_Thymus_2 Vild typ Rep 2 Bräss Parat slut 100 bp
GSM4116223 SRR10261606 WT_Thymus_3 Vild typ Rep 3 Bräss Parat slut 100 bp
3P-Seq GSM1480973 SRR1553129 WT_1 Vild typ (WT) Rep 1 Musembryonala fibroblaster (MEF) Enkel ände 40 punkter
GSM1480974 SRR1553130 WT_2 Vild typ (WT) Rep 2 Musembryonala fibroblaster (MEF) Enkel ände 40 punkter
GSM1480975 SRR1553131 DKO_1 Mbnl 1/2 dubbel knockout (DKO) Rep 1 Musembryonala fibroblaster (MEF) Enkel ände 40 punkter
GSM1480976 SRR1553132 DKO_2 Mbnl 1/2 dubbel knockout (DKO) Rep 2 Musembryonala fibroblaster (MEF) Enkel ände 40 punkter
GSM1480977 SRR1553133 DKOsiRNA_1 Mbnl 1/2 dubbel knockout med Mbnl 3 siRNA (KD) Rep 1 Musembryonala fibroblaster (MEF) Enkel ände 40 punkter
GSM1480978 SRR1553134 DKOsiRNA_2 Mbnl 1/2 dubbel knockout med Mbnl 3 siRNA (KD) Rep 2 Musembryonala fibroblaster (MEF) Enkel ände 36 bp
GSM1480979 SRR1553135 DKONTsiRNA_1 Mbnl 1/2 dubbel knockout med icke-målinriktat siRNA (Ctrl) Rep 1 Musembryonala fibroblaster (MEF) Enkel ände 40 punkter
GSM1480980 SRR1553136 DKONTsiRNA_2 Mbnl 1/2 dubbel knockout med icke-målinriktat siRNA (Ctrl) Rep 2 Musembryonala fibroblaster (MEF) Enkel ände 40 punkter

Tabell 1. Sammanfattning av RNA-Seq- och PolyA-seq-dataset som används för analysen.

Subscription Required. Please recommend JoVE to your librarian.

Protocol

1. Installation av verktyg och R-paket som används i analysen

  1. Conda är en populär och flexibel pakethanterare som möjliggör bekväm installation av paket med deras beroenden på alla plattformar. Använd 'Anaconda' (conda package manager) för att installera 'conda' som kan användas för att installera de verktyg/paket som krävs för analysen.
  2. Ladda ner 'Anaconda' enligt systemkraven från https://www.anaconda.com/products/individual#Downloads och installera det genom att följa anvisningarna i grafiskt installationsprogram. Installera alla nödvändiga paket med "conda" genom att skriva följande på Linux-kommandoraden.
    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. Om du vill ladda ned alla R-paket som används i protokollet skriver du följande kod i R-konsolen (startade på Linux-kommandoraden genom att skriva "R") eller Rstudio-konsolen.
    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)
    }
    OBS: I detta beräkningsprotokoll kommer kommandon att ges som antingen R Notebook-filer (filer med tillägg ". Rmd"), R-kodfiler (filer med tillägg ". R") eller Linux Bash-skalskript (filer med tillägg ".sh"). R Notebook -filer (Rmd) ska öppnas i RStudio med File | Öppna Arkiv ..., och enskilda kodsegment (som kan vara R-kommandon eller Bash-skalkommandon) körs sedan interaktivt genom att klicka på den gröna pilen längst upp till höger. R-kodfiler kan köras genom att öppna i RStudio, eller på Linux-kommandoraden genom att prefacing med "Rscript", t.ex. Rscript example.R. Shell-skript körs på Linux-kommandoraden genom att prefacing skriptet med kommandot "sh" t.ex.sh example.sh.

2. Analys av alternativ skarvning (AS) med RNA-seq

  1. Nedladdning och förbehandling av data
    OBS: Kodavsnitten som kommenteras nedan finns i den kompletterande kodfilen "AS_analysis_RNASeq.Rmd", för att följa de enskilda stegen interaktivt och tillhandahålls också som ett bash-skript som ska köras i batch på Linux-kommandoraden (sh downloading_data_preprocessing.sh).
    1. Ladda ner rådata.
      1. Ladda ned rådata från Sequence Read Archive (SRA) med kommandot 'prefetch' från SRA toolkit (v2.10.8)17. Ge SRA-anslutnings-ID:n i följd i följande kommando för att ladda ner dem parallellt med GNU parallel utility18. Om du vill hämta SRA-filer med anslutnings-ID:n från SRR10261601 till SRR10261606 parallellt använder du följande på Linux-kommandoraden.
        ​seq 10261601 10261606 | parallel prefetch SRR{}
      2. Extrahera fastq-filerna från arkivet med funktionen "fastq-dump" från SRA-verktygslådan. Använd GNU parallellt och ge namnen på alla SRA-filer tillsammans.
        ​parallel -j 3 fastq-dump --gzip --skip-technical --read-filter pass --dumpbase --split-e --clip --origfmt {} :::
      3. Ladda ned referensgenomet och anteckningarna för mus (genomsammansättning GRCm39) från www.ensembl.org med hjälp av följande på Linux-kommandoraden.
        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. Förbehandling och kartläggning läser till genomsammansättningen
      1. Kvalitetskontroll. Bedöm kvaliteten på råläsningar med FASTQC (FASTQ Quality Check v0.11.9)19. Skapa en utdatamapp och kör fastqc med parallell på flera indatafiler. Det här steget genererar en kvalitetsrapport för varje prov. Undersök rapporterna för att säkerställa att läskvaliteten är acceptabel innan du gör ytterligare analys. (Se användarhandboken för att förstå rapporterna på https://www.bioinformatics.babraham.ac.uk/projects/fastqc/)
        mkdir fastqc_out
        parallel "fastqc {} -o fastqc_out" ::: $RAW_DATA/*.fastq.gz
        OBS: Utför vid behov adaptertrimning med 'cutadapt'20 eller 'trimmomatic'21 för att ta bort sekvensering i flankerande adaptrar, som varierar beroende på RNA-fragmentstorlek och läslängd. I den här analysen hoppade vi över det här steget eftersom andelen läsningar som påverkades var minimal.
      2. Läs justering. Nästa steg i förbehandlingen inkluderar att kartlägga läsningarna till referensgenomet. För det första, bygg indexet för referensgenomet med hjälp av funktionen "genomeGenerate" i STAR22 och justera sedan de råa läsningarna till referensen (alternativt finns förbyggda index tillgängliga från STAR-webbplatsen och kan användas direkt för justering). Kör följande kommandon på Linux-kommandoraden.
        #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
        STAR-justeraren genererar och sorterar BAM-filer (Binary Alignment Map) för varje prov efter läsjustering. Bam-filer måste sorteras innan du fortsätter till ytterligare steg.
  2. Förbereda Exon-anteckningar.
    1. Kör den kompletterande kodfilen "prepare_mm_exon_annotation. R" med den nedladdade anteckningen i GTF-format (Gene transfer format) för att förbereda anteckningarna. Om du vill köra skriver du följande på Linux-kommandoraden.
      Rscript prepare_mm_exon_annotation.R annotation.gtf
      OBS: GTF-filen innehåller flera exon-poster för olika isoformer. Den här filen används för att "komprimera" flera transkriptions-ID:n för varje exon. Det är ett viktigt steg att definiera exonräkningsfack.
  3. Räkna läser. Nästa steg är att räkna antalet läsningar som mappats till olika transkriptioner/exoner. Se Kompletterande fil: "AS_analysis_RNASeq.Rmd".
    1. Läs in nödvändiga bibliotek:
      packages<- c("Rsubread","tidyverse", "magrittr", "EnhancedVolcano", "edgeR","openxlsx")
      invisible(lapply(packages, library, character.only=TRUE))
    2. Läs in den bearbetade anteckningsfilen som erhållits från föregående steg (2.2).
      load("mm_exon_anno.RData")
    3. Läs alla bam-filer som erhölls i steg 2.2.2 som indata för 'featureCounts' för att räkna läsningar. Läs mappen som innehåller bam-filer genom att först lista varje fil från katalogen som slutar med .bam. Använd "featureCounts" från Rsubread-paketet som tar bam-filer och bearbetade GTF-anteckningar (referens) som indata för att generera en matris med antal som är associerade med varje funktion med rader som representerar exons(features) och kolumner som representerar exempel.
      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. Utför sedan icke-specifik filtrering för att ta bort lågt uttryckta exoner ("icke-specifik" indikerar att den experimentella tillståndsinformationen inte används i filtreringen för att undvika urvalsfördomar). Transformera data från råskala till antal per miljon (cpm) med hjälp av cpm-funktionen från "edgeR"-paket23 och behåll exons med räkningar som är större än ett inställbart tröskelvärde (för den här datauppsättningen används en cpm) i minst tre exempel. Ta också bort gener med endast 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, ]
      Kontrollera de nödvändiga parametrarna för featureCounts när du använder olika data, till exempel för enkelläsningar, ställ in 'isPairedEnd = FALSE'. Se användarhandboken för RSubread för att välja alternativ för dina data och se avsnittet Diskussion nedan.
  4. Differentiell skarvning och exonanvändningsanalys. Vi beskriver två alternativ för detta steg: DEXSeq och DiffSplice. Endera kan användas och ge liknande resultat. För konsekvens väljer du DEXSeq om du föredrar ett DESeq2-paket för DGE och använder DiffSplice för en Limma-baserad DGE-analys. Se kompletterande fil: "AS_analysis_RNASeq.Rmd".
    1. Använda DEXSeq-paket för differentiell exonanalys.
      1. Läs in biblioteket och skapa en exempeltabell för att definiera experimentdesignen.
        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")))
        Radnamnen ska vara konsekventa med de bam-filnamn som används av featureCounts för att räkna läsningarna. sampleTable består av detaljer om varje prov som innehåller: bibliotekstyp och villkor. Detta krävs för att definiera kontrasterna eller testgruppen för att identifiera differentiell användning.
      2. Förbered exon-informationsfilen. Exoninformation i form av GRanges -objekt (genomiska intervall) (https://bioconductor.org/packages/release/bioc/html/GenomicRanges.html) krävs som indata för att skapa DEXSeq-objektet i nästa steg. Matcha gen-ID:n med antalet läsningar för att skapa 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. Skapa DEXSeq-objekt med funktionen DEXSeqDataSet. DEXSeq-objektet samlar in läsantal, exon-funktionsinformation och exempelinformation. Använd läsantalet som genererades i steg 3 och exoninformationen som erhållits från föregående steg för att skapa DEXSeq-objektet från räkningsmatrisen. Argumentet sampleData tar en datarams indata som definierar exemplen (och deras attribut: bibliotekstyp och villkor), "design" använder sampleData för att generera en designmatris för differentialtestning med hjälp av modellformelnotation. Observera att en signifikant interaktionsterm, condition:exon, indikerar att bråkdelen av läsningar över en gen som faller på en viss exon beror på det experimentella tillståndet, dvs det finns AS. Se DEXSeq-dokumentationen för en fullständig beskrivning av inställningen av modellformeln för mer komplexa experimentella konstruktioner. För funktionsinformation krävs exon-id, motsvarande gen och transkriptioner.
        ​dxd<-DEXSeqDataSet(countData$counts,sampleData=sampleTable, design=~sample+exon+condition:exon,featureID=exoninfo$ExonID,groupID=exoninfo$GeneID,featureRanges=exoninfo, transcripts=transcripts_l)
      4. Normalisering och spridningsuppskattning. Utför sedan normalisering mellan prover och uppskatta variansen av data, på grund av både Poisson-räkningsbrus från den diskreta naturen hos RNA-seq och biologisk variabilitet, med hjälp av följande kommandon.
        ​dxd %<>% estimateSizeFactors %>% estimateDispersions %T>% plotDispEsts
      5. Testa för differentiell användning. Efter uppskattning av variation, testa för differentiell exonanvändning för varje gen och generera resultaten.
        dxd%<>%testForDEU%>%estimateExonFoldChanges(fitExpToVar=
        "condition")#Estimate fold changes
        dxr=DEXSeqResults(dxd)
      6. Visualisera skarvningshändelser för valda gener med hjälp av följande 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" för att generera ytterligare diagram för gener av intresse och för att generera vulkandiagram vid olika trösklar.
    2. Använda diffSplice från Limma för att identifiera differentiell skarvning. Följ R Notebook-filen "AS_analysis_RNASeq.Rmd". Se till att steg 2.1-2.3 har följts för att förbereda indatafiler innan du fortsätter vidare.
      1. Ladda bibliotek
        library(limma)
        library(edgeR)
      2. Icke-specifik filtrering. Extrahera matrisen med läsantal som erhållits i 2.3. Skapa en lista över funktioner med funktionen 'DGEList' från edgeR-paketet, där rader representerar gener och kolumner representerar prover.
        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, "\\,")
        OBS: Som ett icke-specifikt filtreringssteg filtreras antalet med cpm < 1 in x av n prover, där x är det minsta antalet replikat i alla förhållanden. n = 6 och x = 3 för dessa exempeldata.
      3. Normalisera antalet mellan exempel, med funktionen "calcNormFactors" från "edgeR"-paketet med hjälp av trimmade medelvärden av M-värden (TMM-normaliseringsmetod)24 Den beräknar skalningsfaktorer för att justera biblioteksstorlekar.
        ​dge<-calcNormFactors(dge)
      4. Använd sampleTable som genererades i steg 2.4.1.1 och skapa designmatrisen. Designmatrisen kännetecknar designen. Se Limma användarhandbok (https://www.bioconductor.org/packages/devel/bioc/vignettes/limma/inst/doc/usersguide.pdf) kapitel 8 & 9 för detaljer om designmatriser för mer avancerade experimentella mönster.
        Treat<- factor(sampleTable$condition)
        design<- model.matrix(~0+Treat)
        colnames(design) <- levels(Treat)
      5. Montera en linjär modell per exon. Kör "voom" -funktionen i "limma" -paketet för att bearbeta RNA-seq-data för att uppskatta varians och generera precisionsvikter för att korrigera för Poisson-räkningsbrus och omvandla exon-nivåräkningarna till log2-räkningar per miljon (logCPM). Kör sedan linjär modellering med funktionen 'lmfit' för att anpassa linjära modeller till uttrycksdata för varje exon. Beräkna empirisk Bayes-statistik för anpassad modell med hjälp av 'eBayes' -funktionen för att detektera differentiellt exonuttryck. Definiera sedan en kontrastmatris för de experimentella jämförelserna av intresse. Använd 'contrasts.fit' för att få koefficienter och standardfel för varje jämförelsepar.
        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. Differentiell skarvningsanalys. Kör 'diffSplice' på den anpassade modellen för att testa skillnaderna i exonanvändning av gener mellan vildtyp och knockout och utforska de högst rankade resultaten med hjälp av 'topSplice' -funktionen: test = "t" ger en ranking av AS exons, test = "simes" ger en ranking 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. Plotta resultaten med funktionen "plotSplice", vilket ger gen av intresse för geneidargumentet. Spara de bästa resultaten sorterade efter logg Vik ändra till ett objekt och generera ett vulkandiagram för att visa exonerna.
        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. Använda rMATS
      1. Se till att den senaste versionen av rMATS v4.1.1 (även känd som rMATS turbo på grund av den minskade bearbetningstiden och mindre krav på minne) installeras antingen med conda eller github (https://github.com/Xinglab/rmats-turbo/releases/download/v4.1.1/rmats_turbo_v4_1_1.tar.gz) i arbetskatalogen. Följ avsnitt 4.3 i "AS_analysis_RNASeq.Rmd".
      2. Gå till mappen som innehåller bam-filer som erhållits efter mappning och förbered textfiler, enligt kraven i rMATS, för de två villkoren genom att kopiera namnet på bam-filer (tillsammans med sökvägen) åtskilda av ','. Följande kommandon ska köras på Linux-kommandoraden:
        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 två indatafilerna som genererades i föregående steg, tillsammans med GTF-filen som erhölls i 2.1.1.3. Detta genererar en utdatamapp 'rmats_out' som innehåller textfiler som beskriver statistik (p-värden och inklusionsnivåer) för varje skarvningshä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
        OBS: Referensanteckningen i form av en GTF-fil krävs också. Kontrollera parametrarna om data är enkeländade och ändra alternativet -t i enlighet därmed.
      4. Utforska rMATS-resultat. Använd Bioconductor-paketet 'maser'25 för att utforska rMATS-resultaten. Ladda textfilerna Junction och exon counts (JCEC) i "maser" -objektet och filtrera resultatet baserat på täckning genom att inkludera minst fem genomsnittliga läsningar per skarvningshä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 av rMATS-resultat. Välj de signifikanta skarvningshändelserna vid False Discovery Rate (FDR) 10 % och minst 10 % förändring i Procent skarvade in (deltaPSI) med funktionen "topEvents" från "maser"-paketet. Kontrollera sedan genhändelserna för enskilda gener av intresse (provgen-Wnk1) och plotta PSI-värden för varje skarvningshändelse för den genen. Generera ett vulkandiagram genom att ange händelsetypen.
        #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. Generera Sashimi-diagram för skarvningshändelseresultatet som erhållits med rMATS i form av textfiler med hjälp av paketet 'rmats2shahimiplot'. Kör python-skriptet på Linux-kommandoraden.
        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
        OBS: Den här processen kan vara tidskrävande eftersom den genererar Sashimi-diagrammet för alla resultat i händelsefilen. Välj de bästa resultaten (gennamn och exoner) som visas av topEvents-funktionen från 'maser' och visualisera motsvarande Sashimi-diagram.

3. Analys av alternativ polyadenylering (APA) med 3'-slutsekvensering

  1. Nedladdning och förbehandling av data
    OBS: Se den kompletterande R Notebook-filen "APA_analysis_3PSeq_notebook. Rmd" för de fullständiga kommandona för datanedladdning och förbehandlingssteg, eller kör den kompletterande bash-filen "APA_data_downloading_preprocessing.sh" på Linux-kommandoraden.
    1. Ladda ner data från SRA med anslutnings-ID:n (1553129 till 1553136).
    2. Trimadaptrar och omvänd komplement för att få känslasträngsekvensen.
      OBS: Detta steg är specifikt för den PolyA-seq-analys som används.
    3. Karta läser till musgenomsammansättning med bowtie aligner26.
  2. Förbereda anteckningar för pA-webbplatser.
    OBS: Bearbetningen av pA-platsannoteringsfilen utförs först med hjälp av kompletterande R Notebook-fil "APA_analysis_3PSeq_notebook. Rmd" (2.1 - 2.6), och sedan använda bash-filen "APA_annotation_preparation.sh".
    1. Ladda ned pA-webbplatsanteckningar från PolyASite 2.0-databasen6.
    2. Välj pA-platsanteckningar för att behålla 3'-oöversatta region (UTR) pA-platser, som kommenteras som Terminal Exon (TE) eller 1000 nt nedströms om en kommenterad terminal exon (DS) för nedströmsanalys.
    3. Skaffa pA-platstoppar. Ankra vid varje pA-klyvningsplats och visualisera den genomsnittliga lästäckningen med hjälp av sängverktyg och deeptools27,28. Resultaten visade att topparna på de kartlagda avläsningarna huvudsakligen var spridda inom ~60 bp uppströms klyvningsplatserna (figur 5 och kompletterande figur 5). Därför utökades koordinaterna för pA-platser från annoteringsfilen till 60 bp uppströms deras klyvningsplatser. Beroende på det specifika 3'-slutsekvenseringsprotokollet som används måste detta steg optimeras för andra analyser än PolyA-seq.
  3. Räkna läser
    1. Förbered pA-platsanteckningsfilen.
      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. Använd 'featureCounts' för att få råantal. Spara räkningstabellen som RData-filen "APA_countData.Rdata" för APA-analys med olika verktyg.
      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")
      OBS: Var medveten om att ändra någon av parametrarna som anges i funktionen 'featureCounts'. Ändra parametern 'strandSpecific' för att säkerställa att den överensstämmer med sekvenseringsriktningen för den 3'-ändsekvenseringsanalys som används (empiriskt kommer visualisering av data i en genomwebbläsare över gener på plus- och minussträngar att klargöra detta).
    3. Tillämpa icke-specifik filtrering av countData. Filtrering kan avsevärt förbättra den statistiska robustheten i differentiella pA-webbplatsanvändningstester. Först tog vi bort de gener med endast ett pA-ställe, på vilka differentiellt pA-webbplatsanvändning inte kan definieras. För det andra tillämpar vi icke-specifik filtrering baserat på täckning: antalet filtreras med cpm mindre än 1 in x av n prover, där x är det minsta antalet replikat i alla förhållanden. N = 8 och x = 2 för dessa exempeldata.
      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. Differentiell polyadenyleringsplatsanvändningsanalys med hjälp av DEXSeq- och diffSplice-rörledningar.
    1. Använda DEXSeq-paketet
      OBS: Eftersom en kontrastmatris inte kan definieras för DEXSeq-pipelinen måste den differentiella APA-analysen av varje två experimentella förhållanden utföras separat. Den differentiella APA-analysen av tillståndet WT och villkoret DKO utförs som ett exempel för att förklara proceduren. Se tilläggsfilen "APA_analysis_3PSeq_notebook. Rmd" för det stegvisa arbetsflödet i detta avsnitt och den differentiella APA-analysen av andra kontraster.
      1. Läs in biblioteket och skapa en exempeltabell för att definiera experimentdesignen.
        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. Förbered informationsfilen för pA-platser med hjälp av Bioconductor-paketet 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. Använd läsantalet som genererades i steg 3.3 och pA-platsinformationen som hämtades från föregående steg för att skapa 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. Definiera kontrastparet genom att definiera villkorsnivåerna i DEXSeq-objektet.
        dxd1$condition<- factor(dxd1$condition, levels= c("WT", "DKO"))
        # The contrast pair will be "DKO - WT"
      5. Normalisering och spridningsuppskattning. I likhet med RNA-seq-data utför för 3 'slutsekvenseringsdata normalisering mellan prover (kolumnvis median av förhållanden för varje prov) med hjälp av funktionen 'estimateSizeFactors' och uppskattar variationen av data med funktionen 'estimateDispersions' och visualiserar sedan dispersionsuppskattningsresultatet med funktionen 'plotDispEsts'.
        ​dxd1 %<>% estimateSizeFactors %>% estimateDispersions %T>% plotDispEsts
      6. Differentiell pA-platsanvändningstestning för varje gen med funktionen 'testForDEU', uppskatta sedan vikförändringen av pA-platsanvändning med funktionen 'estimateExonFoldChanges'. Kontrollera resultaten med funktionen 'DEXSeqResults' och ställ in 'FDR < 10%' som kriterium för signifikant differentiella pA-platser.
        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 differentiella pA-platsanvändningsresultat med hjälp av differentiella APA-diagram som genereras av funktionen 'plotDEXSeq' och vulkandiagram av 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. Använda diffSplice-paketet. Se den kompletterande R Notebook-filen "APA_analysis_3PSeq_notebook. Rmd" för det stegvisa arbetsflödet i det här avsnittet.
      1. Definiera kontraster av intresse för differentiell pA-användningsanalys.
        OBS: Detta steg bör utföras efter konstruktion och bearbetning av DGEList-objektet, som ingår 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. Visualisera resultatet av kontraster av intresse (här "DKO - WT") med hjälp av differentiella APA-diagram med funktionen 'plotSplice' och vulkandiagram med funktionen 'EnhancedVolcano'. Se R Notebook-filen "APA_analysis_3PSeq_notebook. Rmd" 4.2.7 - 4.2.9 för visualisering av andra 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 att ha kört ovanstående steg-för-steg-arbetsflöde är AS- och APA-analysutdata och representativa resultat i form av tabeller och datadiagram, genererade enligt följande.

SOM:
Huvudresultatet av AS-analysen (tilläggstabell 1 för diffSplice; Tabell 2 för DEXSeq) är en lista över exoner som visar differentiell användning över förhållanden, och en lista över gener som visar signifikant total skarvningsaktivitet för en eller flera av dess beståndsdelar, rangordnade efter statistisk signifikans. Tilläggstabell 1, flik 2 visar signifikanta exoner, med kolumner som visar differentiellt FC av exon kontra vila, per-exon ojusterat p-värde och justerat p-värde (Benjamini-Hockberg-korrigering). Tröskelvärden på de justerade p-värdena ger en uppsättning exoner med definierad FDR. Kompletterande tabell 1, flik 3 visar en rangordnad lista över gener som visar signifikans av den totala skarvningsaktiviteten, med en kolumn som visar gennivåjusterat p-värde beräknat med Simes-metoden. Liknande data visas i tabell 2 för DEXSeq. Kompletterande figur 1 och kompletterande figur 2 visar differentiellt skarvningsmönster i Mbnl1-, Tcf7- och Lef1-gener som har validerats experimentellt i den publicerade artikeln som presenteras med data15. Författarna har visat experimentell validering av fem gener- Mbnl1, Mbnl2, Lef1, Tcf7 och Ncor2. Vårt tillvägagångssätt upptäckte differentiell skarvning i alla dessa gener. Här presenterar vi FDR-nivåerna för varje gen med hjälp av DEXSeq, diffSplice respektive rMATS som erhållits i tilläggstabellerna 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) och Ncor2 (9.2E-11, 1.6E-06, 0).

Figur 2 visar en jämförelse mellan utdata erhållna från tre olika verktyg och illustrerar alternativa skarvningsmönster i Wnk1-genen . Vulkandiagrammen visas i figur 2A (diffSplice) och figur 2B (DEXSeq). Ytterligare tre högt rankade gener visas i kompletterande figur 1 (diffSplice) och kompletterande figur 2 (DEXSeq). Y-axeln visar statistisk signifikans (-log10 P-värden) och x-axeln visar effektstorlek (vikändring). Gener som ligger längst upp till vänster eller höger kvadranter indikerar betydande FC och starka statistiska bevis på verkliga skillnader.

Figure 2
Figur 2. Jämförelse av alternativa skarvningsresultat erhållna från diffSplice, DEXSeq och rMATS. |
(A) Vulkandiagram (vänster) av RNA-Seq från Limma diffSplice-analys: X-axeln visar log exon-vikförändring; Y-axeln visar-log 10 p-värde. Varje punkt motsvarar en exon. Horisontell streckad linje vid p-värde = 1E-5; vertikala streckade linjer vid tvåfaldig förändring (FC). Röda exoner visar betydande FC och statistisk signifikans. Differentiellt skarvningsdiagram (höger): Skarvningsmönster uppvisas för en exempelgen Wnk1 där x-axeln visar exon-id per transkript; y-axeln visar exons relativa loggviktsändring (skillnaden mellan exonens logFC och den totala logFC för alla andra exoner). Exoner markerade med rött visar statistiskt signifikant differentialuttryck (FDR < 0,1).
(B) Vulkandiagram (vänster) och differentiell exonanvändning (höger) av RNA-Seq erhållen från DEXSeq-analys. Wnk1-genen visar signifikant differentiell användning av exoner mellan WT och Mbnl1 knock-out markerad med rosa, vilket motsvarar samma differentiella exoner i (A).
(C) Vulkandiagram (vänster) och Sashimi-plot (höger) för Wnk1 erhållet från rMATS-analys. Vulkandiagram som visar signifikant hoppad (kassett) exon (SE) händelse i vildtyp jämfört med knockout vid 10% FDR med förändring i procent skarvade i (PSI eller ΔΨ) värden > 0,1. X-axeln visar förändringar i PSI-värden över villkor och y-axeln visar logg-P-värde. Sashimi-diagrammet visar en hoppad exonhändelse i Wnk1-genen , vilket motsvarar en signifikant differentiell exon i (A) och (B). Varje rad representerar ett RNA-Seq-prov: tre replikat av vildtyp och Mbnl1-knock-out. Höjden visar lästäckning i RPKM och de anslutande bågarna visar korsning läser över exoner. Kommenterade genmodellalternativ isoformer visas längst ner i diagrammet. Den nedre panelen i C illustrerar korsningsläsningar som används för att beräkna PSI-statistiken.
(D) Venndiagram som jämför antalet signifikanta differentiella exoner erhållna med de olika metoderna. Klicka här för att se en större version av denna siffra.

Figur 2 A (högra panelen) visar en diagrammatisk visning av exonskillnader i en av de högst rankade generna, som visar logFC på y-axeln och exon-numret på x-axeln. Detta exempel visar en kassettexon som varierar mellan förhållanden för gen Wnk1. Den differentiella exonanvändningsdiagrammet från DEXSeq visar bevis på differentiell skarvning vid fem exonplatser nära Wnk1.6.45. De markerade exonerna i rosa kommer sannolikt att skarvas ut i Mbnl1 KO-prover jämfört med WT. Dessa exoner kompletterar resultaten erhållna av diffSplice som visar ett liknande mönster vid den specifika genomiska positionen. Fler exempel visas i kompletterande figur 1 och kompletterande figur 2. En mer detaljerad vy för att bekräfta intressanta resultat kan ges genom att jämföra täckningsspår (vicka) i RPM -enheter (Läser per miljon) i UCSC (University of Santa Cruz) eller IGV (Integrative Genomics Viewer) genomwebbläsare (visas inte), tillsammans med visuell korrelation med andra spår av intresse, såsom kända genmodeller, bevarande och andra genomomfattande analyser.

rMATS-utdatatabellen visar signifikanta alternativa skarvningshändelser kategoriserade efter typ (tilläggstabell 3). Figur 2C visar ett vulkandiagram över gener som är alternativa skarvade, med effektstorleken mätt med differentialstatistiken "procent skarvad i" (PSI eller ΔΨ) på11.

PSI avser procentandelen läsningar som överensstämmer med införandet av en kassettexon (dvs. läser mappning till själva kassettexonen eller korsningsläsningar som överlappar exon) jämfört med läsningar som överensstämmer med exon-uteslutning, dvs. korsning läser över intilliggande uppströms och nedströms exoner (Den nedre panelen i figur 2C). Den högra panelen i figur 2C visar sashimi-plot av Wnk1-genen med differentiell skarvningshändelse överlagrad på täckningsspår för genen, med en hoppad exon i Mbnl1 KO. Bågar som går med i exoner visar antalet korsningsläsningar (läser som korsar en utskuren intron). Olika flikar i tilläggstabell 3 visar signifikanta läsningar av varje typ av händelse som sträcker sig över korsningar med exongränser (korsningsantal och exonantal (JCEC)). Figur 2D jämför de signifikanta differentiellt skarvade exoner som detekteras av de tre verktygen.

Figure 3
Figur 3. Alternativa skarvningshändelser förvärvade genom rMATS-analys. A) Typer av AS-händelser. Denna figur är anpassad från rMATS dokumentation11 som förklarar skarvningshändelsen med konstitutiva och alternativt skarvade exoner. Märkt med antalet för varje händelse vid FDR 10%. B) Sashimi-plot av Add3-gen som uppvisar hoppad exon (SE). C) Sashimi-plot av Baz2b-genen som uppvisar alternativ 5'-skarvplats (A5SS). D) Sashimi-plot av Lsm14b-gen som uppvisar alternativ 3'-skarvplats (A3SS). E) Sashimi-plot av Mta1-gen som uppvisar ömsesidigt exklusiva exoner (MXE). F) Sashimi-plot av Arpp21-gen som uppvisar kvarhållen intron (RI). Röda rader representerar tre replikat av vild typ och orange rader representerar Mbnl1-knock-out-repliker. X-axeln motsvarar genomiska koordinater och stränginformation, y-axeln visar täckning i RPKM. Klicka här för att se en större version av denna siffra.

Figur 3 illustrerar typer av skarvningshändelser SE, A5SS, A3SS, MXE och RI med hjälp av Sashimi-diagram över de viktigaste generna för dessa händelser. Vid en jämförelse av de tre replikaten av både WT och Mbnl1_KO upptäcktes totalt 1272 SE-händelser, 130 A5SS, 116 A3SS, 215 MXE och 313 RI-händelser vid FDR 10%. Sashimi-plot illustrerar typen av händelse med toppgener som exempel.

APA:
Utdata från APA-analysen liknar AS-analysen på exon-nivå. En tabell över de vanligaste generna rangordnade efter differentiell APA-aktivitet i 3'UTR finns (tilläggstabell 4 och tilläggstabell 5). Figur 4A visar vulkandiagrammen för gener genom differentiell APA-aktivitet i 3'UTR som genereras med diffSplice och DEXSeq separat. Figur 4B visar Venn-diagrammet som jämför de signifikant differentiella pA-platsanvändningsresultaten som erhållits från olika rörledningar. Figur 4C och 4D visar den diagrammatiska representationen av differentiell pA-platsanvändning i 3'UTR för genen Fosl2 (figur 4C) och Papola (figur 4D) genererad med både diffSplice och DEXSeq, som är experimentellt validerade för att visa signifikant distal till proximal förskjutning (Fosl2) och signifikant proximal till distal förskjutning (Papola) av pA-platsanvändning i DKO, respektive. Fler exempel visas i kompletterande figur 3 och kompletterande figur 4.

Figure 4
Figur 4. Alternativa polyadenylationsdiagram med diffSplice och DEXSeq. A) Vulkandiagram över PolyA-seq-data som genereras med diffSplice och DEXSeq. X-axeln visar logg pA platsveckändring; y-axeln visar -log10 p-värde. Varje punkt motsvarar en pA-plats. Horisontell streckad linje vid p-värde = 1E-5; vertikala streckade linjer vid 2-faldig FC. Röda exoner visar betydande FC och statistisk signifikans. B) Venndiagram som jämför de signifikanta differentiella pA-platsanvändningsresultaten som erhållits från olika rörledningar. C-D) Differentiella APA-diagram genererade med diffSplice och DEXSeq som visar de proximala, interna och distala pA-platserna för Fosl2 - och Papola-genen . Siffror genereras med samma funktion som figur 2 (B) men med pA-platser som ersätter exoner. Klicka här för att se en större version av denna siffra.

Figur 5 är ett diagnostiskt diagram för att bekräfta den förväntade läsfördelningen runt kommenterade pA-klyvningsställen för den PolyA-seq-analys som används. Den visar den genomsnittliga täckningen i flankerande regioner förankrade vid kända pA-klyvningsställen på genomomfattande nivå. I det här fallet visualiseras den förväntade högen av läsningar uppströms om webbplatserna. De läsfördelningar som är förankrade på pA-platser för alla PolyA-seq-prover visas i tilläggsfigur 5.

Figure 5
Figur 5. Genomsnittlig täckningsplot runt pA-klyvningsplatser. Klyvningsstället för PolyA-seq-data visas med lodrät streckad linje. X-axeln visar baspositionen i förhållande till pA-klyvningsställen, upp till 100 nukleotider uppströms och nedströms; y-axeln visar den genomsnittliga lästäckningen över alla pA-klyvningsplatser, normaliserad efter biblioteksstorlek i CPM. Klicka här för att se en större version av denna siffra.

De differentiella APA-resultaten av olika kontraster som genereras av samma pipeline kan jämföras och verifieras genom att visualisera lästäckningen för representativa signifikant differentiella pA-platser i genomwebbläsaren. Figur 6A är Venn-diagrammet som jämför den signifikant differentiella pA-platsanvändningen av olika kontraster som förvärvats från diffSplice. Figur 6B-D är IGV-ögonblicksbilderna av lästäckningen vid pA-platser för olika gener, som visar mönstren som överensstämmer med de som upptäcktes i APA-analysen med diffSplice. Figur 6B validerar den signifikanta proximala till distala förskjutningen av pA-platsanvändning för gen Paip2, som unikt detekteras i kontrasten DKO vs WT, men inte i andra två kontraster KD vs WT, och Ctr vs WT. Figur 6C validerar den signifikanta distala till proximala förskjutningen av pA-platsanvändning för gen Ccl2 detekterad unikt i kontrasten KD vs WT, medan figur 6D validerar den signifikanta differentiella pA-användningen av alla kontraster för genen Cacna2d1.

Figure 6
Figur 6. Kontrastjämförelse och verifiering av diffSplice-resultat. A) Venndiagram som jämför signifikanta differentiella pA-platsanvändningsresultat av olika kontraster som förvärvats från diffSplice. B-D) IGV-ögonblicksbild som visualiserar pA toppar täckningen av gener Paip2, Ccl2 och Cacna2d1 över förhållanden. Varje spår representerar lästäckningen i ett specifikt villkor. Klicka här för att se en större version av denna siffra.

Kompletterande figur 1. RNA-Seq-analys av differentialsplitsning med Limma diffSplice. (A) Vulkandiagram över RNA-Seq från Limma diffSplice-analys: X-axeln visar log exon-vikförändring; Y-axeln visar-log 10 p-värde. Varje punkt motsvarar en exon. Horisontell streckad linje vid p-värde = 1E-5; vertikala streckade linjer vid tvåfaldig förändring (FC). Röda exoner visar betydande FC och statistisk signifikans. (B-D) Differentiella skarvningsytor: Skarvningsmönster uppvisas till exempel generna Mbnl1, Tcf7 respektive Lef1, där x-axeln visar exon-id per transkript; y-axeln visar exons relativa loggviktsändring (skillnaden mellan exonens logFC och den totala logFC för alla andra exoner). Exoner markerade med rött visar statistiskt signifikant differentialuttryck (FDR < 0,1). Klicka här för att ladda ner den här filen.

Kompletterande figur 2. RNA-Seq-analys av differentiell exonanvändning med DEXSeq. (A) Vulkan tomt. (B-D) Differentiell exonanvändning av RNA-Seq erhållen från DEXSeq-analys. Generna Mbnl1, Tcf7 respektive Lef1 visar signifikant differentiell användning av exoner mellan WT och Mbnl1 knock-out markerade med rosa, vilket motsvarar samma differentiella exoner i kompletterande figur 1. Klicka här för att ladda ner den här filen.  

Kompletterande figur 3. Alternativa polyadenylationsplottar av diffSplice. A) Vulkandiagram av PolyA-seq-data som genereras med diffSplice i tre kontrastpar erhållna från musen PolyA-seq-data, inklusive dubbel knockout (DKO) vs vildtyp (WT), knock-down (KD) vs WT och kontroll (Ctrl) vs WT. X-axeln visar log pA platsveckförändring; y-axeln visar -log10 p-värde. Varje punkt motsvarar en pA-plats. Horisontell streckad linje vid p-värde = 1E-5; vertikala streckade linjer vid 2-faldig FC. Röda pA-platser visar betydande FC och statistisk signifikans. B) Differentiella APA-diagram genererade med diffSplice som visar de proximala, interna och distala pA-platserna för de högt rankade generna S100a7a, Tpm1 och Smc6. Klicka här för att ladda ner den här filen.  

Kompletterande figur 4. Differentiell pA-användningsanalys av DEXSeq-pipeline. A) Vulkandiagram av PolyA-seq-data som genereras med DEXSeq i tre par erhållna från musen PolyA-seq-data, inklusive dubbel knockout (DKO) vs vildtyp (WT), knock-down (KD) vs WT och kontroll (Ctrl) vs WT. X-axeln visar log pA-platsveckändring; y-axeln visar -log10 p-värde. Varje punkt motsvarar en pA-plats. Horisontell streckad linje vid p-värde = 1E-5; vertikala streckade linjer vid 2-faldig FC. Röda pA-platser visar betydande FC och statistisk signifikans. B) Differentiella APA-diagram som genereras med hjälp av DEXSeq som visar de proximala, interna och distala pA-platserna för de högt rankade generna S100a7a, Tpm1 och Smc6. Klicka här för att ladda ner den här filen.  

Kompletterande figur 5. Genomsnittlig täckningsplot och värmekartor runt pA-klyvningsplatser.  Täckning visas för fyra förhållanden, med gener på framåt- och bakåtsträngar som visas separat. X-axeln visar baspositionen i förhållande till pA-klyvningsställen, upp till 100 nukleotider uppströms och nedströms; y-axeln avser medeltäckningen vid motsvarande relativa baspositioner över alla pA-klyvningsställen. Värmekartor ger en alternativ vy, där varje pA-klyvningsplats visas som en separat rad på x-axeln, ordnad efter täckning. Intensitet visar lästäckning (se förklaring). Klicka här för att ladda ner den här filen.

Kompletterande tabell 1. diffSplice-utdata från AS-analysen. Den första fliken definierar kolumnnamnen för de ursprungliga utdata som presenteras på andra (exon-nivå) och tredje (gennivå) flikar. Klicka här för att ladda ner den här tabellen.

Kompletterande tabell 2. DEXSeq-utdata från AS-analysen. Första fliken definierar kolumnnamnen för de ursprungliga utdata som presenteras på andra (exon-nivå) och tredje (gennivå) flikar. Klicka här för att ladda ner den här tabellen.

Kompletterande tabell 3. rMATS-utdata från AS-analysen. Första fliken definierar kolumnnamnen för sammanfattningsfilen (flik 2) och JCEC-filerna för varje händelse (flik 3-7). Klicka här för att ladda ner den här tabellen.

Kompletterande tabell 4. diffSplice-utdata från APA-analysen. Första fliken definierar kolumnnamnen för de ursprungliga utdata som presenteras på andra (pA-platsnivå) och tredje (gennivå) flikar. Klicka här för att ladda ner den här tabellen.

Kompletterande tabell 5. DEXSeq-utdata från APA-analysen. Första fliken definierar kolumnnamnen för de ursprungliga utdata som presenteras på andra (pA-platsnivå) och tredje (gennivå) flikar. Klicka här för att ladda ner den här tabellen.

Kompletterande tabell 6. En sammanfattning av antalet väsentligt förändrade exoner för AS- och pA-platser för APA. Klicka här för att ladda ner den här tabellen.

Kompletterande tabell 7. En sammanfattning av verktyg och paket som används i AS/APA-analysen. Klicka här för att ladda ner den här tabellen.

Subscription Required. Please recommend JoVE to your librarian.

Discussion

I denna studie utvärderade vi exonbaserade och händelsebaserade metoder för att upptäcka AS och APA i bulk RNA-Seq och 3 'slutsekvenseringsdata. De exonbaserade AS-metoderna ger både en lista över differentiellt uttryckta exoner och en rangordning på gennivå ordnad efter den statistiska signifikansen av den totala differentialsplitsaktiviteten på gennivå (tabellerna 1-2, 4-5). För diffSplice-paketet bestäms differentiell användning genom att anpassa viktade linjära modeller på en exon-nivå för att uppskatta den differentiella log-vikningsförändringen för en exon mot den genomsnittliga log-vikningsförändringen för de andra exonerna inom samma gen (kallad per exon FC). Den statistiska signifikansen på gennivå beräknas genom att kombinera individuella signifikanstester på exonnivå till ett genmässigt test med Simes metod10. Denna rangordning efter differentiell skarvningsaktivitet på gennivå kan därefter användas för att utföra en genuppsättningsberikningsanalys (GSEA) av viktiga vägar som är involverade10. DEXSeq använder en liknande strategi genom att anpassa en generaliserad linjär modell för att mäta differentiell exonanvändning, men skiljer sig åt i vissa steg som filtrering, normalisering och dispersionsuppskattning. När vi jämförde de 500 högst rankade exonerna som visade AS-aktivitet och APA med DEXSeq och DiffSplice, fann vi en överlappning av 310 exoner respektive 300 pA-platser, vilket visar överensstämmelsen mellan de två exonbaserade metoderna, vilket också visades i en tidigare studie 29. Vi rekommenderar att du använder en kombination av både en exon-baserad (antingen DEXSeq eller diffSplice) och en händelsebaserad metod för omfattande detektering och klassificering av AS. För APA kan användare välja antingen DEXSeq eller diffSplice: båda metoderna har visat sig fungera bra över ett brett spektrum av transkriptomikexperiment29.

När man förbereder RNA-seq-biblioteket för en AS-analys är det viktigt att använda ett strängspecifikt bulk-RNA-seq-protokoll8, eftersom en stor del av generna i ryggradsdjurens genom överlappar varandra på olika strängar, och ett icke-strängspecifikt protokoll kan inte skilja dessa överlappande regioner, vilket förvirrar slutlig exondetektering. Ett annat övervägande är läsdjup, med skarvningsanalyser som kräver djupare sekvensering än DGE, t.ex. 30-60 miljoner läsningar per prov, jämfört med 5-25 miljoner läsningar per prov för DGE (https://sapac.support.illumina.com/bulletins/2017/04/considerations-for-rna-seq-read-length-and-coverage-.html). Alla verktyg som demonstreras i protokollet stöder både single-end och paired-end sekvenseringsdata. Om endast kända genanteckningar används för att detektera korsningsläsningar kan enstaka kortare läsningar (≥ 50 bp) användas, även om de novo-detektering av nya skarvkorsningar drar nytta av parad ände och längre (≥ 100bp) läser30,31. Valet av RNA-extraktionsprotokoll - antingen polyA-selektion eller rRNA-utarmning - beror på kvaliteten på RNA och den experimentella frågan - se31 för en diskussion. Beroende på detaljerna i bibliotekskonstruktionen krävs ändringar av exempelskripten som ges här för parametrarna för läsjustering, funktionsräkning och rMATS. Vid beräkning av de initiala läsräkningarna på exon-nivå med featureCounts eller liknande metoder måste man se till att konfigurera funktionsalternativen korrekt för räkningar och stränglighet: i featureCounts ställer vi in argumentet "strandSpecific" på lämpligt sätt för det strängspecifika RNA-seq-protokollet som används; och för kvantifiering på exon-nivå förväntas en läsning mappas över intilliggande exoner, så vi ställer in parametern allowMultiOverlap på TRUE. För APA finns det olika 3'-slutsekvenseringsprotokoll6 som varierar i den exakta placeringen av toppar i förhållande till pA-platsen. För våra exempeldata bestämmer vi att toppen är 60 bp uppströms pA-platsen som visas i figur 5, och denna analys måste anpassas för andra 3'-slutsekvenseringsprotokoll.

I detta protokoll begränsar vi omfattningen till diskussionen om differentialanalyser på nivån för enskilda exoner och skarvningshändelser som består av intilliggande exon-intronkombinationer. Vi diskuterar inte klassen av analyser baserade på isoform de novo-rekonstruktion såsom manschettknappar, Cuffdiff32, RSEM 33, Kallisto34 som syftar till att detektera och kvantifiera det absoluta och relativa uttrycket av hela alternativa isoformer. Exon- och händelsebaserade metoder är känsligare för att detektera enskilda skarvningshändelser30 och ger i många fall all information som behövs för vidare analys, utan behov av kvantifiering på isoformnivå.

Den senaste versionen av källfilerna i detta protokoll finns på https://github.com/jiayuwen/AS_APA_JoVE

Subscription Required. Please recommend JoVE to your librarian.

Disclosures

Författarna har inget att avslöja.

Acknowledgments

Denna studie stöddes av ett Australian Research Council (ARC) Future Fellowship (FT16010043) och 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 utgåva 172
Identifiering av alternativ skarvning och 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