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

Samenvatting

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

Samenvatting

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

Inleiding

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

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

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

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

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

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

Protocol

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

1. Conda is een populaire en flexibele pakketbeheerder die een gemakkelijke installatie van pakketten met hun afhankelijkheden op alle platforms mogelijk maakt. Gebruik 'Anaconda' (conda package manager) om 'conda' te installeren, wat kan worden gebruikt om de tools/pakketten te installeren die nodig zijn voor de analyse.
2. Download 'Anaconda' volgens de systeemvereisten van https://www.anaconda.com/products/individual#Downloads en installeer het door de aanwijzingen in het grafische installatieprogramma te volgen. Installeer alle vereiste pakketten met behulp van 'conda' door het volgende te typen op de Linux-opdrachtregel.
   conda install -c daler sratoolkit
conda install -c conda-forge parallel
conda install -c bioconda star bowtie fastqc rmats rmats2sashimiplot samtools fasterq-dump cutadapt bedtools deeptools
3. Als u alle R-pakketten wilt downloaden die in het protocol worden gebruikt, typt u de volgende code in de R-console (gestart op de Linux-opdrachtregel door 'R') of Rstudio-console te typen.
   bioc_packages<- c("DEXSeq", "Rsubread", "EnhancedVolcano", "edgeR", "limma", "maser","GenomicRanges")
packages<- c("magrittr", "rtracklayer", "tidyverse", "openxlsx", "BiocManager")
#Install if not already installed
installed_packages<-packages%in% rownames(installed.packages())
installed_bioc_packages<-bioc_packages%in% rownames(installed.packages())
if(any(installed_packages==FALSE)) {
install.packages(packages[!installed_packages],dependencies=TRUE)
BiocManager::install(packages[!installed_bioc_packages], dependencies=TRUE)
}
   OPMERKING: In dit computationele protocol worden opdrachten gegeven als R Notebook-bestanden (bestanden met de extensie ". Rmd"), R code bestanden (bestanden met de extensie ". R"), of Linux Bash shell scripts (bestanden met de extensie ".sh"). R Notebook (Rmd)-bestanden moeten in RStudio worden geopend met File| Open Bestand... en afzonderlijke codebrokken (wat R-opdrachten of Bash-shellopdrachten kunnen zijn) en voer vervolgens interactief uit door op de groene pijl rechtsboven te klikken. R-codebestanden kunnen worden uitgevoerd door te openen in RStudio, of op de Linux-opdrachtregel door vooraf te gaan met "Rscript", bijvoorbeeld Rscript example.R. Shell-scripts worden uitgevoerd op de Linux-opdrachtregel door het script vooraf te maken met het commando "sh", bijvoorbeeld .sh example.sh.

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

1. Gegevens downloaden en voorbewerking
   OPMERKING: De onderstaande codefragmenten zijn beschikbaar in het aanvullende codebestand "AS_analysis_RNASeq.Rmd", om de afzonderlijke stappen interactief te volgen, en worden ook geleverd als een bash-script om in batch op de Linux-opdrachtregel te worden uitgevoerd (SH downloading_data_preprocessing.sh).
   1. Het downloaden van de onbewerkte gegevens.
      1. Download de onbewerkte gegevens uit Sequence Read Archive (SRA) met de opdracht 'prefetch' uit SRA toolkit (v2.10.8)¹⁷. Geef de SRA-toetredings-ID's in volgorde in het volgende commando om ze parallel te downloaden met behulp van GNU parallel utility¹⁸. Als u SRA-bestanden van toetredings-id's van SRR10261601 naar SRR10261606 parallel wilt downloaden, gebruikt u het volgende op de Linux-opdrachtregel.
         ​seq 10261601 10261606 | parallel prefetch SRR{}
      2. Pak de fastq-bestanden uit het archief met behulp van de 'fastq-dump'-functie uit de SRA-toolkit. Gebruik GNU parallel en geef de namen van alle SRA-bestanden samen.
         ​parallel -j 3 fastq-dump --gzip --skip-technical --read-filter pass --dumpbase --split-e --clip --origfmt {} :::
      3. Download het referentiegenoom en de annotaties voor muis (Genome assembly GRCm39) van www.ensembl.org met behulp van het volgende op de Linux-opdrachtregel.
         wget -nv -O annotation.gtf.gz http://ftp.ensembl.org/pub/release-103/gtf/mus_musculus/Mus_musculus.GRCm39.103.gtf.gz \ && gunzip -f annotation.gtf.gz
wget -nv -O genome.fa.gz http://ftp.ensembl.org/pub/release-103/fasta/mus_musculus/dna/Mus_musculus.GRCm39.dna.primary_assembly.fa.gz \ && gunzip -f genome.fa.gz
GTF=$(readlink -f annotation.gtf)
GENOME=$(readlink -f genome.fa)
   2. Voorbewerking en mapping leest voor aan de genoomassemblage
      1. Kwaliteitscontrole. Beoordeel de kwaliteit van onbewerkte reads met FASTQC (FASTQ Quality Check v0.11.9)¹⁹. Maak een uitvoermap en voer fastqc parallel uit op meerdere invoer fasta-bestanden. Deze stap genereert een kwaliteitsrapport voor elk monster. Onderzoek de rapporten om ervoor te zorgen dat de kwaliteit van de gelezen teksten acceptabel is voordat u verdere analyse uitvoert. (Raadpleeg de gebruikershandleiding voor meer informatie over de rapporten op https://www.bioinformatics.babraham.ac.uk/projects/fastqc/)
         mkdir fastqc_out
parallel "fastqc {} -o fastqc_out" ::: $RAW_DATA/*.fastq.gz
         OPMERKING: Voer indien nodig adapter trimming uit met '^cutadapt'20 of 'trimmomatic'21 om sequencing in flankerende adapters te verwijderen, wat varieert op basis van de grootte van RNA-fragmenten en de leeslengte. In deze analyse hebben we deze stap overgeslagen omdat de fractie van de getroffen reads minimaal was.
      2. Lees uitlijning. De volgende stap in de voorbewerking omvat het in kaart brengen van de reads aan het referentiegenoom. Bouw eerst de index voor het referentiegenoom met behulp van de 'genomeGenerate'-functie van STAR²² en stem vervolgens de onbewerkte reads af op de referentie (als alternatief zijn vooraf gebouwde indexen beschikbaar op de STAR-website en kunnen rechtstreeks worden gebruikt voor uitlijning). Voer de volgende opdrachten uit op de Linux-opdrachtregel.
         #Build STAR index
GDIR=STAR_indices
mkdir $GDIR
STAR --runMode genomeGenerate --genomeFastaFiles $GENOME --sjdbGTFfile $GTF --runThreadN 8 --genomeDir $GDIR
ODIR=results/mapping
mkdir -p $ODIR
#Align reads to the genome
for fq1 in $RAW_DATA/*R1.fastq.gz;
do
fq2=$(echo $fq1| sed 's/1.fastq.gz/2.fastq.gz/g');
OUTPUT=$(basename ${fq1}| sed 's/R1.fastq.gz//g');
STAR --genomeDir $GDIR \
--runThreadN 12 \
--readFilesCommand zcat \
--readFilesIn ${fq1}${fq2}\
--outFileNamePrefix $ODIR\/${OUTPUT} \
--outSAMtype BAM SortedByCoordinate \
--outSAMunmapped Within \
--outSAMattributes Standard
Done
         OPMERKING: De STAR-aligner genereert en sorteert BAM-bestanden (Binary Alignment Map) voor elk voorbeeld na uitlijning. Bam-bestanden moeten worden gesorteerd voordat ze doorgaan met verdere stappen.
2. Exon-annotaties voorbereiden.
   1. Voer het aanvullende codebestand "prepare_mm_exon_annotation. R" met de gedownloade annotatie in GTF-indeling (Gene transfer format) om de annotaties voor te bereiden. Als u wilt uitvoeren, typt u het volgende op de Linux-opdrachtregel.
      Rscript prepare_mm_exon_annotation.R annotation.gtf
      OPMERKING: Het GTF-bestand bevat meerdere exon-vermeldingen voor verschillende isovormen. Dit bestand wordt gebruikt om de meerdere transcript-ID's voor elke exon "samen te vouwen". Het is een belangrijke stap om exon telbakken te definiëren.
3. Tellen Leest. De volgende stap is het tellen van het aantal reads dat is toegewezen aan verschillende transcripten / exonen. Zie Aanvullend bestand: "AS_analysis_RNASeq.Rmd".
   1. Laad vereiste bibliotheken:
      packages<- c("Rsubread","tidyverse", "magrittr", "EnhancedVolcano", "edgeR","openxlsx")
invisible(lapply(packages, library, character.only=TRUE))
   2. Laad het verwerkte annotatiebestand dat is verkregen uit de vorige stap (2.2).
      load("mm_exon_anno.RData")
   3. Lees alle bam-bestanden verkregen in stap 2.2.2 als invoer voor 'featureCounts' om reads te tellen. Lees de map met bam-bestanden door eerst elk bestand uit de map te vermelden dat eindigt op .bam. Gebruik 'featureCounts' uit het Rsubread-pakket dat bam-bestanden en verwerkte GTF-annotatie (referentie) als invoer gebruikt om een matrix van tellingen te genereren die aan elke functie zijn gekoppeld met rijen die exonen (functies) vertegenwoordigen en kolommen die voorbeelden vertegenwoordigen.
      countData <- dir("bams", pattern=".bam$", full.names=T) %>%
featureCounts(annot.ext=anno,
isGTFAnnotationFile=FALSE,
minMQS=0,useMetaFeatures=FALSE,
allowMultiOverlap=TRUE,
largestOverlap=TRUE,
countMultiMappingReads=FALSE,
primaryOnly=TRUE,
isPairedEnd=TRUE,
nthreads=12)
   4. Voer vervolgens niet-specifieke filtering uit om laag uitgedrukte exonen te verwijderen ("niet-specifiek" geeft aan dat de informatie over de experimentele toestand niet wordt gebruikt in de filtering, om selectiebias te voorkomen). Transformeer de gegevens van ruwe schaal naar tellingen per miljoen (cpm) met behulp van de cpm-functie uit 'edgeR'-pakket²³ en bewaar exonen met tellingen groter dan een instelbare drempel (voor deze dataset wordt één cpm gebruikt) in ten minste drie steekproeven. Verwijder ook genen met slechts één exon.
      # Non-specific filtering: Remove the exons with low counts
isexpr<- rownames(countData$counts)[rowSums(cpm(countData$counts)>1) >=3]
countData$counts<-countData$counts[rownames(countData$counts) %in%isexpr, ]
anno<-anno%>% filter(GeneID%in% rownames(countData$counts))
# Remove genes with only 1 site and NA in geneIDs
dn<-anno%>%group_by(GeneID)%>%summarise(nsites=n())%>% filter(nsites>1&!is.na(GeneID))
anno<-anno%>% filter(GeneID%in%dn$GeneID)
countData$counts<-countData$counts[rownames(countData$counts) %in%anno$GeneID, ]
      OPMERKING: Controleer de vereiste parameters voor featureCounts bij het gebruik van verschillende gegevens, bijvoorbeeld voor single-end reads, stel 'isPairedEnd = FALSE' in. Raadpleeg de RSubread-gebruikershandleiding om de opties voor uw gegevens te kiezen en raadpleeg het gedeelte Discussie hieronder.
4. Differentiële splicing en exon gebruiksanalyse. We beschrijven twee alternatieven voor deze stap: DEXSeq en DiffSplice. Beide kunnen worden gebruikt en vergelijkbare resultaten geven. Voor consistentie selecteert u DEXSeq als u de voorkeur geeft aan een DESeq2-pakket voor DGE en DiffSplice gebruikt voor een op Limma gebaseerde DGE-analyse. Zie Aanvullend dossier: "AS_analysis_RNASeq.Rmd".
   1. Het DEXSeq-pakket gebruiken voor differentiële exonanalyse.
      1. Laad de bibliotheek en maak een voorbeeldtabel om het experimentele ontwerp te definiëren.
         library(DEXSeq)
sampleTable<-data.frame(row.names= c("Mbnl1KO_Thymus_1", "Mbnl1KO_Thymus_2", "Mbnl1KO_Thymus_3", "WT_Thymus_1", "WT_Thymus_2", "WT_Thymus_3"), condition= rep(c("Mbnl1_KO", "WT"),c(3,3)), libType= rep(c("paired-end")))
         OPMERKING: De rijnamen moeten consistent zijn met de bam-bestandsnamen die door featureCounts worden gebruikt om de reads te tellen. sampleTable bestaat uit details van elk voorbeeld, waaronder: bibliotheektype en conditie. Dit is nodig om de contrasten of testgroep te definiëren voor het detecteren van differentieel gebruik.
      2. Bereid het exon-informatiebestand voor. Exon-informatie in de vorm van GRanges-objecten (genomische bereiken) (https://bioconductor.org/packages/release/bioc/html/GenomicRanges.html) is vereist als invoer om het DEXSeq-object in de volgende stap te maken. Vergelijk het gen Ids met de gelezen tellingen om een exoninfo-object te maken.
         exoninfo<-anno[anno$GeneID%in% rownames(countData$counts),]
exoninfo<-GRanges(seqnames=anno$Chr,
ranges=IRanges(start=anno$Start, end=anno$End, width=anno$Width),strand=Rle(anno$Strand))
mcols(exoninfo)$TranscriptIDs<-anno$TranscriptIDs
mcols(exoninfo)$Ticker<-anno$Ticker
mcols(exoninfo)$ExonID<-anno$ExonID
mcols(exoninfo)$n<-anno$n
mcols(exoninfo)$GeneID<-anno$GeneID
transcripts_l= strsplit(exoninfo$TranscriptIDs, "\\,")
save(countData, sampleTable, exoninfo, transcripts_l, file="AS_countdata.RData")
      3. Maak een DEXSeq-object met de functie DEXSeqDataSet. Het DEXSeq-object verzamelt leestellingen, exon-functiegegevens en voorbeeldinformatie. Gebruik de leestellingen die in stap 3 zijn gegenereerd en de exon-informatie die is verkregen uit de vorige stap om het DEXSeq-object te maken op basis van de tellingsmatrix. Het argument sampleData neemt een gegevensframe-invoer die de voorbeelden definieert (en hun kenmerken: bibliotheektype en -voorwaarde), 'ontwerp' gebruikt sampleData om een ontwerpmatrix voor de differentiële tests te genereren met behulp van modelformulenotatie. Merk op dat een significante interactieterm, condition:exon, aangeeft dat de fractie van reads over een gen dat op een bepaald exon valt, afhankelijk is van de experimentele conditie, d.w.z. er is AS. Zie de DEXSeq-documentatie voor een volledige beschrijving van het instellen van de modelformule voor complexere experimentele ontwerpen. Voor functie-informatie zijn exon-id's, overeenkomstige genen en transcripten vereist.
         ​dxd<-DEXSeqDataSet(countData$counts,sampleData=sampleTable, design=~sample+exon+condition:exon,featureID=exoninfo$ExonID,groupID=exoninfo$GeneID,featureRanges=exoninfo, transcripts=transcripts_l)
      4. Normalisatie en dispersieschatting. Voer vervolgens normalisatie uit tussen monsters en schat de variantie van de gegevens, vanwege zowel Poisson-tellingsruis van de discrete aard van RNA-seq als biologische variabiliteit, met behulp van de volgende commando's.
         ​dxd %<>% estimateSizeFactors %>% estimateDispersions %T>% plotDispEsts
      5. Test op differentieel gebruik. Test na de schatting van de variatie op differentieel exongebruik voor elk gen en genereer de resultaten.
         dxd%<>%testForDEU%>%estimateExonFoldChanges(fitExpToVar=
"condition")#Estimate fold changes
dxr=DEXSeqResults(dxd)
      6. Visualiseer splicinggebeurtenissen voor geselecteerde genen met behulp van de volgende opdracht.
         plotDEXSeq(dxr,"Wnk1", displayTranscripts=TRUE, splicing=TRUE,legend
=TRUE,cex.axis=1.2,cex=1.3,lwd=2)
         Onderzoek het R Notebook-bestand "AS_analysis_RNASeq.Rmd" om extra plots te genereren voor genen van belang en om vulkaanplots op verschillende drempels te genereren.
   2. Met behulp van diffSplice van Limma om differentiële splicing te identificeren. Volg het R Notebook bestand "AS_analysis_RNASeq.Rmd". Zorg ervoor dat de stappen 2.1-2.3 zijn gevolgd om invoerbestanden voor te bereiden voordat u verder gaat.
      1. Bibliotheken laden
         library(limma)
library(edgeR)
      2. Niet-specifieke filtering. Extraheer de matrix van afgelezen tellingen verkregen in 2.3. Maak een lijst met functies met behulp van de functie 'DGEList' van het edgeR-pakket, waarbij rijen genen vertegenwoordigen en kolommen monsters vertegenwoordigen.
         mycounts=countData$counts
#Change the rownames of the countdata to exon Ids instead of genes for unique rownames.
rownames(mycounts) = exoninfo$ExonID
dge<-DGEList(counts=mycounts)
#Filtering
isexpr<- rowSums(cpm(dge)>1) >=3
dge<-dge[isexpr,,keep.lib.sizes=FALSE]
#Extract the exon annotations for only transcripts meeting non-specific filter
exoninfo=anno%>% filter(ExonID%in% rownames(dge$counts))
#Convert the exoninfo into GRanges object
exoninfo1<-GRanges(seqnames=exoninfo$Chr,
ranges=IRanges(start=exoninfo$Start, end=exoninfo$End, width=exoninfo$Width),strand=Rle(exoninfo$Strand))
mcols(exoninfo1)$TranscriptIDs<-exoninfo$TranscriptIDs
mcols(exoninfo1)$Ticker<-exoninfo$Ticker
mcols(exoninfo1)$ExonID<-exoninfo$ExonID
mcols(exoninfo1)$n<-exoninfo$n
mcols(exoninfo1)$GeneID<-exoninfo$GeneID
transcripts_l= strsplit(exoninfo1$TranscriptIDs, "\\,")
         OPMERKING: Als niet-specifieke filterstap worden tellingen gefilterd door cpm < 1 in x uit n monsters, waarbij x het minimale aantal replicaties is in elke omstandigheid. n = 6 en x = 3 voor dit voorbeeld gegevens.
      3. Normaliseer de tellingen in monsters, met de functie 'calcNormFactors' uit het 'edgeR'-pakket met behulp van bijgesneden gemiddelde van M-waarden (TMM-normalisatiemethode)²⁴ Het berekent schaalfactoren om bibliotheekgrootten aan te passen.
         ​dge<-calcNormFactors(dge)
      4. Gebruik sampleTable zoals gegenereerd in stap 2.4.1.1 en maak de ontwerpmatrix. De ontwerpmatrix kenmerkt het ontwerp. Zie de Limma User Guide (https://www.bioconductor.org/packages/devel/bioc/vignettes/limma/inst/doc/usersguide.pdf) hoofdstukken 8 &9 voor meer informatie over ontwerpmatrices voor meer geavanceerde experimentele ontwerpen.
         Treat<- factor(sampleTable$condition)
design<- model.matrix(~0+Treat)
colnames(design) <- levels(Treat)
      5. Pas een lineair model per exon aan. Voer de 'voom'-functie van het 'limma'-pakket uit om RNA-seq-gegevens te verwerken om variantie te schatten en precisiegewichten te genereren om te corrigeren voor Poisson-tellingsruis, en transformeer de exon-niveautellingen naar log2-tellingen per miljoen (logCPM). Voer vervolgens lineaire modellering uit met behulp van de 'lmfit'-functie om lineaire modellen aan te passen aan de expressiegegevens voor elk exon. Bereken empirische Bayes-statistieken voor aangepast model met behulp van de 'eBayes'-functie om differentiële exonexpressie te detecteren. Definieer vervolgens een contrastmatrix voor de experimentele vergelijkingen van belang. Gebruik 'contrasten.fit' om coëfficiënten en standaardfouten voor elk vergelijkingspaar te verkrijgen.
         v<-voom(dge,design,plot=FALSE)
fit<-lmFit(v,design)
fit<-eBayes(fit)
colnames(fit)
cont.matrix<-makeContrasts(
Mbnl1_KO_WT=Mbnl1_KO-WT,
levels=design)
fit2<-contrasts.fit(fit,cont.matrix)
      6. Differentiële splicing analyse. Voer 'diffSplice' uit op het aangepaste model om de verschillen in exongebruik van genen tussen wild-type en knock-out te testen en de best gerangschikte resultaten te verkennen met behulp van de 'topSplice'-functie: test = "t" geeft een rangschikking van AS-exonen, test = "simes" geeft een rangschikking van genen.
         ex<-diffSplice(fit2,geneid=exoninfo$GeneID,exonid=exoninfo$ExonID)
ts<-topSplice(ex,n=Inf,FDR=0.1, test="t", sort.by="logFC")
​tg<-topSplice(ex,n=Inf,FDR=0.1, test="simes")
      7. Visualisatie. Plot de resultaten met de 'plotSplice'-functie, waardoor het gen van belang is voor het geneid-argument. Sla de topresultaten gesorteerd op log Vouw verandering in een object en genereer een vulkaanplot om de exonen tentoon te stellen.
         plotSplice(ex,geneid="Wnk1", FDR=0.1)
#Volcano plot
EnhancedVolcano(ts,lab=ts$ExonID,selectLab= head((ts$ExonID),2000), xlab= bquote(~Log[2]~'fold change'), x='logFC', y='P.Value', title='Volcano Plot', subtitle='Mbnl1_KO vs WT (Limma_diffSplice)', FCcutoff=2, labSize=4,legendPosition="right", caption= bquote(~Log[2]~"Fold change cutoff, 2; FDR 10%"))
   3. rMATS gebruiken
      1. Zorg ervoor dat de nieuwste versie van rMATS v4.1.1 (ook bekend als rMATS turbo vanwege de kortere verwerkingstijd en minder vereisten van het geheugen) is geïnstalleerd met behulp van conda of github (https://github.com/Xinglab/rmats-turbo/releases/download/v4.1.1/rmats_turbo_v4_1_1.tar.gz) in de werkmap. Volg paragraaf 4.3 in "AS_analysis_RNASeq.Rmd".
      2. Ga naar de map met bam-bestanden die zijn verkregen na toewijzing en bereid tekstbestanden voor, zoals vereist door rMATS, voor de twee voorwaarden door de naam van bam-bestanden (samen met het pad) gescheiden door ',' te kopiëren. De volgende opdrachten moeten worden uitgevoerd op de Linux-opdrachtregel:
         mkdir rMATS_analysis
cd bams/
ls -pd "$PWD"/*| grep "WT"| tr '\n'','> Wt.txt
ls -pd "$PWD"/*| grep "Mb"| tr '\n'','> KO.txt
mv *.txt ../rMATS_analysis
      3. Voer rmats.py uit met de twee invoerbestanden die in de vorige stap zijn gegenereerd, samen met het GTF-bestand dat is verkregen in 2.1.1.3. Dit genereert een uitvoermap 'rmats_out' met tekstbestanden die statistieken (p-waarden en inclusieniveaus) beschrijven voor elke splicinggebeurtenis afzonderlijk.
         python rmats-turbo/rmats.py --b1 KO.txt --b2 Wt.txt --gtf annotation.gtf -t paired --readLength 50 --nthread 8 --od rmats_out/ --tmp rmats_tmp --task pos
         OPMERKING: De referentieannotatie in de vorm van een GTF-bestand is ook vereist. Controleer de parameters als de gegevens single-end zijn en wijzig de optie -t dienovereenkomstig.
      4. RMATS-resultaten verkennen. Gebruik Bioconductor pakket '^maser'25 om de rMATS resultaten te verkennen. Laad de JCEC-tekstbestanden (Junction and Exon Counts) in het 'maser'-object en filter het resultaat op basis van dekking door ten minste vijf gemiddelde leeswaarden per splicinggebeurtenis op te nemen.
         library(maser)
mbnl1<-maser("/rmats_out/", c("WT","Mbnl1_KO"), ftype="JCEC")
#Filtering out events by coverage
mbnl1_filt<-filterByCoverage(mbnl1,avg_reads=5)
      5. Visualisatie van rMATS-resultaten. Selecteer de significante splicinggebeurtenissen bij False Discovery Rate (FDR) 10% en minimaal 10% verandering in Percent Spliced In (deltaPSI) met behulp van de functie 'topEvents' uit het 'maser'-pakket. Controleer vervolgens de gengebeurtenissen op individuele genen van belang (monstergen-Wnk1) en plot PSI-waarden voor elke splicinggebeurtenis van dat gen. Genereer een vulkaanplot door het gebeurtenistype op te geven.
         #Top splicing events at 10% FDR
mbnl1_top<-topEvents(mbnl1_filt,fdr=0.1, deltaPSI=0.1)
mbnl1_top
#Check the gene events for a particular gene
mbnl1_wnk1<-geneEvents(mbnl1_filt,geneS="Wnk1", fdr=0.1, deltaPSI=0.1)
maser::display(mbnl1_wnk1,"SE")
plotGenePSI(mbnl1_wnk1,type="SE", show_replicates
=TRUE)
​volcano(mbnl1_filt,fdr=0.1, deltaPSI=0.1,type="SE")
+xlab("deltaPSI")+ylab("Log10 Adj. Pvalue")+ggtitle("Volcano Plot of exon skipping events")
      6. Genereer Sashimi-plots voor het splicinggebeurtenissenresultaat verkregen met rMATS in de vorm van tekstbestanden met behulp van het 'rmats2shahimiplot'-pakket. Voer het python-script uit op de Linux-opdrachtregel.
         python ./src/rmats2sashimiplot/rmats2sashimiplot.py --b1 ../bams/WT_Thymus_1.bam,../bams/WT_Thymus_2.bam,../bams/WT_Thymus_3.bam --b2 ../bams/Mbnl1KO_Thymus_1.bam,../bams/Mbnl1KO_Thymus_2.bam,../bams/Mbnl1KO_Thymus_3.bam -t SE -e ../rMATS_analysis/rmats_out/SE.MATS.JC.txt --l1 WT --l2 Mbnl1_KO --exon_s 1 --intron_s 5 -o ../rMATS_analysis/rmats2shasmi_output
         OPMERKING: Dit proces kan tijdrovend zijn omdat het de Sashimi-plot genereert voor alle resultaten in het gebeurtenissenbestand. Kies de topresultaten (gennamen en exonen) zoals weergegeven door de topEvents-functie van 'maser' en visualiseer de bijbehorende Sashimi-plot.

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

1. Gegevens downloaden en voorbewerking
   OPMERKING: Raadpleeg het aanvullende R Notebook-bestand "APA_analysis_3PSeq_notebook. Rmd" voor de volledige opdrachten voor het downloaden en voorbewerkingsstappen van gegevens, of voer het aanvullende bash-bestand "APA_data_downloading_preprocessing.sh" uit op de Linux-opdrachtregel.
   1. Download gegevens van SRA met de Toetredings-Ids (1553129 tot 1553136).
   2. Trim adapters en omgekeerde aanvulling om de zinstrengvolgorde te verkrijgen.
      OPMERKING: Deze stap is specifiek voor de gebruikte PolyA-seq-test.
   3. Kaart leest naar muis genoom assemblage met behulp van bowtie aligner²⁶.
2. Het voorbereiden van pA-sitesannotaties.
   OPMERKING: De verwerking van het pA-siteannotatiebestand wordt eerst uitgevoerd met behulp van het aanvullende R Notebook-bestand "APA_analysis_3PSeq_notebook. Rmd" (2.1 - 2.6), en dan met behulp van bash bestand "APA_annotation_preparation.sh".
   1. Download de annotatie van pA-sites uit de PolyASite 2.0-database⁶.
   2. Selecteer pA-siteannotaties om 3'-untranslated region (UTR) pA-sites te behouden, die zijn geannoteerd als Terminal Exon (TE) of 1000 nt downstream van een geannoteerde terminal exon (DS) voor downstream-analyse.
   3. Verkrijg pA-sitepieken. Anker op elke pA-splitsingsplaats en visualiseer de gemiddelde leesdekking met bedgereedschappen en deeptools^27,28. De resultaten toonden aan dat de pieken van de in kaart gebrachte reads voornamelijk verspreid waren binnen ~60 bp stroomopwaarts van de splitsingsplaatsen (figuur 5 en aanvullende figuur 5). Daarom werden de coördinaten van pA-sites uitgebreid van het annotatiebestand tot 60 bp stroomopwaarts van hun splitsingsplaatsen. Afhankelijk van het specifieke 3'-eindsequencingprotocol dat wordt gebruikt, moet deze stap worden geoptimaliseerd voor andere assays dan PolyA-seq.
3. Tellen leest
   1. Bereid het annotatiebestand van pA-sites voor.
      anno<- read.table(file= "flanking60added.pA_annotation.bed",
stringsAsFactors=FALSE, check.names=FALSE, header=FALSE, sep="")
colnames(anno) <- c("chrom", "chromStart", "chromEnd", "name", "score", "strand", "rep", "annotation", "gene_name", "gene_id")
anno<- dplyr::select(anno,name,chrom, chromStart,chromEnd, strand,gene_id,gene_name,rep)
colnames(anno) <- c("GeneID", "Chr", "Start", "End", "Strand", "Ensembl", "Symbol", "repID")
   2. Pas 'featureCounts' toe om onbewerkte tellingen te verkrijgen. Sla de teltabel op als het RData-bestand "APA_countData.Rdata" voor APA-analyse met behulp van verschillende hulpprogramma's.
      countData<- dir("bamfiles", pattern="sorted.bam$", full.names=TRUE) %>%
# Read all bam files as input for featureCounts
featureCounts(annot.ext=anno, isGTFAnnotationFile= FALSE,minMQS=0,useMetaFeatures= TRUE,allowMultiOverlap=TRUE, largestOverlap= TRUE,strandSpecific=1, countMultiMappingReads =TRUE,primaryOnly= TRUE,isPairedEnd= FALSE,nthreads=12)%T>%
save(file="APA_countData.Rdata")
      OPMERKING: Wees bewust van het wijzigen van een van de parameters die worden vermeld in de functie 'featureCounts'. Wijzig de parameter 'strandSpecific' om ervoor te zorgen dat deze consistent is met de sequencingrichting van de gebruikte 3'-eindsequencingtest (empirisch gezien zal het visualiseren van de gegevens in een genoombrowser over genen op plus- en minstrengen dit verduidelijken).
   3. Pas niet-specifieke filtering van countData toe. Filtering kan de statistische robuustheid in differentiële pA-sitegebruikstests aanzienlijk verbeteren. Ten eerste hebben we die genen verwijderd met slechts één pA-site, waarop differentieel pA-sitegebruik niet kan worden gedefinieerd. Ten tweede passen we niet-specifieke filtering toe op basis van dekking: tellingen worden gefilterd door cpm minder dan 1 in x uit n monsters, waarbij x het minimale aantal replicaties is in elke omstandigheid. N = 8 en x = 2 voor dit voorbeeld gegevens.
      load(file= "APA_countData.Rdata")# Skip this step if already loaded
# Non-specific filtering: Remove the pA sites not differentially expressed in the samples
countData<-countData$counts%>%as.data.frame%>% .[rowSums(edgeR::cpm(.)>1) >=2, ]
anno%<>% .[.$GeneID%in% rownames(countData), ]
# Remove genes with only 1 site and NA in geneIDs
dnsites<-anno%>%group_by(Symbol)%>%summarise(nsites=n())%>% filter(nsites>1&!is.na(Symbol))
anno<-anno%>% filter(Symbol%in%dnsites$Symbol)
countData<-countData[rownames(countData) %in%anno$GeneID, ]
4. Differentiële polyadenylation site gebruiksanalyse met behulp van DEXSeq en diffSplice pipelines.
   1. DeXSeq-pakket gebruiken
      OPMERKING: Aangezien een contrastmatrix niet kan worden gedefinieerd voor de DEXSeq-pijplijn, moet de differentiële APA-analyse van elke twee experimentele omstandigheden afzonderlijk worden uitgevoerd. De differentiële APA-analyse van de aandoening WT en de voorwaarde DKO wordt als voorbeeld uitgevoerd om de procedure uit te leggen. Raadpleeg het aanvullend dossier "APA_analysis_3PSeq_notebook. Rmd" voor de stapsgewijze workflow van deze sectie en de differentiële APA-analyse van andere contrasten.
      1. Laad de bibliotheek en maak een voorbeeldtabel om het experimentele ontwerp te definiëren.
         c("DEXSeq", "GenomicRanges") %>% lapply(library, character.only=TRUE) %>%invisible
sampleTable1<- data.frame(row.names= c("WT_1","WT_2","DKO_1","DKO_2"),
condition= c(rep("WT", 2), rep("DKO", 2)),
​libType= rep("single-end", 4))
      2. Bereid het informatiebestand van pA-sites voor met behulp van Bioconductor-pakket GRanges.
         # Prepare the GRanges object for DEXSeqDataSet object construction
PASinfo <- GRanges(seqnames = anno$Chr,
ranges = IRanges(start = anno$Start, end = anno$End),strand = Rle(anno$Strand))
mcols(PASinfo)$PASID<-anno$repID
mcols(PASinfo)$GeneEns<-anno$Ensembl
mcols(PASinfo)$GeneID<-anno$Symbol
# Prepare the new feature IDs, replace the strand information with letters to match the current pA site clusterID
new.featureID <- anno$Strand %>% as.character %>% replace(. %in% "+", "F") %>% replace(. %in% "-", "R") %>% paste0(as.character(anno$repID), .)
      3. Gebruik de leestellingen die zijn gegenereerd in stap 3.3 en de pA-sitegegevens die zijn verkregen uit de vorige stap om het DEXSeq-object te maken.
         # Select the read counts of the condition WT and DKO
countData1<- dplyr::select(countData, SRR1553129.sorted.bam, SRR1553130.sorted.bam, SRR1553131.sorted.bam, SRR1553132.sorted.bam)
# Rename the columns of countData using sample names in sampleTable
colnames(countData1) <- rownames(sampleTable1)
dxd1<-DEXSeqDataSet(countData=countData1,
sampleData=sampleTable1,
design=~sample+exon+condition:exon,
featureID=new.featureID,
groupID=anno$Symbol,
featureRanges=PASinfo)
      4. Definieer het contrastpaar door de niveaus van voorwaarden in het DEXSeq-object te definiëren.
         dxd1$condition<- factor(dxd1$condition, levels= c("WT", "DKO"))
# The contrast pair will be "DKO - WT"
      5. Normalisatie en dispersieschatting. Net als bij RNA-seq-gegevens, voeren 3'-eindsequencinggegevens normalisatie uit tussen monsters (kolomgewijze mediaan van verhoudingen voor elk monster) met behulp van de functie 'estimateSizeFactors' en schatten de variatie van de gegevens met behulp van de functie 'estimateDispersions', en visualiseer vervolgens het dispersieschattingsresultaat met behulp van de functie 'plotDispEsts'.
         ​dxd1 %<>% estimateSizeFactors %>% estimateDispersions %T>% plotDispEsts
      6. Differentiële pA-sitegebruikstests voor elk gen met behulp van de functie 'testForDEU', en schatten vervolgens de vouwverandering van pA-sitegebruik met behulp van de functie 'estimateExonFoldChanges'. Controleer de resultaten met behulp van de functie 'DEXSeqResults' en stel 'FDR < 10%' in als criterium voor significant differentiële pA-sites.
         dxd1 %<>% testForDEU %>% estimateExonFoldChanges(fitExpToVar = "condition")
dxr1 <- DEXSeqResults(dxd1)
dxr1
mcols(dxr1)$description
table(dxr1$padj<0.1) # Check the number of differential pA sites (FDR < 0.1)
table(tapply(dxr1$padj<0.1, dxr1$groupID, any)) # Check the number of gene overlapped with differential pA site
      7. Visualisatie van de differentiële pA-sitegebruiksresultaten met behulp van differentiële APA-plots gegenereerd door de functie 'plotDEXSeq' en vulkaanplot door de functie 'EnhancedVolcano'.
         # Select the top 100 significant differential pA sites ranked by FDR
topdiff.PAS<- dxr1%>%as.data.frame%>%rownames_to_column%>%arrange(padj)%$%groupID[1:100]

# Apply plotDEXSeq for the visualization of differential polyA usage
plotDEXSeq(dxr1,"S100a7a", legend=TRUE, expression=FALSE,splicing=TRUE, cex.axis=1.2, cex=1.3,lwd=2)

# Apply perGeneQValue to check the top genes with differential polyA site usage
dxr1%<>% .[!is.na(.$padj), ]
dgene<- data.frame(perGeneQValue= perGeneQValue (dxr1)) %>%rownames_to_column("groupID")

dePAS_sig1<-dxr1%>% data.frame() %>%
dplyr::select(-matches("dispersion|stat|countData|genomicData"))%>%
inner_join(dgene)%>%arrange(perGeneQValue)%>%distinct()%>%
filter(padj<0.1)

# Apply EnhancedVolcano package to visualise differential polyA site usage
"EnhancedVolcano"%>% lapply(library, character.only=TRUE) %>%invisible
EnhancedVolcano(dePAS_sig1, lab=dePAS_sig1$groupID, x='log2fold_DKO_WT',
y='pvalue',title='Volcano Plot',subtitle='DKO vs WT',
FCcutoff=1,labSize=4, legendPosition="right",
caption= bquote(~Log[2]~"Fold change cutoff, 1; FDR 10%"))
   2. Met behulp van diffSplice pakket. Raadpleeg het aanvullende R Notebook-bestand "APA_analysis_3PSeq_notebook. Rmd" voor de stapsgewijze workflow van deze sectie.
      1. Definieer interessante contrasten voor differentiële pA-gebruiksanalyse.
         OPMERKING: Deze stap moet worden uitgevoerd na de constructie en verwerking van het DGEList-object, dat is opgenomen in het R Notebook-bestand "APA_analysis_3PSeq_notebook. Rmd".
         contrast.matrix<-makeContrasts(DKO_vs_WT=DKO-
WT,Ctrl_vs_DKO=Ctrl-DKO,
KD_vs_Ctrl=KD-Ctrl,KD_vs_DKO=KD-DKO,levels=design)
fit2<-fit%>%contrasts.fit(contrast.matrix)%>%eBayes
summary(decideTests(fit2))
ex<-diffSplice(fit2,geneid=anno$Symbol,exonid=new.featureID)
topSplice(ex) #Check the top significant results with topSplice
      2. Visualiseer het resultaat van interessante contrasten (hier "DKO - WT") met behulp van differentiële APA-plots door de functie 'plotSplice' en vulkaanplots met de functie 'EnhancedVolcano'. Raadpleeg het R Notebook-bestand "APA_analysis_3PSeq_notebook. Rmd" 4.2.7 - 4.2.9 voor de visualisatie van andere contrastparen.
         sig1<-topSplice(ex,n=Inf,FDR=0.1,coef=1, test="t", sort.by="logFC")
sig1.genes<-topSplice(ex,n=Inf,FDR=0.1,coef=1, test="simes")
plotSplice(ex, coef=1,geneid="S100a7a", FDR = 0.1)
plotSplice(ex,coef=1,geneid="Tpm1", FDR = 0.1)
plotSplice(ex,coef=1,geneid="Smc6", FDR = 0.1)
EnhancedVolcano(sig1, lab=sig1$GeneID,xlab= bquote(~Log[2]~'fold change'),
x='logFC', y='P.Value', title='Volcano Plot', subtitle='DKO vs WT',
FCcutoff=1, labSize=6, legendPosition="right")

Resultaten

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

ALS:
De belangrijkste output van de AS-analyse (aanvullende tabel 1 voor diffSplice; Tabel 2 voor DEXSeq) is een lijst van exonen die differentieel gebruik tussen omstandigheden laten zien, en een lijst van genen die significante...

Discussie

In deze studie evalueerden we exon-gebaseerde en event-gebaseerde benaderingen om AS en APA te detecteren in bulk RNA-Seq en 3'-end sequencinggegevens. De exon-gebaseerde AS-benaderingen produceren zowel een lijst van differentieel tot expressie gebrachte exonen als een rangschikking op genniveau, geordend op basis van de statistische significantie van de totale differentiële splicingactiviteit op genniveau (tabellen 1-2, 4-5). Voor het diffSplice-pakket wordt het differentiële gebruik bepaald door gewog...

Materialen

Lijst van materialen gebruikt in dit artikel

Naam	Bedrijf	Catalogusnummer	Opmerkingen
Niet relevant voor computationeel onderzoek