Identification d’un épissage alternatif et d’une polyadénylation dans les données de séquençage de l’ARN

Résumé

L’épissage alternatif (AS) et la polyadénylation alternative (APA) élargissent la diversité des isoformes de transcrits et de leurs produits. Ici, nous décrivons des protocoles bioinformatiques pour analyser les tests de séquençage d’ARN en vrac et de séquençage à 3' pour détecter et visualiser l’AS et l’APA variant selon les conditions expérimentales.

Résumé

En plus de l’analyse typique du RNA-Seq pour mesurer l’expression génique différentielle (DGE) dans des conditions expérimentales / biologiques, les données RNA-seq peuvent également être utilisées pour explorer d’autres mécanismes de régulation complexes au niveau de l’exon. L’épissage alternatif et la polyadénylation jouent un rôle crucial dans la diversité fonctionnelle d’un gène en générant différentes isoformes pour réguler l’expression génique au niveau post-transcriptionnel, et limiter les analyses à l’ensemble du niveau du gène peut manquer cette couche régulatrice importante. Ici, nous démontrons des analyses détaillées étape par étape pour l’identification et la visualisation de l’utilisation différentielle des exons et des sites de polyadénylation dans toutes les conditions, en utilisant Bioconductor et d’autres packages et fonctions, y compris DEXSeq, diffSplice du package Limma et rMATS.

Introduction

Le séquençage de l’ARN a été largement utilisé au fil des ans, généralement pour estimer l’expression différentielle des gènes et la découverte de gènes¹. En outre, il peut également être utilisé pour estimer l’utilisation variable du niveau d’exon en raison de l’expression de gènes différentes isoformes, contribuant ainsi à une meilleure compréhension de la régulation des gènes au niveau post-transcriptionnel. La majorité des gènes eucaryotes génèrent différentes isoformes par épissage alternatif (AS) pour augmenter la diversité de l’expression de l’ARNm. Les événements AS peuvent être divisés en différents modèles: saut d’exons complets (SE) où un exon (« cassette ») est complètement retiré de la transcription avec ses introns flanquants; sélection alternative (donneur) du site d’épissage 5' (A5SS) et alternative 3' (accepteur) sélection du site d’épissage (A3SS) lorsque deux sites d’épissage ou plus sont présents à chaque extrémité d’un exon; rétention d’introns (RI) lorsqu’un intron est retenu dans le transcrit d’ARNm mature et exclusion mutuelle de l’utilisation d’exons (MXE) où un seul des deux exons disponibles peut être retenu à la fois ^2,3. La polyadénylation alternative (APA) joue également un rôle important dans la régulation de l’expression génique en utilisant des sites poly (A) alternatifs pour générer plusieurs isoformes d’ARNm à partir d’un seul transcrit⁴. La plupart des sites de polyadénylation (pA) sont situés dans la région 3' non traduite (3' UTR), générant des isoformes d’ARNm avec diverses longueurs 3' UTR. Comme l’UTR 3' est le centre central pour la reconnaissance des éléments régulateurs, différentes longueurs 3' UTR peuvent affecter la localisation, la stabilité et la traduction de l’ARNm⁵. Il existe une classe de tests de séquençage d’extrémité 3' optimisés pour détecter l’APA qui diffèrent dans les détails du protocole⁶. Le pipeline décrit ici est conçu pour PolyA-seq, mais peut être adapté à d’autres protocoles comme décrit.

Dans cette étude, nous présentons un pipeline de méthodes d’analyse différentielle^des exons ^7,8 (Figure 1), qui peuvent être divisées en deux grandes catégories : basées sur les exons (DEXSeq^9, diffSplice¹⁰) et basées sur les événements (analyse multivariée répliquée de l’épissage des transcriptions (rMATS)¹¹). Les méthodes basées sur les exons comparent le changement de pli entre les conditions des exons individuels, à une mesure du changement global du pli des gènes pour appeler l’utilisation d’exons exprimée différentiellement, et à partir de là, calculez une mesure au niveau du gène de l’activité SA. Les méthodes basées sur les événements utilisent des lectures de jonction couvrant exon-intron pour détecter et classer des événements d’épissage spécifiques tels que le saut d’exon ou la rétention d’introns, et distinguer ces types AS dans la sortie³. Ainsi, ces méthodes fournissent des points de vue complémentaires pour une analyse complète de la SA^12,13. Nous avons sélectionné DEXSeq (basé sur le package DESeq2¹⁴ DGE) et diffSplice (basé sur le package Limma¹⁰ DGE) pour l’étude car ils sont parmi les packages les plus largement utilisés pour l’analyse d’épissage différentiel. rMATS a été choisi comme méthode populaire pour l’analyse basée sur les événements. Une autre méthode populaire basée sur les événements est MISO (Mixture of Isoforms)¹. Pour l’APA, nous adaptons l’approche basée sur les exons.

Graphique 1. Pipeline d’analyse. Organigramme des étapes utilisées dans l’analyse. Les étapes comprennent : l’obtention des données, l’exécution de contrôles de qualité et d’alignement des lectures, suivis du comptage des lectures à l’aide d’annotations pour les exons, les introns et les sites pA connus, le filtrage pour supprimer les faibles nombres et la normalisation. Les données PolyA-seq ont été analysées pour d’autres sites pA à l’aide des méthodes diffSplice/DEXSeq, le RNA-Seq en vrac a été analysé pour l’épissage alternatif au niveau de l’exon avec les méthodes diffSplice/DEXseq et les événements AS analysés avec rMATS. Veuillez cliquer ici pour voir une version agrandie de cette figure.

Les données RNA-seq utilisées dans cette enquête ont été acquises à partir de Gene Expression Omnibus (GEO) (GSE138691)¹⁵. Nous avons utilisé les données de séquençage de l’ARN de souris de cette étude avec deux groupes de conditions : le type sauvage (WT) et le knockout de type 1 de type Muscleblind (Mbnl1 KO) avec trois réplications chacun. Pour démontrer l’analyse différentielle de l’utilisation du site de polyadénylation, nous avons obtenu des données PolyA-seq sur des fibroblastes embryonnaires de souris (MEF) (GEO Accession GSE60487)¹⁶. Les données comportent quatre groupes de conditions : Wild-type (WT), Muscleblind-like type1/type 2 double knockout (Mbnl1/2 DKO), Mbnl 1/2 DKO avec Mbnl3 knockdown (KD) et Mbnl1/2 DKO avec contrôle Mbnl3 (Ctrl). Chaque groupe de conditions se compose de deux répétitions.

	Adhésion GEO	Numéro d’exécution SRA	Nom de l’échantillon	Condition	Répliquer	Tissu	Séquençage	Longueur de lecture
RNA-Seq	GSM4116218	SRR10261601	Mbnl1KO_Thymus_1	Mbnl1 knockout	Rép. 1	Thymus	Extrémité jumelée	100 pb
	GSM4116219	SRR10261602	Mbnl1KO_Thymus_2	Mbnl1 knockout	Rép. 2	Thymus	Extrémité jumelée	100 pb
	GSM4116220	SRR10261603	Mbnl1KO_Thymus_3	Mbnl1 knockout	Rép. 3	Thymus	Extrémité jumelée	100 pb
	GSM4116221	SRR10261604	WT_Thymus_1	Type sauvage	Rép. 1	Thymus	Extrémité jumelée	100 pb
	GSM4116222	SRR10261605	WT_Thymus_2	Type sauvage	Rép. 2	Thymus	Extrémité jumelée	100 pb
	GSM4116223	SRR10261606	WT_Thymus_3	Type sauvage	Rép. 3	Thymus	Extrémité jumelée	100 pb
3P-Seq	GSM1480973	SRR1553129	WT_1	Type sauvage (WT)	Rép. 1	Fibroblastes embryonnaires de souris (MEF)	Extrémité unique	40 pb
	GSM1480974	SRR1553130	WT_2	Type sauvage (WT)	Rép. 2	Fibroblastes embryonnaires de souris (MEF)	Extrémité unique	40 pb
	GSM1480975	SRR1553131	DKO_1	Mbnl 1/2 double knockout (DKO)	Rép. 1	Fibroblastes embryonnaires de souris (MEF)	Extrémité unique	40 pb
	GSM1480976	SRR1553132	DKO_2	Mbnl 1/2 double knockout (DKO)	Rép. 2	Fibroblastes embryonnaires de souris (MEF)	Extrémité unique	40 pb
	GSM1480977	SRR1553133	DKOsiRNA_1	Mbnl 1/2 double knockout avec Mbnl 3 siRNA (KD)	Rép. 1	Fibroblastes embryonnaires de souris (MEF)	Extrémité unique	40 pb
	GSM1480978	SRR1553134	DKOsiRNA_2	Mbnl 1/2 double knockout avec Mbnl 3 siRNA (KD)	Rép. 2	Fibroblastes embryonnaires de souris (MEF)	Extrémité unique	36 pb
	GSM1480979	SRR1553135	DKONTsiRNA_1	Mbnl 1/2 double knockout avec siRNA non ciblé (Ctrl)	Rép. 1	Fibroblastes embryonnaires de souris (MEF)	Extrémité unique	40 pb
	GSM1480980	SRR1553136	DKONTsiRNA_2	Mbnl 1/2 double knockout avec siRNA non ciblé (Ctrl)	Rép. 2	Fibroblastes embryonnaires de souris (MEF)	Extrémité unique	40 pb

Tableau 1. Résumé des ensembles de données RNA-Seq et PolyA-seq utilisés pour l’analyse.

Protocole

1. Installation des outils et des packages R utilisés dans l’analyse

1. Conda est un gestionnaire de paquets populaire et flexible qui permet une installation pratique des paquets avec leurs dépendances sur toutes les plates-formes. Utilisez 'Anaconda' (gestionnaire de paquets conda) pour installer 'conda' qui peut être utilisé pour installer les outils/paquets requis pour l’analyse.
2. Téléchargez 'Anaconda' selon la configuration système requise de https://www.anaconda.com/products/individual#Downloads et installez-le en suivant les instructions de l’installateur graphique. Installez tous les paquets requis en utilisant 'conda' en tapant ce qui suit sur la ligne de commande Linux.
   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. Pour télécharger tous les paquets R utilisés dans le protocole, tapez le code suivant dans la console R (démarrée sur la ligne de commande Linux en tapant 'R') ou la console Rstudio.
   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)
}
   REMARQUE: Dans ce protocole de calcul, les commandes seront données sous forme de fichiers R Notebook (fichiers avec extension « . Rmd »), Fichiers de code R (fichiers avec extension « . R »), ou des scripts shell Linux Bash (fichiers avec l’extension « .sh »). Les fichiers R Notebook (Rmd) doivent être ouverts dans RStudio à l’aide de Fichier| Ouvrez le fichier..., et les morceaux de code individuels (qui peuvent être des commandes R ou des commandes shell Bash) s’exécutent ensuite de manière interactive en cliquant sur la flèche verte en haut à droite. Les fichiers de code R peuvent être exécutés en ouvrant dans RStudio, ou sur la ligne de commande Linux en préfaçant avec « Rscript », par exemple Rscript example.R. Les scripts Shell sont exécutés sur la ligne de commande Linux en préfaceant le script avec la commande « sh », par exemple.sh example.sh.

2. Analyse d’épissage alternatif (AS) à l’aide du séquençage de l’ARN

1. Téléchargement et prétraitement des données
   Remarque : Les extraits de code annotés ci-dessous sont disponibles dans le fichier de code supplémentaire »AS_analysis_RNASeq.Rmd« , pour suivre les étapes individuelles de manière interactive, et sont également fournis sous forme de script bash à exécuter par lots sur la ligne de commande Linux (sh downloading_data_preprocessing.sh).
   1. Téléchargement des données brutes.
      1. Téléchargez les données brutes à partir de Sequence Read Archive (SRA) à l’aide de la commande 'prefetch' de la boîte à outils SRA (v2.10.8)¹⁷. Donnez les ID d’accession SRA dans l’ordre dans la commande suivante pour les télécharger en parallèle en utilisant l’utilitaire parallèle^{GNU 18}. Pour télécharger des fichiers SRA d’ID d’accession de SRR10261601 à SRR10261606 en parallèle, utilisez la commande suivante sur la ligne de commande Linux.
         ​seq 10261601 10261606 | parallel prefetch SRR{}
      2. Extrayez les fichiers fastq de l’archive à l’aide de la fonction 'fastq-dump' de la boîte à outils SRA. Utilisez GNU en parallèle et donnez les noms de tous les fichiers SRA ensemble.
         ​parallel -j 3 fastq-dump --gzip --skip-technical --read-filter pass --dumpbase --split-e --clip --origfmt {} :::
      3. Téléchargez le génome de référence et les annotations pour Mouse (assemblage du génome GRCm39) à partir de www.ensembl.org en utilisant ce qui suit sur la ligne de commande Linux.
         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. Prétraitement et mappage des lectures à l’assemblage du génome
      1. Contrôle qualité. Évaluez la qualité des lectures brutes avec FASTQC (FASTQ Quality Check v0.11.9)¹⁹. Créez un dossier de sortie et exécutez fastqc en parallèle sur plusieurs fichiers fasta d’entrée. Cette étape générera un rapport de qualité pour chaque échantillon. Examinez les rapports pour vous assurer que la qualité des lectures est acceptable avant de procéder à une analyse plus approfondie. (Reportez-vous au manuel de l’utilisateur pour comprendre les rapports à https://www.bioinformatics.babraham.ac.uk/projects/fastqc/)
         mkdir fastqc_out
parallel "fastqc {} -o fastqc_out" ::: $RAW_DATA/*.fastq.gz
         REMARQUE: Si nécessaire, effectuez le découpage de l’adaptateur avec 'cutadapt'20 ou 'trimmomatic'21 pour supprimer le séquençage dans les adaptateurs latéraux, qui varie en fonction de la taille du fragment d’ARN et de la longueur de lecture. Dans cette analyse, nous avons sauté cette étape car la fraction de lectures affectées était minime.
      2. Alignement de lecture. L’étape suivante du prétraitement comprend la mise en correspondance des lectures avec le génome de référence. Tout d’abord, construisez l’index pour le génome de référence en utilisant la fonction 'genomeGenerate' de STAR²² , puis alignez les lectures brutes sur la référence (alternativement, des index prédéfinis sont disponibles sur le site Web de STAR et peuvent être utilisés directement pour l’alignement). Exécutez les commandes suivantes sur la ligne de commande Linux.
         #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
         Remarque : L’aligneur STAR génère et trie les fichiers BAM (Binary Alignment Map) pour chaque échantillon après l’alignement en lecture. Les fichiers Bam doivent être triés avant de passer aux étapes suivantes.
2. Préparation des annotations Exon.
   1. Exécutez le fichier de code supplémentaire « prepare_mm_exon_annotation. R » avec l’annotation téléchargée au format GTF (Gene transfer format) pour préparer les annotations. Pour exécuter, tapez ce qui suit sur la ligne de commande Linux.
      Rscript prepare_mm_exon_annotation.R annotation.gtf
      Remarque : Le fichier GTF contient plusieurs entrées d’exon pour différentes isoformes. Ce fichier est utilisé pour « réduire » les multiples ID de transcription pour chaque exon. C’est une étape importante pour définir les bacs de comptage d’exons.
3. Comptage lit. L’étape suivante consiste à compter le nombre de lectures mappées à différentes transcriptions / exons. Voir Fichier supplémentaire: « AS_analysis_RNASeq.Rmd ».
   1. Chargez les bibliothèques requises :
      packages<- c("Rsubread","tidyverse", "magrittr", "EnhancedVolcano", "edgeR","openxlsx")
invisible(lapply(packages, library, character.only=TRUE))
   2. Chargez le fichier d’annotation traité obtenu à partir de l’étape précédente (2.2).
      load("mm_exon_anno.RData")
   3. Lisez tous les fichiers bam obtenus à l’étape 2.2.2 en entrée pour 'featureCounts' pour compter les lectures. Lisez le dossier contenant les fichiers bam en répertoriant d’abord chaque fichier du répertoire qui se termine par .bam. Utilisez 'featureCounts' du paquet Rsubread qui prend les fichiers bam et l’annotation GTF traitée (référence) comme entrée pour générer une matrice de comptages associés à chaque entité avec des lignes représentant des exons (entités) et des colonnes représentant des échantillons.
      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. Ensuite, effectuez un filtrage non spécifique pour éliminer les exons faiblement exprimés (« non spécifique » indique que l’information sur l’état expérimental n’est pas utilisée dans le filtrage, afin d’éviter les biais de sélection). Transformez les données de l’échelle brute en nombres par million (cpm) à l’aide de la fonction cpm du package 'edgeR'²³ et conservez les exons avec des nombres supérieurs à un seuil définissable (pour cet ensemble de données, un cpm est utilisé) dans au moins trois échantillons. Supprimez également les gènes avec un seul 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, ]
      Remarque : Vérifiez les paramètres requis pour featureCounts lors de l’utilisation de différentes données, par exemple, pour les lectures à une extrémité, définissez 'isPairedEnd = FALSE'. Reportez-vous au guide de l’utilisateur RSubread pour choisir les options pour vos données, et consultez la section Discussion ci-dessous.
4. Épissage différentiel et analyse de l’utilisation des exons. Nous décrivons deux alternatives pour cette étape : DEXSeq et DiffSplice. L’un ou l’autre peut être utilisé et donner des résultats similaires. Pour des raisons de cohérence, sélectionnez DEXSeq si vous préférez un package DESeq2 pour DGE et utilisez DiffSplice pour une analyse DGE basée sur Limma. Voir dossier complémentaire : »AS_analysis_RNASeq.Rmd".
   1. Utilisation du package DEXSeq pour l’analyse différentielle des exons.
      1. Chargez la bibliothèque et créez un exemple de table pour définir le plan expérimental.
         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")))
         Remarque : les noms de ligne doivent être cohérents avec les noms de fichier bam utilisés par featureCounts pour compter les lectures. sampleTable se compose des détails de chaque échantillon, y compris : le type de bibliothèque et l’état. Ceci est nécessaire pour définir les contrastes ou le groupe de test pour détecter l’utilisation différentielle.
      2. Préparez le fichier d’informations sur l’exon. Les informations d’exon sous forme d’objets GRanges (plages génomiques) (https://bioconductor.org/packages/release/bioc/html/GenomicRanges.html) sont requises en tant qu’entrée pour créer l’objet DEXSeq à l’étape suivante. Faites correspondre les ID de gène avec le nombre de lectures pour créer l’objet exoninfo.
         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. Créez un objet DEXSeq à l’aide de la fonction DEXSeqDataSet. L’objet DEXSeq collecte ensemble le nombre de lectures, les informations sur les entités d’exon et les exemples d’informations. Utilisez les comptes de lecture générés à l’étape 3 et les informations d’exon obtenues à l’étape précédente pour créer l’objet DEXSeq à partir de la matrice de comptage. L’argument sampleData prend une entrée de bloc de données définissant les exemples (et leurs attributs : type et condition de bibliothèque), 'design' utilise sampleData pour générer une matrice de conception pour les tests différentiels à l’aide de la notation de formule de modèle. Notez qu’un terme d’interaction significatif, condition:exon, indique que la fraction de lectures sur un gène tombant sur un exon particulier dépend de la condition expérimentale, c’est-à-dire qu’il y a AS. Consultez la documentation DEXSeq pour une description complète de la définition de la formule du modèle pour les plans expérimentaux plus complexes. Pour les informations sur les caractéristiques, les identifiants d’exon, le gène correspondant et les transcriptions sont requis.
         ​dxd<-DEXSeqDataSet(countData$counts,sampleData=sampleTable, design=~sample+exon+condition:exon,featureID=exoninfo$ExonID,groupID=exoninfo$GeneID,featureRanges=exoninfo, transcripts=transcripts_l)
      4. Normalisation et estimation de la dispersion. Ensuite, effectuez une normalisation entre les échantillons et estimez la variance des données, due à la fois au bruit de comptage de Poisson dû à la nature discrète du séquençage de l’ARN et à la variabilité biologique, à l’aide des commandes suivantes.
         ​dxd %<>% estimateSizeFactors %>% estimateDispersions %T>% plotDispEsts
      5. Testez l’utilisation différentielle. Après l’estimation de la variation, tester l’utilisation différentielle de l’exon pour chaque gène et générer les résultats.
         dxd%<>%testForDEU%>%estimateExonFoldChanges(fitExpToVar=
"condition")#Estimate fold changes
dxr=DEXSeqResults(dxd)
      6. Visualisez les événements d’épissage pour les gènes sélectionnés à l’aide de la commande suivante.
         plotDEXSeq(dxr,"Wnk1", displayTranscripts=TRUE, splicing=TRUE,legend
=TRUE,cex.axis=1.2,cex=1.3,lwd=2)
         Examinez le fichier R Notebook « AS_analysis_RNASeq.Rmd » pour générer des tracés supplémentaires pour les gènes d’intérêt et pour générer des tracés de volcans à différents seuils.
   2. Utilisation de diffSplice de Limma pour identifier l’épissage différentiel. Suivez le fichier R Notebook »AS_analysis_RNASeq.Rmd". Assurez-vous que les étapes 2.1 à 2.3 ont été suivies pour préparer les fichiers d’entrée avant de poursuivre.
      1. Charger des bibliothèques
         library(limma)
library(edgeR)
      2. Filtrage non spécifique. Extraire la matrice des comptages de lecture obtenue en 2.3. Créez une liste d’entités à l’aide de la fonction 'DGEList' à partir du package edgeR, où les lignes représentent les gènes et les colonnes représentent les échantillons.
         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, "\\,")
         REMARQUE : En tant qu’étape de filtrage non spécifique, les comptages sont filtrés par cpm < 1 dans x sur n échantillons, où x est le nombre minimum de répétitions dans n’importe quelle condition. n = 6 et x = 3 pour cet exemple de données.
      3. Normalisez les comptages entre les échantillons, avec la fonction 'calcNormFactors' du package 'edgeR' en utilisant la moyenne tronquée des valeurs M (méthode de normalisation TMM)²⁴ Il calculera les facteurs d’échelle pour ajuster la taille des bibliothèques.
         ​dge<-calcNormFactors(dge)
      4. Utilisez sampleTable tel que généré à l’étape 2.4.1.1 et créez la matrice de conception. La matrice de conception caractérise la conception. Voir les chapitres 8 et 9 du Guide de l’utilisateur de Limma (https://www.bioconductor.org/packages/devel/bioc/vignettes/limma/inst/doc/usersguide.pdf) pour plus de détails sur les matrices de conception pour des conceptions expérimentales plus avancées.
         Treat<- factor(sampleTable$condition)
design<- model.matrix(~0+Treat)
colnames(design) <- levels(Treat)
      5. Ajuster un modèle linéaire par exon. Exécutez la fonction 'voom' du package 'limma' pour traiter les données RNA-seq afin d’estimer la variance et de générer des poids de précision pour corriger le bruit de comptage de Poisson, et transformez les comptages au niveau des exons en log2-counts par million (logCPM). Exécutez ensuite la modélisation linéaire à l’aide de la fonction 'lmfit' pour ajuster les modèles linéaires aux données d’expression pour chaque exon. Calculez des statistiques bayésiennes empiriques pour un modèle ajusté à l’aide de la fonction 'eBayes' pour détecter l’expression différentielle des exons. Ensuite, définissez une matrice de contraste pour les comparaisons expérimentales d’intérêt. Utilisez 'contrasts.fit' pour obtenir des coefficients et des erreurs-types pour chaque paire de comparaison.
         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. Analyse d’épissage différentiel. Exécutez 'diffSplice' sur le modèle ajusté pour tester les différences dans l’utilisation des gènes d’exons entre le type sauvage et le knockout et explorez les résultats les mieux classés en utilisant la fonction 'topSplice': test="t » donne un classement des exons AS, test="simes » donne un classement des gènes.
         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. Visualisation. Tracez les résultats avec la fonction 'plotSplice', en donnant le gène d’intérêt dans l’argument généide. Enregistrez les meilleurs résultats triés par log Pliez le changement dans un objet et générez un tracé de volcan pour exposer les exons.
         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. Utilisation de rMATS
      1. Assurez-vous que la dernière version de rMATS v4.1.1 (également appelée rMATS turbo en raison du temps de traitement réduit et des besoins en mémoire) est installée à l’aide de conda ou github (https://github.com/Xinglab/rmats-turbo/releases/download/v4.1.1/rmats_turbo_v4_1_1.tar.gz) dans le répertoire de travail. Suivez la section 4.3 dans « AS_analysis_RNASeq.Rmd ».
      2. Accédez au dossier contenant les fichiers bam obtenus après le mappage et préparez les fichiers texte, comme requis par rMATS, pour les deux conditions en copiant le nom des fichiers bam (ainsi que le chemin) séparés par ','. Les commandes suivantes doivent être exécutées sur la ligne de commande Linux :
         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. Exécutez rmats.py avec les deux fichiers d’entrée générés à l’étape précédente, ainsi que le fichier GTF obtenu dans 2.1.1.3. Cela générera un dossier de sortie 'rmats_out' contenant des fichiers texte décrivant les statistiques (valeurs p et niveaux d’inclusion) pour chaque événement d’épissage séparément.
         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
         REMARQUE : L’annotation de référence sous la forme d’un fichier GTF est également requise. Vérifiez les paramètres si les données sont à une extrémité unique et modifiez l’option -t en conséquence.
      4. Exploration des résultats rMATS. Utilisez le package Bioconductor '^maser'25 pour explorer les résultats rMATS. Chargez les fichiers texte JCEC (Junction and exon counts) dans l’objet « maser » et filtrez le résultat en fonction de la couverture en incluant au moins cinq lectures moyennes par événement d’épissage.
         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. Visualisation des résultats rMATS. Sélectionnez les événements d’épissage significatifs à un taux de faux découvert (FDR) de 10 % et une variation minimale de 10 % du pourcentage d’épissage entrant (deltaPSI) à l’aide de la fonction « topEvents » du package « maser ». Ensuite, vérifiez les événements génétiques pour les gènes d’intérêt individuels (échantillon de gène - Wnk1) et tracez les valeurs PSI pour chaque événement d’épissage de ce gène. Générez un tracé de volcan en spécifiant le type d’événement.
         #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. Générez des tracés Sashimi pour le résultat des événements d’épissage obtenu avec rMATS sous forme de fichiers texte en utilisant le package 'rmats2shahimiplot'. Exécutez le script python sur la ligne de commande Linux.
         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
         REMARQUE: Ce processus peut prendre beaucoup de temps car il générera le tracé Sashimi pour tous les résultats du fichier d’événements. Choisissez les meilleurs résultats (noms de gènes et exons) tels qu’affichés par la fonction topEvents de 'maser' et visualisez le diagramme de Sashimi correspondant.

3. Analyse alternative de polyadénylation (APA) utilisant le séquençage de l’extrémité 3'

1. Téléchargement et prétraitement des données
   REMARQUE : reportez-vous au fichier R Notebook supplémentaire « APA_analysis_3PSeq_notebook. Rmd » pour les commandes complètes pour les étapes de téléchargement et de prétraitement des données, ou exécutez le fichier bash supplémentaire « APA_data_downloading_preprocessing.sh » sur la ligne de commande Linux.
   1. Téléchargez les données de SRA avec les ID d’accession (1553129 à 1553136).
   2. Ajuster les adaptateurs et inverser le complément pour obtenir la séquence de brin de détection.
      REMARQUE: Cette étape est spécifique au test PolyA-seq utilisé.
   3. La carte se lit à l’assemblage du génome de la souris à l’aide de l’aligneur nœud papillon²⁶.
2. Préparation des annotations des sites pA.
   REMARQUE: Le traitement du fichier d’annotation de site pA est effectué d’abord à l’aide du fichier R Notebook supplémentaire « APA_analysis_3PSeq_notebook. Rmd » (2.1 - 2.6), puis en utilisant le fichier bash « APA_annotation_preparation.sh ».
   1. Téléchargez l’annotation des sites pA à partir de la base de données PolyASite 2.0⁶.
   2. Sélectionnez les annotations de site pA pour conserver les sites pA de région 3' non traduite (UTR), qui sont annotés comme Terminal Exon (TE) ou 1000 nt en aval d’un exon terminal annoté (DS) pour l’analyse en aval.
   3. Obtenez les pics de site pA. Ancrer à chaque site de clivage pA et visualiser la couverture moyenne de lecture à l’aide de bedtools et deeptools^27,28. Les résultats ont montré que les pics des lectures cartographiées étaient principalement dispersés à ~60 pb en amont des sites de clivage (figure 5 et figure supplémentaire 5). Par conséquent, les coordonnées des sites pA ont été étendues du fichier d’annotation à 60 pb en amont de leurs sites de clivage. En fonction du protocole de séquençage d’extrémité 3' spécifique utilisé, cette étape devra être optimisée pour les tests autres que PolyA-seq.
3. Comptage des lectures
   1. Préparez le fichier d’annotation des sites pA.
      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. Appliquez 'featureCounts' pour acquérir des nombres bruts. Enregistrez la table de comptage en tant que fichier RData « APA_countData.Rdata » pour l’analyse APA à l’aide de différents outils.
      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")
      REMARQUE: Soyez conscient de modifier l’un des paramètres répertoriés dans la fonction 'featureCounts'. Modifier le paramètre 'strandSpecific' pour s’assurer qu’il est cohérent avec la direction de séquençage du test de séquençage de l’extrémité 3' utilisé (empiriquement, la visualisation des données dans un navigateur de génome sur les gènes sur les brins plus et moins clarifiera cela).
   3. Appliquez un filtrage non spécifique de countData. Le filtrage peut améliorer considérablement la robustesse statistique dans les tests différentiels d’utilisation du site pA. Tout d’abord, nous avons supprimé les gènes avec un seul site pA, sur lequel l’utilisation différentielle du site pA ne peut pas être définie. Deuxièmement, nous appliquons un filtrage non spécifique basé sur la couverture: les comptages sont filtrés par cpm inférieur à 1 dans x sur n échantillons, où x est le nombre minimum de répétitions dans n’importe quelle condition. N = 8 et x = 2 pour cet exemple de données.
      load(file= "APA_countData.Rdata")# Skip this step if already loaded
# Non-specific filtering: Remove the pA sites not differentially expressed in the samples
countData<-countData$counts%>%as.data.frame%>% .[rowSums(edgeR::cpm(.)>1) >=2, ]
anno%<>% .[.$GeneID%in% rownames(countData), ]
# Remove genes with only 1 site and NA in geneIDs
dnsites<-anno%>%group_by(Symbol)%>%summarise(nsites=n())%>% filter(nsites>1&!is.na(Symbol))
anno<-anno%>% filter(Symbol%in%dnsites$Symbol)
countData<-countData[rownames(countData) %in%anno$GeneID, ]
4. Analyse différentielle de l’utilisation du site de polyadénylation à l’aide de pipelines DEXSeq et diffSplice.
   1. Utilisation du package DEXSeq
      REMARQUE : Comme une matrice de contraste ne peut pas être définie pour le pipeline DEXSeq, l’analyse différentielle APA de chacune des deux conditions expérimentales doit être effectuée séparément. L’analyse APA différentielle de la condition WT et de la condition DKO est effectuée à titre d’exemple pour expliquer la procédure. Reportez-vous au dossier complémentaire »APA_analysis_3PSeq_notebook. RMD« pour le flux de travail étape par étape de cette section et l’analyse APA différentielle des autres contrastes.
      1. Chargez la bibliothèque et créez un exemple de table pour définir le plan expérimental.
         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. Préparez le fichier d’information des sites pA à l’aide du package Bioconductor 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. Utilisez les comptes de lecture générés à l’étape 3.3 et les informations de site pA obtenues à l’étape précédente pour créer l’objet DEXSeq.
         # 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. Définissez la paire de contrastes en définissant les niveaux de conditions dans l’objet DEXSeq.
         dxd1$condition<- factor(dxd1$condition, levels= c("WT", "DKO"))
# The contrast pair will be "DKO - WT"
      5. Normalisation et estimation de la dispersion. Semblable aux données de séquençage de l’ARN, pour les données de séquençage d’extrémité 3', effectuez une normalisation entre les échantillons (médiane des rapports par colonne pour chaque échantillon) en utilisant la fonction 'estimateSizeFactors', et estimez la variation des données en utilisant la fonction 'estimateDispersions', puis visualisez le résultat de l’estimation de la dispersion à l’aide de la fonction 'plotDispEsts'.
         ​dxd1 %<>% estimateSizeFactors %>% estimateDispersions %T>% plotDispEsts
      6. Test différentiel d’utilisation du site pA pour chaque gène à l’aide de la fonction 'testForDEU', puis estimation du changement de pli de l’utilisation du site pA à l’aide de la fonction 'estimateExonFoldChanges'. Vérifiez les résultats à l’aide de la fonction 'DEXSeqResults' et définissez 'FDR < 10%' comme critère pour les sites pA significativement différentiels.
         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. Visualisation des résultats différentiels d’utilisation du site pA à l’aide de tracés APA différentiels générés par la fonction 'plotDEXSeq' et de tracé de volcan par la fonction '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. Utilisation du package diffSplice. Reportez-vous au fichier R Notebook supplémentaire « APA_analysis_3PSeq_notebook. Rmd » pour le flux de travail étape par étape de cette section.
      1. Définir les contrastes d’intérêt pour l’analyse différentielle de l’utilisation des pA.
         Remarque : Cette étape doit être effectuée après la construction et le traitement de l’objet DGEList, qui est inclus dans le fichier R Notebook « 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. Visualisez le résultat des contrastes d’intérêt (ici « DKO - WT ») en utilisant des tracés APA différentiels par la fonction 'plotSplice' et des tracés de volcans avec la fonction 'EnhancedVolcano'. Reportez-vous au fichier R Notebook « APA_analysis_3PSeq_notebook. Rmd » 4.2.7 - 4.2.9 pour la visualisation d’autres paires de contraste.
         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")

Résultats

Après avoir exécuté le flux de travail étape par étape ci-dessus, les résultats d’analyse AS et APA et les résultats représentatifs se présentent sous la forme de tableaux et de diagrammes de données, générés comme suit.

COMME:
Le principal résultat de l’analyse AS (tableau supplémentaire 1 pour diffSplice; Le tableau 2 pour DEXSeq) est une liste d’exons montrant une utilisation différentielle selon les conditions, et ...

Discussion

Dans cette étude, nous avons évalué des approches basées sur les exons et les événements pour détecter l’AS et l’APA dans les données de séquençage massif de l’ARN-Seq et de l’extrémité 3'. Les approches AS basées sur les exons produisent à la fois une liste d’exons exprimés différentiellement et un classement au niveau du gène ordonné en fonction de la signification statistique de l’activité globale d’épissage différentiel au niveau du gène (tableaux 1-2, 4-5). Pour l...

matériels

Name	Company	Catalog Number	Comments
Not relevent for computational study