Identyfikacja alternatywnego splicingu i poliadenylacji w danych sekwencyjnych RNA

Podsumowanie

Alternatywny splicing (AS) i alternatywna poliadenylacja (APA) rozszerzają różnorodność izoform transkryptu i ich produktów. W tym miejscu opisujemy protokoły bioinformatyczne do analizy masowych testów sekwencjonowania RNA i sekwencjonowania końca 3' w celu wykrycia i wizualizacji AS i APA różniących się w różnych warunkach eksperymentalnych.

Streszczenie

Oprócz typowej analizy RNA-Seq do pomiaru różnicowej ekspresji genów (DGE) w warunkach eksperymentalnych/biologicznych, dane sekwencyjne RNA mogą być również wykorzystane do zbadania innych złożonych mechanizmów regulacyjnych na poziomie eksonu. Alternatywny splicing i poliadenylacja odgrywają kluczową rolę w funkcjonalnej różnorodności genu, generując różne izoformy w celu regulacji ekspresji genów na poziomie potranskrypcyjnym, a ograniczenie analiz do całego poziomu genu może pominąć tę ważną warstwę regulacyjną. W tym miejscu demonstrujemy szczegółowe analizy krok po kroku w celu identyfikacji i wizualizacji różnicowego wykorzystania eksonów i miejsc poliadenylacji w różnych warunkach, przy użyciu Bioconductor oraz innych pakietów i funkcji, w tym DEXSeq, diffSplice z pakietu Limma i rMATS.

Wprowadzenie

sekwencja RNA była przez lata szeroko stosowana, zazwyczaj do szacowania różnicowej ekspresji genów i odkrywania genów¹. Ponadto można go również wykorzystać do oszacowania zróżnicowanego wykorzystania poziomu eksonów ze względu na ekspresję różnych izoform genów, przyczyniając się w ten sposób do lepszego zrozumienia regulacji genów na poziomie potranskrypcyjnym. Większość genów eukariotycznych generuje różne izoformy poprzez alternatywny splicing (AS) w celu zwiększenia różnorodności ekspresji mRNA. Zdarzenia AS można podzielić na różne wzorce: pomijanie pełnych eksonów (SE), w których ekson ("kasetowy") jest całkowicie usuwany z transkryptu wraz z jego flankującymi intronami; alternatywny (donorowy) wybór miejsca splicingu 5' (A5SS) i alternatywny wybór miejsca splicingu 3' (akceptorowego) (A3SS), gdy na obu końcach eksonu obecne są dwa lub więcej miejsc splicingu; zachowanie intronów (RI), gdy intron jest zachowany w dojrzałym transkryptie mRNA i wzajemne wykluczenie użycia eksonów (MXE), gdzie tylko jeden z dwóch dostępnych eksonów może być zachowany w tym samym czasie^2,3. Alternatywna poliadenylacja (APA) odgrywa również ważną rolę w regulacji ekspresji genów przy użyciu alternatywnych miejsc poli(A) do generowania wielu izoform mRNA z jednego transkryptu⁴. Większość miejsc poliadenylacji (pAs) znajduje się w regionie 3' bez translacji (3' UTR), generując izoformy mRNA o różnych długościach UTR 3'. Ponieważ 3' UTR jest centralnym węzłem do rozpoznawania elementów regulatorowych, różne długości 3' UTR mogą wpływać na lokalizację, stabilność i translację mRNA⁵. Istnieje klasa testów sekwencjonowania końca 3' zoptymalizowanych pod kątem wykrywania APA, które różnią się szczegółami protokołu⁶. Opisany tutaj potok jest przeznaczony dla PolyA-seq, ale może być dostosowany do innych protokołów zgodnie z opisem.

W tym badaniu prezentujemy potok metod różnicowej analizy eksonów^{7,8 (<}strong class="xfig">Rysunek 1), które można podzielić na dwie szerokie kategorie: oparte na eksonach (DEXSeq⁹, diffSplice¹⁰) i oparte na zdarzeniach (replikowane Wielowymiarowe Analizy Splicingu Transkryptów (rMATS)¹¹). Metody oparte na eksonach porównują zmianę fałdowania w różnych warunkach poszczególnych eksonów z miarą ogólnej zmiany fałdowania genu, aby wywołać użycie eksonów o zróżnicowanej ekspresji, a na tej podstawie obliczają miarę aktywności AS na poziomie genu. Metody oparte na zdarzeniach wykorzystują odczyty połączeń obejmujących ekson-intron do wykrywania i klasyfikowania określonych zdarzeń splicingu, takich jak pomijanie eksonów lub zachowanie intronów, oraz rozróżniania tych typów AS w output³. W związku z tym metody te zapewniają komplementarne widoki dla pełnej analizy AS^12,13. Do badania wybraliśmy DEXSeq (oparty na DESeq2¹⁴ DGE package) i diffSplice (oparty na pakiecie Limma¹⁰ DGE), ponieważ są to jedne z najczęściej używanych pakietów do różnicowej analizy splicingu. rMATS został wybrany jako popularna metoda analizy opartej na zdarzeniach. Inną popularną metodą opartą na zdarzeniach jest MISO (Mieszanina Izoform)¹. W przypadku APA stosujemy podejście oparte na eksonach.

Rysunek 1. Potok analizy. Schemat blokowy kroków wykorzystanych w analizie. Kroki obejmują: uzyskanie danych, przeprowadzenie kontroli jakości i wyrównanie odczytu, a następnie zliczanie odczytów przy użyciu adnotacji dla znanych eksonów, intronów i miejsc pA, filtrowanie w celu usunięcia niskich liczb i normalizację. Dane PolyA-seq analizowano pod kątem alternatywnych miejsc pA przy użyciu metod diffSplice/DEXSeq, masowe sekwencje RNA-Seq analizowano pod kątem alternatywnego splicingu na poziomie eksonów za pomocą metod diffSplice/DEXseq, a zdarzenia AS analizowano za pomocą rMATS. Kliknij tutaj, aby zobaczyć większą wersję tego rysunku.

Dane sekwencyjne RNA użyte w tym badaniu zostały pozyskane z Gene Expression Omnibus (GEO) (GSE138691)¹⁵. Wykorzystaliśmy dane sekwencyjne RNA myszy z tego badania z dwiema grupami schorzeń: nokautem typu dzikiego (WT) i nokautem typu 1 przypominającym mięśnie (Mbnl1 KO) z trzema powtórzeniami każda. Aby zademonstrować różnicową analizę użycia miejsca poliadenylacji, uzyskaliśmy dane PolyA-seq fibroblastów zarodków myszy (MEF) (GEO Accession GSE60487)¹⁶. Dane mają cztery grupy warunków: typ dziki (WT), podwójny nokaut typu 1/typ 2 (Mbnl1/2 DKO), Mbnl 1/2 DKO z powaleniem Mbnl3 (KD) i Mbnl1/2 DKO z kontrolą Mbnl3 (Ctrl). Każda grupa warunków składa się z dwóch powtórzeń.

	Przystąpienie GEO	Numer uruchomienia	Przykładowa nazwa	warunek	symulować	tkanka	Sekwencjonowanie	Długość odczytu
Sekwencja RNA	GSM4116218	SRR10261601	Mbnl1KO_Thymus_1	Nokaut Mbnl1	Przedstawiciel 1	grasica	Sparowany koniec	100 punktów bazowych
	GSM4116219	SRR10261602	Mbnl1KO_Thymus_2	Nokaut Mbnl1	Przedstawiciel 2	grasica	Sparowany koniec	100 punktów bazowych
	GSM4116220	SRR10261603	Mbnl1KO_Thymus_3	Nokaut Mbnl1	Przedstawiciel 3	grasica	Sparowany koniec	100 punktów bazowych
	GSM4116221	SRR10261604	WT_Thymus_1	Typ dziki	Przedstawiciel 1	grasica	Sparowany koniec	100 punktów bazowych
	GSM4116222	SRR10261605	WT_Thymus_2	Typ dziki	Przedstawiciel 2	grasica	Sparowany koniec	100 punktów bazowych
	GSM4116223	SRR10261606	WT_Thymus_3	Typ dziki	Przedstawiciel 3	grasica	Sparowany koniec	100 punktów bazowych
3P-Seq	GSM1480973	SRR1553129	WT_1	Typ dziki (WT)	Przedstawiciel 1	Mysie fibroblasty embrionalne (MEF)	Pojedynczy koniec	40 punktów bazowych
	GSM1480974	SRR1553130	WT_2	Typ dziki (WT)	Przedstawiciel 2	Mysie fibroblasty embrionalne (MEF)	Pojedynczy koniec	40 punktów bazowych
	GSM1480975	SRR1553131	DKO_1	Mbnl 1/2 podwójny nokaut (DKO)	Przedstawiciel 1	Mysie fibroblasty embrionalne (MEF)	Pojedynczy koniec	40 punktów bazowych
	GSM1480976	SRR1553132	DKO_2	Mbnl 1/2 podwójny nokaut (DKO)	Przedstawiciel 2	Mysie fibroblasty embrionalne (MEF)	Pojedynczy koniec	40 punktów bazowych
	GSM1480977	SRR1553133	DKOsiRNA_1	Podwójny nokaut Mbnl 1/2 z Mbnl 3 siRNA (KD)	Przedstawiciel 1	Mysie fibroblasty embrionalne (MEF)	Pojedynczy koniec	40 punktów bazowych
	GSM1480978	SRR1553134	DKOsiRNA_2	Podwójny nokaut Mbnl 1/2 z Mbnl 3 siRNA (KD)	Przedstawiciel 2	Mysie fibroblasty embrionalne (MEF)	Pojedynczy koniec	36 pz
	GSM1480979	SRR1553135	DKONTsiRNA_1	Podwójny nokaut Mbnl 1/2 z niecelującym siRNA (Ctrl)	Przedstawiciel 1	Mysie fibroblasty embrionalne (MEF)	Pojedynczy koniec	40 punktów bazowych
	GSM1480980	SRR1553136	DKONTsiRNA_2	Podwójny nokaut Mbnl 1/2 z niecelującym siRNA (Ctrl)	Przedstawiciel 2	Mysie fibroblasty embrionalne (MEF)	Pojedynczy koniec	40 punktów bazowych

Tabela 1. Podsumowanie zestawów danych RNA-Seq i PolyA-seq wykorzystanych do analizy.

Protokół

1. Instalacja narzędzi i pakietów R używanych w analizie

1. Conda to popularny i elastyczny menedżer pakietów, który umożliwia wygodną instalację pakietów wraz z ich zależnościami na wszystkich platformach. Użyj "Anaconda" (menedżer pakietów conda), aby zainstalować "conda", którego można użyć do zainstalowania narzędzi/pakietów wymaganych do analizy.
2. Pobierz program "Anaconda" zgodnie z wymaganiami systemowymi z https://www.anaconda.com/products/individual#Downloads i zainstaluj go, postępując zgodnie z instrukcjami instalatora graficznego. Zainstaluj wszystkie wymagane pakiety za pomocą polecenia "conda", wpisując następujące polecenie w wierszu poleceń systemu Linux.
   conda install -c daler sratoolkit
conda install -c conda-forge równoległy
conda install -c bioconda star bowtie fastqc rmats rmats2sashimiplot samtools fasterq-dump cutdostosuj bedtools deeptools
3. Aby pobrać wszystkie pakiety języka R używane w protokole, wpisz następujący kod w konsoli języka R (uruchamianym w wierszu polecenia systemu Linux przez wpisanie "R") lub konsoli Rstudio.
   bioc_packages<- c("DEXSeq", "Rsubread", "EnhancedVolcano", "edgeR", "limma", "maser","GenomicRanges")
pakiety<- c("magrittr", "rtracklayer", "tidyverse", "openxlsx", "BiocManager")
#Install, jeśli nie jest jeszcze zainstalowany
installed_packages<-packages%in% nazwy wierszy(installed.packages())
installed_bioc_packages<-bioc_packages%w% nazw wierszy(installed.packages())
if(any(installed_packages==FALSE)) {
install.packages(pakiety[!installed_packages],dependencies=PRAWDA)
BiocManager::install(pakiety[!installed_bioc_packages], dependencies=PRAWDA)
}
   UWAGA: W tym protokole obliczeniowym polecenia będą podawane jako pliki notesu języka R (pliki z rozszerzeniem ". Rmd"), pliki z kodem R (pliki z rozszerzeniem ". R") lub skrypty powłoki Linux Bash (pliki z rozszerzeniem ".sh"). Pliki notesu języka R (Rmd) należy otworzyć w programie RStudio przy użyciu polecenia Plik|Otwórz plik..., a poszczególne fragmenty kodu (które mogą być poleceniami języka R lub poleceniami powłoki powłoki Bash), a następnie uruchom je interaktywnie, klikając zieloną strzałkę w prawym górnym rogu. Pliki kodu języka R można uruchamiać, otwierając je w programie RStudio lub w wierszu polecenia systemu Linux, poprzedzając je ciągiem "Rscript", np. przykład Rscript.R. Skrypty powłoki są uruchamiane w wierszu polecenia systemu Linux, poprzedzając skrypt poleceniem "sh" e.g. sh example.sh.

2. Alternatywna analiza splicingu (AS) przy użyciu sekwencji RNA

1. Pobieranie i wstępne przetwarzanie danych
   UWAGA: Fragmenty kodu oznaczone poniżej są dostępne w dodatkowym pliku kodu "AS_analysis_RNASeq.Rmd", aby interaktywnie wykonywać poszczególne kroki, a także są dostarczane jako skrypt bash do uruchamiania wsadowego w wierszu poleceń systemu Linux (sh downloading_data_preprocessing.sh).
   1. Pobieranie surowych danych.
      1. Pobierz nieprzetworzone dane z archiwum odczytu sekwencyjnego () za pomocą polecenia "prefetch" z zestawu narzędzi (wersja 2.10.8)¹⁷. Podaj identyfikatory dostępu do w kolejności w poniższym poleceniu, aby pobrać je równolegle za pomocą narzędzia GNU parallel¹⁸. Aby pobrać pliki z identyfikatorami dostępu z SRR10261601 do SRR10261606 równolegle, użyj następujących poleceń w wierszu poleceń systemu Linux.
         seq 10261601 10261606 | parallel prefetch SRR{}
      2. Rozpakuj pliki fastq z archiwum za pomocą funkcji "fastq-dump" z zestawu narzędzi. Użyj GNU równolegle i podaj nazwy wszystkich plików razem.
         parallel -j 3 fastq-dump --gzip --skip-technical --read-filter pass --dumpbase --split-e --clip --origfmt {} :::
      3. Pobierz genom referencyjny i adnotacje dla myszy (zestaw genomu GRCm39) z www.ensembl.org, korzystając z następujących elementów w wierszu poleceń systemu 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 adnotacja.gtf)
GENOME=$(readlink -f genome.fa)
   2. Wstępne przetwarzanie i mapowanie odczytów do zestawu genomu
      1. Kontrola jakości. Oceń jakość surowych odczytów za pomocą FASTQC (FASTQ Quality Check v0.11.9)¹⁹. Utwórz folder wyjściowy i uruchom fastqc z równoległym na wielu wejściowych plikach fasta. Ten krok spowoduje wygenerowanie raportu jakości dla każdej próbki. Przed przystąpieniem do dalszej analizy należy sprawdzić raporty, aby upewnić się, że jakość odczytów jest akceptowalna. (Zapoznaj się z instrukcją obsługi, aby zapoznać się z raportami pod adresem https://www.bioinformatics.babraham.ac.uk/projects/fastqc/)
         mkdir fastqc_out
parallel "fastqc {} -o fastqc_out" ::: $RAW_DATA/*.fastq.gz
         UWAGA: Jeśli to konieczne, wykonaj przycinanie adaptera za pomocą 'cutadapt'²⁰ lub 'trimmomatic'²¹, aby usunąć sekwencjonowanie w adapterach flankujących, które różni się w zależności od rozmiaru fragmentu RNA i długości odczytu. W tej analizie pominęliśmy ten krok, ponieważ odsetek odczytów, których dotyczy problem, był minimalny.
      2. Odczytywanie wyrównania. Kolejnym krokiem w przetwarzaniu wstępnym jest mapowanie odczytów do genomu referencyjnego. Po pierwsze, zbuduj indeks dla genomu referencyjnego za pomocą funkcji "genomeGenerate" STAR²², a następnie wyrównaj surowe odczyty do odniesienia (alternatywnie wstępnie utworzone indeksy są dostępne na stronie internetowej STAR i mogą być używane bezpośrednio do wyrównania). Uruchom następujące polecenia w wierszu poleceń systemu Linux.
         #Build STAR indeks
GDIR=STAR_indices
mkdir $GDIR
STAR --runMode genomeGenerate --genomeFastaFiles $GENOME --sjdbGTFfile $GTF --runThreadN 8 --genomeDir $GDIR
ODIR=wyniki/mapowanie
mkdir -p $ODIR
#Align odczytuje do genomu
dla fq1 w $RAW_DATA/*R1.fastq.gz;
do
fq2=$(echo $fq 1| sed 's/1.fastq.gz/2.fastq.gz/g');
OUTPUT=$(nazwa podstawowa ${fq1}| sed 's/R1.fastq.gz//g');
STAR --genomeDir $GDIR \
--runThreadN 12 \
--readFilesCommand zcat \
--readFilesIn ${fq1}${fq2}\
--outFileNamePrefix $ODIR\/${OUTPUT} \
--outSAMtype BAM PosortowaneWedługWspółrzędnych \
--outSAMunmappped W \
--outSAMattributes Standard
Zrobione
         UWAGA: Wyrównywacz STAR wygeneruje i posortuje pliki BAM (Binary Alignment Map) dla każdej próbki po odczytanym wyrównaniu. Pliki Bam należy posortować przed przejściem do dalszych kroków.
2. Przygotowanie adnotacji eksonów.
   1. Uruchom plik kodu uzupełniającego "prepare_mm_exon_annotation. R" z pobraną adnotacją w formacie GTF (Gene transfer format) w celu przygotowania adnotacji. Aby uruchomić, wpisz następujące polecenie w wierszu poleceń systemu Linux.
      Rscript prepare_mm_exon_annotation. Adnotacja R.gtf
      UWAGA: Plik GTF zawiera wiele wpisów eksonów dla różnych izoform. Ten plik służy do "zwijania" wielu identyfikatorów transkrypcji dla każdego eksonu. Jest to ważny krok w celu zdefiniowania pojemników do liczenia eksonów.
3. Liczenie odczytów. Następnym krokiem jest policzenie liczby odczytów zmapowanych na różne transkrypty/eksony. Patrz plik uzupełniający: "AS_analysis_RNASeq.Rmd".
   1. Załaduj wymagane biblioteki:
      pakiety<- c("Rsubread","tidyverse", "magrittr", "EnhancedVolcano", "edgeR","openxlsx")
invisible(lapply(pakiety, biblioteka, character.only=PRAWDA))
   2. Załaduj przetworzony plik adnotacji uzyskany w poprzednim kroku (2.2).
      load("mm_exon_anno. RData")
   3. Odczytaj wszystkie pliki bam uzyskane w kroku 2.2.2 jako dane wejściowe dla "featureCounts" w celu zliczania odczytów. Odczytaj folder zawierający pliki bam, najpierw wyświetlając listę wszystkich plików z katalogu, który kończy się na .bam. Użyj "featureCounts" z pakietu Rsubread, który pobiera pliki bam i przetworzoną adnotację GTF (odniesienie) jako dane wejściowe, aby wygenerować macierz liczb powiązanych z każdą funkcją z wierszami reprezentującymi eksony (cechy) i kolumnami reprezentującymi próbki.
      countData <- dir("bams", pattern=".bam$", full.names=T) %>%
featureCounts(annot.ext=anno,
isGTFAnnotationFile=FAŁSZ,
minMQS=0,useMetaFeatures=FALSE,
allowMultiOverlap=PRAWDA,
largestOverlap=PRAWDA,
countMultiMappingReads=FAŁSZ,
primaryOnly=PRAWDA,
isPairedEnd=PRAWDA,
nthreads=12)
   4. Następnie należy przeprowadzić filtrowanie niespecyficzne, aby usunąć eksony o niskiej ekspresji ("niespecyficzne" oznacza, że informacje o warunkach eksperymentalnych nie są używane w filtrowaniu, aby uniknąć błędów selekcji). Przekształć dane ze skali surowej w liczbę na milion (cpm) za pomocą funkcji cpm z pakietu "edgeR" package²³ i zachowaj eksony z liczbami większymi niż ustawiony próg (dla tego zestawu danych używany jest jeden cpm) w co najmniej trzech próbkach. Usuń również geny z tylko jednym eksonem.
      # Filtrowanie niespecyficzne: Usuń eksony o niskich liczbach
isexpr<- rownames(countData$counts)[rowSums(cpm(countData$counts)>1) >=3]
countData$counts<-countData$counts[nazwy wierszy(countData$counts) %in%isexpr, ]
anno<-anno%>% filter(GeneID%w% nazwach wierszy(countData$counts))
# Usuń geny z tylko 1 miejscem i NA w 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[nazwy wierszy(countData$counts) %in%anno$GeneID, ]
      UWAGA: Sprawdź wymagane parametry dla featureCounts w przypadku korzystania z różnych danych, na przykład w przypadku odczytów jednokońcowych ustaw wartość "isPairedEnd = FALSE". Zapoznaj się z podręcznikiem użytkownika RSubread, aby wybrać opcje dla swoich danych, i zapoznaj się z sekcją Dyskusja poniżej.
4. Spicing różnicowy i analiza użycia eksonów. Opisujemy dwie alternatywy dla tego kroku: DEXSeq i DiffSplice. Każdy z nich może być używany i dawać podobne wyniki. Aby zapewnić spójność, wybierz DEXSeq, jeśli wolisz pakiet DESeq2 dla DGE i użyj DiffSplice do analizy DGE opartej na Limma. Patrz plik uzupełniający: "AS_analysis_RNASeq.Rmd".
   1. Korzystanie z pakietu DEXSeq do różnicowej analizy eksonów.
      1. Załaduj bibliotekę i utwórz przykładową tabelę, aby zdefiniować projekt eksperymentalny.
         biblioteka(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("spared-end"))))<br /> UWAGA: Nazwy wierszy powinny być spójne z nazwami plików bam używanymi przez featureCounts do zliczania odczytów. sampleTable składa się ze szczegółów każdej próbki, które obejmują: typ biblioteki i warunek. Jest to wymagane do zdefiniowania kontrastów lub grupy testowej do wykrywania użycia różnicowego.
      2. Przygotuj plik z informacjami o eksonach. Informacje eksonowe w postaci obiektów GRanges (zakresy genomu) (https://bioconductor.org/packages/release/bioc/html/GenomicRanges.html) są wymagane jako dane wejściowe do utworzenia obiektu DEXSeq w następnym kroku. Dopasuj identyfikatory genów do liczby odczytów, aby utworzyć exoninfo object.
         exoninfo<-anno[anno$GeneID%in% nazwy wierszy(countData$counts),]
exoninfo<-GRanges(seqnames=anno$Chr,
ranges=IRanges(start=anno$Start, end=anno$Koniec, width=anno$Szerokość),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. Utwórz obiekt DEXSeq za pomocą funkcji DEXSeqDataSet. Obiekt DEXSeq zbiera razem liczbę odczytów, informacje o cechach eksonów i informacje o próbce. Użyj zliczeń odczytu wygenerowanych w kroku 3 i informacji o eksonach uzyskanych w poprzednim kroku, aby utworzyć obiekt DEXSeq na podstawie macierzy zliczania. Argument sampleData przyjmuje dane wejściowe ramki danych definiujące próbki (i ich atrybuty: typ biblioteki i warunek), "projekt" używa sampleData do wygenerowania macierzy projektowej dla testów różnicowych przy użyciu notacji formuły modelu. Zauważ, że istotny termin interakcji, warunek:ekson, wskazuje, że ułamek odczytów genu przypadającego na dany ekson zależy od warunków eksperymentalnych, tj. występuje AS. Zapoznaj się z dokumentacją DEXSeq, aby uzyskać pełny opis ustawiania formuły modelu dla bardziej złożonych projektów eksperymentalnych. Aby uzyskać informacje o cechach, wymagane są identyfikatory eksonów, odpowiadające im geny i transkrypty.
         dxd<-DEXSeqDataSet(countData$counts,sampleData=sampleTable, design=~sample+exon+condition:exon,featureID=exoninfo$ExonID,groupID=exoninfo$GeneID,featureRanges=exoninfo, transkrypcje=transcripts_l)
      4. Normalizacja i oszacowanie dyspersji. Następnie wykonaj normalizację między próbkami i oszacuj wariancję danych, ze względu na szum liczby Poissona wynikający z dyskretnej natury sekwencji RNA, jak i zmienności biologicznej, używając następujących poleceń.
         dxd %<>% estimateSizeFactors %>% estimateDyspersje %T>% plotDispEsts
      5. Przetestuj użycie różnicowe. Po oszacowaniu zmienności przetestuj użycie eksonów różnicowych dla każdego genu i wygeneruj wyniki.
         dxd%<>%testForDEU%>%estimateExonFoldChanges(fitExpToVar=
"warunek")#Estimate krotne zmiany
dxr=DEXSeqWyniki(dxd)
      6. Wizualizacja zdarzeń splicingu dla wybranych genów za pomocą następującego polecenia.
         plotDEXSeq(dxr,"Wnk1", displayTranscripts=TRUE, splicing=TRUE,legenda
=PRAWDA,oś cex=1.2,cex=1.3,lwd=2)
         Zapoznaj się z plikiem notesu języka R "AS_analysis_RNASeq.Rmd", aby wygenerować dodatkowe wykresy dla genów będących przedmiotem zainteresowania i wygenerować wykresy wulkanów przy różnych progach.
   2. Używanie dyfuzji z Limmy do identyfikacji łączenia różnicowego. Postępuj zgodnie z plikiem notesu języka R "AS_analysis_RNASeq.Rmd". Upewnij się, że kroki 2.1-2.3 zostały wykonane w celu przygotowania plików wejściowych, zanim przejdziesz dalej.
      1. Załaduj biblioteki
         biblioteka(limma)
biblioteka(krawędźR)
      2. Filtrowanie niespecyficzne. Wyodrębnić macierz liczby odczytów uzyskaną w ppkt 2.3. Utwórz listę cech za pomocą funkcji 'DGEList' z pakietu edgeR, gdzie wiersze reprezentują geny, a kolumny reprezentują próbki.
         mycounts=countData$counts
#Change nazwy wierszy count data do identyfikatorów eksonów zamiast genów dla unikatowych nazw wierszy.
nazwy wierszy(mojeliczby) = exoninfo$ExonID
dge<-DGEList(liczby=mojeliczne)
#Filtering
isexpr<- rowSums(cpm(dge)>1) >=3
dge<-dge[isexpr,,keep.lib.sizes=FALSE]
#Extract adnotacje eksonów tylko dla transkrypcji spełniających niespecyficzny filter
exoninfo=anno%>% filter(ExonID%w% nazwach wierszy(dge$counts))
>#Convert exoninfo do 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$Identyfikatory transkrypcji, "\\,")
         UWAGA: Jako niespecyficzny krok filtrowania, liczby są filtrowane przez cpm < 1 in x z n próbek, gdzie x jest minimalną liczbą powtórzeń w dowolnym stanie. n = 6 i x = 3 dla tych przykładowych danych.
      3. Znormalizuj zliczenia w próbkach za pomocą funkcji "calcNormFactors" z pakietu "edgeR" przy użyciu wartości Trimmed Mean of M (metoda normalizacji TMM)²⁴ Obliczy współczynniki skalowania, aby dostosować rozmiary bibliotek.
         dge<-calcNormFactors(dge)
      4. Użyj tabeli sampleTable wygenerowanej w kroku 2.4.1.1 i utwórz macierz projektu. Matryca projektowa charakteryzuje projekt. Zobacz rozdziały 8 i 9 Podręcznika użytkownika Limma (https://www.bioconductor.org/packages/devel/bioc/vignettes/limma/inst/doc/usersguide.pdf), aby uzyskać szczegółowe informacje na temat matryc projektowych dla bardziej zaawansowanych projektów eksperymentalnych.
         Treat<- factor(sampleTable$condition)
design<- model.matrix(~0+Treat)
colnames (projekt) <- poziomy (leczyć)
      5. Dopasuj model liniowy dla każdego eksonu. Uruchom funkcję "voom" pakietu "limma", aby przetworzyć dane sekwencyjne RNA w celu oszacowania wariancji i wygenerowania precyzyjnych wag w celu skorygowania szumu liczby Poissona i przekształcenia zliczeń na poziomie eksonów na log2 zliczenia na milion (logCPM). Następnie uruchom modelowanie liniowe za pomocą funkcji "lmfit", aby dopasować modele liniowe do danych wyrażeń dla każdego eksonu. Obliczanie empirycznych statystyk Bayesa dla dopasowanego modelu przy użyciu funkcji "eBayes" w celu wykrycia różnicowego wyrażenia eksonów. Następnie zdefiniuj macierz kontrastu dla interesujących nas porównań eksperymentalnych. Użyj 'contrasts.fit', aby uzyskać współczynniki i błędy standardowe dla każdej pary porównania.
         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=projekt)
fit2<-contrasts.fit(dopasowanie;matryca ciągła)
      6. Różnicowa analiza spawów. Uruchom 'diffSplice' na dopasowanym modelu, aby przetestować różnice w użyciu eksonów genów między typem dzikim a nokautem i zbadać najwyżej oceniane wyniki za pomocą funkcji 'topSplice': test="t" daje ranking eksonów AS, test="simes" daje ranking genów.
         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. wizualizacja. Wykreśl wyniki za pomocą funkcji "plotSplice", podając interesujący gen w argumencie geneid. Zapisz najlepsze wyniki posortowane według logarytmu Złóż zmianę w obiekt i wygeneruj wykres wulkanu, aby wyświetlić eksony.
         plotSplice(ex,geneid="Wnk1", FDR=0.1)
#Volcano działka
EnhancedVolcano(ts,lab=ts$ExonID,selectLab= head((ts$ExonID),2000), xlab= bquote(~Log[2]~'zmiana złożenia'), x='logFC', y='P.Value', title='Wykres wulkanu', subtitle='Mbnl1_KO vs WT (Limma_diffSplice)', FCcutoff=2, labSize=4,legendPosition="right", caption= bquote(~Log[2]~"Fold change cutoff, 2; FDR 10%"))
   3. Korzystanie z rMATS
      1. Upewnij się, że najnowsza wersja rMATS v4.1.1 (znana również jako rMATS turbo ze względu na skrócony czas przetwarzania i mniejsze wymagania dotyczące pamięci) jest zainstalowana przy użyciu conda lub github (https://github.com/Xinglab/rmats-turbo/releases/download/v4.1.1/rmats_turbo_v4_1_1.tar.gz) w katalogu roboczym. Postępuj zgodnie z sekcją 4.3 w "AS_analysis_RNASeq.Rmd".
      2. Przejdź do folderu zawierającego pliki bam uzyskane po mapowaniu i przygotuj pliki tekstowe, zgodnie z wymaganiami rMATS, dla dwóch warunków, kopiując nazwę plików bam (wraz ze ścieżką) oddzieloną znakiem ','. W wierszu poleceń systemu Linux należy uruchomić następujące polecenia:
         mkdir rMATS_analysis
Płyta 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. Uruchom rmats.py z dwoma plikami wejściowymi wygenerowanymi w poprzednim kroku wraz z plikiem GTF uzyskanym w 2.1.1.3. Spowoduje to wygenerowanie folderu wyjściowego "rmats_out" zawierającego pliki tekstowe opisujące statystyki (wartości p i poziomy włączenia) dla każdego zdarzenia splicingu osobno.
         python rmats-turbo/rmats.py --b1 KO.txt --b2 Wt.txt --gtf adnotacja.gtf -t sparowany --readLength 50 --nwątek 8 --od rmats_out/ --tmp rmats_tmp --task pos
         UWAGA: Wymagana jest również adnotacja referencyjna w postaci pliku GTF. Sprawdź parametry, jeśli dane są jednokońcowe, i odpowiednio zmień opcję -t.
      4. Eksplorowanie wyników rMATS. Użyj pakietu Bioconductor 'maser'²⁵, aby zapoznać się z wynikami rMATS. Załaduj pliki tekstowe JCEC (Junction and Exon count) do obiektu "maser" i przefiltruj wynik na podstawie zasięgu, uwzględniając co najmniej pięć średnich odczytów na zdarzenie splicingu.
         biblioteka(maser)
mbnl1<-maser("/rmats_out/", c("WT","Mbnl1_KO"), ftype="JCEC")
#Filtering wydarzeń według zasięgu
mbnl1_filt<-filterByCoverage(mbnl1,avg_reads=5)
      5. Wizualizacja wyników rMATS. Wybierz znaczące zdarzenia splicingu przy False Discovery Rate (FDR) 10% i minimalnej 10% zmianie procentu splicingu (deltaPSI) za pomocą funkcji 'topEvents' z pakietu 'maser'. Następnie sprawdź zdarzenia genowe dla poszczególnych genów będących przedmiotem zainteresowania (gen próbki-Wnk1) i wykreśl wartości PSI dla każdego zdarzenia splicingu tego genu. Wygeneruj wykres wulkanu, określając typ zdarzenia.
         #Top zdarzenia splicingu przy 10% FDR
mbnl1_top<-topEvents(mbnl1_filt,fdr=0.1, deltaPSI=0.1)
mbnl1_top
#Check zdarzenia genowe dla konkretnego genu
mbnl1_wnk1<-geneEvents(mbnl1_filt,geneS="Wnk1", fdr=0.1, deltaPSI=0.1)
maser::d isplay(mbnl1_wnk1,"SE")
plotGenePSI(mbnl1_wnk1,type="SE", show_replicates
=PRAWDA)
wulkan(mbnl1_filt,fdr=0.1, deltaPSI=0.1,type="SE")
+xlab("deltaPSI")+ylab("Log10 Adj. Pvalue")+ggtitle("Wykres wulkanu zdarzeń pomijania eksonów")
      6. Wygeneruj wykresy Sashimi dla wyniku zdarzeń splicingu uzyskanego za pomocą rMATS w postaci plików tekstowych za pomocą pakietu 'rmats2shahimiplot'. Uruchom skrypt języka Python w wierszu poleceń systemu 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. MATY. JC.txt --l1 WT --l2 Mbnl1_KO --exon_s 1 --intron_s 5 -o .. /rMATS_analysis/rmats2shasmi_output
         UWAGA: Ten proces może być czasochłonny, ponieważ wygeneruje wykres Sashimi dla wszystkich wyników w pliku zdarzeń. Wybierz najlepsze wyniki (nazwy genów i eksony) wyświetlane przez funkcję topEvents z "maser" i zwizualizuj odpowiadający im wykres Sashimi.

3. Alternatywna analiza poliadenylacji (APA) przy użyciu sekwencjonowania końca 3'

1. Pobieranie i wstępne przetwarzanie danych
   UWAGA: Zapoznaj się z dodatkowym plikiem notesu języka R "APA_analysis_3PSeq_notebook. Rmd", aby uzyskać pełne polecenia dotyczące kroków pobierania danych i wstępnego przetwarzania, lub uruchom dodatkowy plik bash "APA_data_downloading_preprocessing.sh" w wierszu poleceń systemu Linux.
   1. Pobierz dane z wraz z identyfikatorami przystąpienia (od 1553129 do 1553136).
   2. Adaptery trymowania i odwrotne uzupełnienie, aby uzyskać sekwencję pasm sensownych.
      UWAGA: Ten krok jest specyficzny dla zastosowanego testu PolyA-seq.
   3. Mapa jest odczytywana do montażu genomu myszy za pomocą muszki aligner²⁶.
2. Przygotowanie adnotacji do witryn pA.
   UWAGA: Przetwarzanie pliku adnotacji witryny pA odbywa się najpierw przy użyciu dodatkowego pliku R Notebook "APA_analysis_3PSeq_notebook. Rmd" (2.1 - 2.6), a następnie za pomocą pliku bash "APA_annotation_preparation.sh".
   1. Pobierz adnotację o witrynach pA z bazy danych PolyASite 2.0⁶.
   2. Wybierz adnotacje lokalizacji pA, aby zachować miejsca pA 3'-nieulegające translacji (UTR), które są oznaczone jako ekson końcowy (TE) lub 1000 nt za adnotowanym eksonem końcowym (DS) do dalszej analizy.
   3. Uzyskaj szczyty lokalizacji pA. Zakotwicz w każdym miejscu rozszczepienia pA i zwizualizuj średni zasięg odczytu za pomocą bedtools i deeptools^27,28. Wyniki pokazały, że piki zmapowanych odczytów były głównie rozproszone w ciągu ~60 pz przed miejscami rozszczepienia (ryc. 5 i ryc. uzupełniający 5). W związku z tym współrzędne miejsc pA zostały rozszerzone z pliku adnotacji do 60 pz przed ich miejscami rozszczepienia. W zależności od konkretnego zastosowanego protokołu sekwencjonowania końca 3', ten krok będzie musiał zostać zoptymalizowany pod kątem testów innych niż sekwencja PolyA.
3. Zliczanie odczytów
   1. Przygotuj plik adnotacji do witryn pA.
      anno<- read.table(file= "flanking60added.pA_annotation.łóżko",
stringsAsFactors=FALSE, check.names=FALSE, header=FALSE, sep="")
colnames(anno) <- c("chrom", "chromStart", "chromEnd", "name", "score", "strand", "rep", "adnotacja", "gene_name", "gene_id")
anno<- dplyr::select(anno,nazwa,chrom, chromStart,chromEnd, pasmo,gene_id,gene_name,rep)
colnames(anno) <- c("GeneID", "Chr", "Start", "Koniec", "Pasmo", "Ensembl", "Symbol", "repID")
   2. Zastosuj wartość "featureCounts", aby uzyskać nieprzetworzone liczby. Zapisz tabelę zliczania jako plik RData "APA_countData.Rdata" do analizy APA przy użyciu różnych narzędzi.
      countData<- dir("bamfiles", pattern="sorted.bam$", full.names=TRUE) %>%
# Odczytuje wszystkie pliki bam jako dane wejściowe dla 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")
      UWAGA: Należy pamiętać, aby zmienić dowolny z parametrów wymienionych w funkcji "featureCounts". Zmodyfikuj parametr "strandSpecific", aby upewnić się, że jest on zgodny z kierunkiem sekwencjonowania użytego testu sekwencjonowania końca 3' (empirycznie, wizualizacja danych w przeglądarce genomu nad genami na niciach dodatnich i ujemnych wyjaśni to).
   3. Zastosuj niespecyficzne filtrowanie countData. Filtrowanie może znacznie poprawić niezawodność statystyczną w różnicowych testach użycia witryny pA. Po pierwsze, usunęliśmy te geny, które miały tylko jedno miejsce pA, na którym nie można określić różnicowo użycia miejsca pA. Po drugie, stosujemy niespecyficzne filtrowanie na podstawie pokrycia: liczby są filtrowane przez cpm mniejsze niż 1 do x z n próbek, gdzie x to minimalna liczba powtórzeń w dowolnym stanie. N = 8 i x = 2 dla tego przykładu data.
      load(file= "APA_countData.Rdata")# Pomiń ten krok, jeśli już załadowany
# Filtrowanie niespecyficzne: Usuń miejsca pA, które nie są wyrażone różnicowo w próbkach
countData<-countData$counts%>%as.data.frame%>% . [rowSums(krawędźR::cpm(.)>1) >=2, ]
anno%<>% . [.$GeneID%in% nazwy wierszy(liczbadanych), ]
# Usuń geny z tylko 1 miejscem i NA w geneIDs
dnsites<-anno%>%group_by(Symbol)%>%summarise(nsites=n())%>% filter(nsites>1&!is.na(Symbol))
anno<-anno%>% filter(Symbol%w%dnsites$Symbol)
countData<-countData[nazwy wierszy(countData) %in%anno$GeneID, ]
4. Różnicowa analiza wykorzystania miejsca poliadenylacji przy użyciu potoków DEXSeq i diffSplice.
   1. Korzystanie z pakietu DEXSeq
      UWAGA: Ponieważ nie można zdefiniować macierzy kontrastu dla potoku DEXSeq, różnicowa analiza APA dla każdych dwóch warunków eksperymentalnych musi być przeprowadzona oddzielnie. Różnicowa analiza APA warunku WT i warunku DKO jest wykonywana jako przykład w celu wyjaśnienia procedury. Patrz plik uzupełniający "APA_analysis_3PSeq_notebook. Rmd" dla przepływu pracy krok po kroku w tej sekcji i różnicowej analizy APA innych kontrastów.
      1. Załaduj bibliotekę i utwórz przykładową tabelę, aby zdefiniować projekt eksperymentalny.
         c("DEXSeq", "ZakresyGenomu") %>% lapply(biblioteka, character.only=PRAWDA) %>%niewidoczny
sampleTable1<- data.frame(row.names= c("WT_1","WT_2","DKO_1","DKO_2"),
warunek= c(rep("WT", 2), rep("DKO", 2)),
libType= rep("pojedynczy koniec", 4))
      2. Przygotuj plik informacyjny o miejscach pA za pomocą pakietu Bioconductor GRanges.
         # Przygotowanie obiektu GRanges do budowy obiektu DEXSeqDataSet
PASinfo <- GRanges(seqnames = anno$Chr,
zakresy = IRanges(początek = anno$Początek, koniec = anno$Koniec),pasmo = Rle(anno$Pasmo))
mcols(PASinfo)$PASID<-anno$repID
mcols(PASinfo)$GeneEns<-anno$Ensembl
mcols(PASinfo)$GeneID<-anno$Symbol
# Przygotuj nowe identyfikatory funkcji, zastąp informacje o nici literami, aby pasowały do bieżącego identyfikatora klastra pA lokacji
new.featureID <- anno$Strand %>% as.character %>% replace(. %in% "+", "F") %>% replace(. %in% "-", "R") %>% paste0(as.character(anno$repID), .)
      3. Użyj liczników odczytów wygenerowanych w kroku 3.3 i informacji o lokacji pA uzyskanych w poprzednim kroku, aby utworzyć obiekt DEXSeq.
         # Wybierz liczbę odczytów warunku WT i DKO
countData1<- dplyr::select(countData, SRR1553129.sorted.bam, SRR1553130.sorted.bam, SRR1553131.sorted.bam, SRR1553132.sorted.bam)
# Zmień nazwy kolumn countData przy użyciu przykładowych nazw w sampleTable
colnames(countData1) <- rownames(sampleTable1)
dxd1<-DEXSeqDataSet(countData=countData1,
sampleData=sampleTable1,
design=~próbka+ekson+warunek:ekson,
featureID=new.featureID,
groupID=anno$Symbol,
featureRanges=PASinfo)
      4. Zdefiniuj parę kontrastu, definiując poziomy warunków w DEXSeq object.
         dxd1$warunek<- czynnik(dxd1$warunek, poziomy= c("WT", "DKO"))
# Para kontrastowa to "DKO - WT"
      5. Normalizacja i oszacowanie dyspersji. Podobnie jak w przypadku danych sekwencjonowania RNA, w przypadku danych sekwencjonowania na końcu 3' wykonaj normalizację między próbkami (mediana proporcji w kolumnach dla każdej próbki) za pomocą funkcji "estimateSizeFactors" i oszacuj zmienność danych za pomocą funkcji "estimateDispersions", a następnie zwizualizuj wynik oszacowania dyspersji za pomocą funkcji "plotDispEsts".
         dxd1 %<>% estimateSizeFactors %>% estimateDyspersje %T>% plotDispEsts
      6. Różnicowe badanie użycia miejsca pA dla każdego genu za pomocą funkcji "testForDEU", a następnie oszacuj zmianę krotności użycia miejsca pA za pomocą funkcji "estimateExonFoldChanges". Sprawdź wyniki za pomocą funkcji "DEXSeqResults" i ustaw "FDR < 10%" jako kryterium dla znacznie zróżnicowanych miejsc pA.
         dxd1 %<>% testForDEU %>% estimateExonFoldChanges(fitExpToVar = "warunek")
dxr1 <- DEXSeqResults(dxd1)
DXR1
mcols(dxr1)$description
table(dxr1$padj<0.1) # Sprawdź liczbę miejsc różnicowych pA (FDR < 0.1)
table(tapply(dxr1$padj<0.1, dxr1$groupID, any)) # Sprawdź liczbę genów, które pokrywają się z różnicowym miejscem pA
      7. Wizualizacja różnicowych wyników wykorzystania terenu pA za pomocą różnicowych wykresów APA generowanych przez funkcję "plotDEXSeq" i wykresu wulkanicznego za pomocą funkcji "EnhancedVolcano".
         # Wybierz 100 najlepszych witryn o znaczącej różnicy pA uszeregowanych według FDR
górna różnica. PAS<- dxr1%>%as.data.frame%>%rownames_to_column%>%arrange(padj)%$%groupID[1:100]

# Zastosuj plotDEXSeq do wizualizacji użycia różniczkowej poliA
plotDEXSeq(dxr1,"S100a7a", legend=PRAWDA, wyrażenie=FAŁSZ,splicing=PRAWDA, oś cex=1.2, cex=1.3,lwd=2)

# Zastosuj perGeneQValue, aby sprawdzić najważniejsze geny z różnicowym użyciem miejsca polyA
dxr1%<>% . [!is.na(.$padj), ]
dgene<- data.frame(perGeneQValue= perGeneQValue (dxr1)) %>%rownames_to_column("groupID")

dePAS_sig1<-DXR1%>% data.frame() %>%
dplyr::select(-matches("dyspersja|statystyka|liczbaData|genomicData"))%>%
inner_join(dgene)%>%arrange(perGeneQValue)%>%distinct()%>%
filtr(padj<0.1)

# Zastosuj pakiet EnhancedVolcano, aby zobrazować użycie strony z różnicą polyA
"EnhancedVolcano"%>% lapply(library, character.only=TRUE) %>%invisible
EnhancedVolcano(dePAS_sig1, lab=dePAS_sig1$identyfikator_grupy, x='log2fold_DKO_WT',
y='pvalue',title='Działka wulkaniczna',subtitle='DKO vs WT',
FCcutoff=1,labSize=4, legendPosition="right",
caption= bquote(~Log[2]~"Odcięcie zmiany składania, 1; FDR 10%"))
   2. Korzystanie z pakietu diffSplice. Zapoznaj się z dodatkowym plikiem notesu języka R "APA_analysis_3PSeq_notebook. Rmd" dla przepływu pracy krok po kroku w tej sekcji.
      1. Definiowanie interesujących kontrastów dla różnicowej analizy użycia pA.
         UWAGA: Ten krok należy wykonać po skonstruowaniu i przetworzeniu obiektu DGEList, który znajduje się w pliku 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=projekt)
fit2<-fit%>%contrasts.fit(matryca.kontrastu)%>%eBayes
summary(decideTests(fit2))
ex<-diffSplice(fit2,geneid=anno$Symbol,exonid=new.featureID)
topSplice(ex) #Check najważniejsze wyniki z topSplice
      2. Zwizualizuj wynik interesujących nas kontrastów (tutaj "DKO - WT") za pomocą różnicowych wykresów APA według funkcji "plotSplice" i wykresów wulkanicznych z funkcją "EnhancedVolcano". Zapoznaj się z plikiem notesu języka R "APA_analysis_3PSeq_notebook. Rmd" 4.2.7 - 4.2.9 do wizualizacji innych par kontrastu.
         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")
plotSplec(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]~'zmiana złożenia'),
x='logFC', y='P.Value', title='Działka wulkanu', subtitle='DKO vs WT',
FCcutoff=1, labSize=6, legendPosition="prawo")

Wyniki

Po uruchomieniu powyższego przepływu pracy krok po kroku, wyniki analizy AS i APA oraz reprezentatywne wyniki są w formie tabel i wykresów danych, generowanych w następujący sposób.

AS:
Główny wynik analizy AS (tabela uzupełniająca 1 dla dyfuzji; Tabela 2 dla DEXSeq) zawiera listę eksonów wykazujących zróżnicowane użycie w różnych warunkach oraz listę genów wykazujących znaczącą ogólną aktywność splicingową jednego lub więcej eksonów skład...

Dyskusja

W tym badaniu oceniliśmy podejścia oparte na eksonach i zdarzeniach do wykrywania AS i APA w masowych danych sekwencjonowania RNA-Seq i końca 3'. Podejścia oparte na eksonach AS tworzą zarówno listę eksonów o zróżnicowanej ekspresji, jak i ranking na poziomie genu uporządkowany według istotności statystycznej ogólnej zróżnicowanej aktywności splicingu na poziomie genu (tabele 1-2, 4-5). W przypadku pakietu diffSplice zróżnicowane użycie określa się poprzez dopasowanie ważonych modeli liniowych na poziomi...

Materiały

Lista materiałów użytych w tym artykule

Nazwa	Firma	Numer katalogowy	Komentarze
Nie ma znaczenia dla badań