Identyfikacja alternatywnego splicingu i poliadenylacji w danych sekwencyjnych RNA

Gunjan Dixit; Ying Zheng; Brian Parker; Jiayu Wen

doi:10.3791/62636

Subskrypcja JoVE jest wymagana do oglądania tego materiału. Zaloguj się lub rozpocznij bezpłatny okres próbny.

Artykuł metodologiczny

Identyfikacja alternatywnego splicingu i poliadenylacji w danych sekwencyjnych RNA

6.2K wyświetleń

⸱

DOI:

10.3791/62636

⸱

24 czerwca 2021

Gunjan Dixit¹ , Ying Zheng¹ , Brian Parker² , Jiayu Wen¹

¹Department of Genome Sciences, The John Curtin School of Medical Research, The Australian National University, ²Department of Biology, New York University

W tym artykule

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

Metoda RNA-seq jest od lat powszechnie stosowana, zazwyczaj do szacowania różnicowej ekspresji genów oraz odkrywania nowych genów¹. Ponadto może być ona wykorzystywana do szacowania zmienności poziomu wykorzystania eksonów wynikającej z ekspresji różnych izoform genu, co przyczynia się do lepszego zrozumienia regulacji genów na poziomie potranskrypcyjnym. Większość genów eukariotycznych generuje różne izoformy poprzez alternatywny splicing (AS), aby zwiększyć różnorodność ekspresji mRNA. Zdarzenia AS można podzielić na różne wzorce: pominięcie całych eksonów (SE), w którym ekson („kasetowy”) wraz z flankującymi go intronami zostaje całkowicie usunięty z transkryptu; alternatywny wybór miejsca splicingowego 5' (donora) (A5SS) oraz alternatywny wybór miejsca splicingowego 3' (akceptora) (A3SS), gdy na dowolnym końcu eksonu obecnych jest dwa lub więcej miejsc splicingowych; retencję intronów (RI), gdy intron zostaje zachowany w dojrzałym transkrypcie mRNA, oraz wzajemne wykluczanie eksonów (MXE), w którym w danym czasie może zostać zachowany tylko jeden z dwóch dostępnych eksonów²^,³. Alternatywna poliadenylacja (APA) również odgrywa ważną rolę w regulacji ekspresji genów, wykorzystując alternatywne miejsca poli(A) do generowania wielu izoform mRNA z pojedynczego transkryptu⁴. Większość miejsc poliadenylacji (pAs) znajduje się w 3' niekodującej regionie (3' UTRs), co prowadzi do powstania izoform mRNA o różnej długości 3' UTR. Ponieważ 3' UTR jest centralnym ośrodkiem rozpoznawania elementów regulacyjnych, różna długość 3' UTR może wpływać na lokalizację, stabilność i translację mRNA⁵. Istnieje klasa testów sekwencjonowania końców 3', zoptymalizowanych pod kątem wykrywania APA, które różnią się szczegółami protokołu⁶. Opisany tutaj potok przetwarzania danych został zaprojektowany dla PolyA-seq, ale może zostać dostosowany do innych protokołów zgodnie z opisem.

W niniejszym badaniu przedstawiamy potok metod analizy różnicowej eksonów⁷^,⁸ (Rysunek 1), które można podzielić na dwie szerokie kategorie: metody oparte na eksonach (DEXSeq⁹, diffSplice¹⁰) oraz metody oparte na zdarzeniach (replicate Multivariate Analysis of Transcript Splicing (rMATS)¹¹). Metody oparte na eksonach porównują krotność zmiany (fold change) poszczególnych eksonów między warunkami z miarą ogólnej krotności zmiany genu, aby zidentyfikować różnicowo wyrażone wykorzystanie eksonów, a następnie na tej podstawie obliczają miarę aktywności AS na poziomie genu. Metody oparte na zdarzeniach wykorzystują odczyty z połączeń ekson-intron do wykrywania i klasyfikowania specyficznych zdarzeń splicingowych, takich jak pominięcie eksonu lub retencja intronów, oraz rozróżniają te typy AS w wynikach³. W ten sposób metody te zapewniają komplementarne spojrzenie dla pełnej analizy AS¹²^,¹³. Do badania wybraliśmy DEXSeq (oparty na pakiecie DGE DESeq2¹⁴) oraz diffSplice (oparty na pakiecie DGE Limma¹⁰), ponieważ należą one do najczęściej stosowanych pakietów w analizie różnicowego splicingu. rMATS wybrano jako popularną metodę analizy opartej na zdarzeniach. Inną popularną metodą opartą na zdarzeniach jest MISO (Mixture of Isoforms)¹. W przypadku APA dostosowaliśmy podejście oparte na eksonach.

Schemat procesu analizy alternatywnego splicingu/poliadenylacji, RNA-Seq, przetwarzanie danych, analiza.
Rycina 1. Potok analizy. Schemat blokowy kroków wykorzystanych w analizie. Kroki obejmują: pozyskanie danych, przeprowadzenie kontroli jakości i mapowanie odczytów, a następnie zliczanie odczytów z wykorzystaniem adnotacji dla znanych egzonów, intronów i miejsc pA, filtrowanie w celu usunięcia niskich wartości zliczeń oraz normalizację. Dane PolyA-seq zostały przeanalizowane pod kątem alternatywnych miejsc pA przy użyciu metod diffSplice/DEXSeq, dane bulk RNA-Seq przeanalizowano pod kątem alternatywnego splicingu na poziomie egzonów metodami diffSplice/DEXseq, a zdarzenia AS przeanalizowano za pomocą rMATS. Kliknij tutaj, aby wyświetlić powiększoną wersję tej ryciny.

Dane RNA-seq wykorzystane w niniejszym badaniu zostały pobrane z bazy Gene Expression Omnibus (GEO) (GSE138691)¹⁵. Wykorzystano dane RNA-seq myszy z tego badania, obejmujące dwie grupy warunków: typ dziki (WT) oraz knockout Muscleblind-like type 1 (Mbnl1 KO), z trzema powtórzeniami dla każdej z grup. Aby zaprezentować analizę różnicowego wykorzystania miejsc poliadenylacji, pozyskano dane PolyA-seq z fibroblastów embrionalnych myszy (MEFs) (numer dostępu GEO GSE60487)¹⁶. Dane te obejmują cztery grupy warunków: typ dziki (WT), podwójny knockout Muscleblind-like type 1/type 2 (Mbnl1/2 DKO), Mbnl1/2 DKO z wyciszeniem Mbnl3 (KD) oraz Mbnl1/2 DKO z kontrolą Mbnl3 (Ctrl). Każda grupa warunków składa się z dwóch powtórzeń.

	Numer akcesyjny GEO	Numer serii SRA	Nazwa próbki	Warunki	Powtórzyć	Tkanka	Sekwencjonowanie	Długość odczytu
Sekwencjonowanie RNA (RNA-Seq)	GSM4116218	SRR10261601	Mbnl1KO_grasica_1	knockout genu Mbnl1	Powtórzenie 1	Grasica	Parne końce	100 bp
	GSM4116219	SRR10261602	Mbnl1KO_Gruczoł żylny_2	Nokaut Mbnl1	Powtórzenie 2	Grasica	odczyty parzyste	100 bp
	GSM4116220	SRR10261603	Mbnl1KO_Thymus_3	Nokaut genu Mbnl1	Powtórzenie 3	grasica	Odczyty parzyste	100 bp
	GSM4116221	SRR10261604	WT_Thymus_1	Typ dziki	Powtórzenie 1	grasica	Odczyty parzyste	100 bp
	GSM4116222	SRR10261605	WT_grasica_2	Typ dziki	Powtórzenie 2	grasica	Odczyty parzyste	100 bp
	GSM4116223	SRR10261606	WT_Thymus_3	Typ dziki	Powtórzenie 3	grasica	Odczyty parzyste	100 bp
3P-Seq	GSM1480973	SRR1553129	WT_1	Typ dziki (WT)	Powtórzenie 1	embrionalne fibroblasty mysie (MEF)	Jednostronny	40 pz
	GSM1480974	SRR1553130	WT_2	Typ dziki (WT)	Powtórzenie 2	embrionalne fibroblasty mysie (MEF)	Jednokońcowy	40 bp
	GSM1480975	SRR1553131	DKO_1	podwójny knockout (DKO) Mbnl 1/2	Powtórzenie 1	embrionalne fibroblasty mysie (MEFs)	Jednostronne	40 bp
	GSM1480976	SRR1553132	DKO_2	podwójny knockout (DKO) Mbnl 1/2	Powtórzenie 2	embrionalne fibroblasty mysie (MEF)	Jednoniciowy	40 bp
	GSM1480977	SRR1553133	DKOsiRNA_1	Podwójny nokaut Mbnl 1/2 z siRNA dla Mbnl 3 (KD)	Powtórzenie 1	embrionalne fibroblasty mysie (MEF)	Jednostronny	40 bp
	GSM1480978	SRR1553134	DKOsiRNA_2	Podwójny knockout Mbnl 1/2 z siRNA przeciwko Mbnl 3 (KD)	Powtórzenie 2	embrionalne fibroblasty mysie (MEFs)	Jednostronne	36 pz
	GSM1480979	SRR1553135	DKONTsiRNA_1	Podwójny nokaut Mbnl 1/2 z niecelującym siRNA (Ctrl)	Powtórzenie 1	embrionalne fibroblasty mysie (MEFs)	Jednostronny	40 bp
	GSM1480980	SRR1553136	DKONTsiRNA_2	Podwójny knockout Mbnl 1/2 z zastosowaniem niecelującego siRNA (Ctrl)	Powtórzenie 2	embrionalne fibroblasty mysie (MEFs)	Jednostronna	40 bp

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

Dostęp ograniczony. Zaloguj się lub rozpocznij wersję próbną, aby wyświetlić tę treść.

Protokół

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

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.
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
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

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.
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.
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.
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'

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²⁶.
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²⁷^,²⁸. 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.
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, ]
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")

Dostęp ograniczony. Zaloguj się lub rozpocznij wersję próbną, aby wyświetlić tę treść.

Wyniki

Po przeprowadzeniu powyższego procesu krok po kroku, wyniki analizy AS oraz APA oraz reprezentatywne rezultaty są przedstawione w formie tabel i wykresów danych, wygenerowanych w następujący sposób.

AS:
Głównym wynikiem analizy AS (Supplementary Table 1 dla diffSplice; Table 2 dla DEXSeq) jest lista eksonów wykazujących różnicowe wykorzystanie w różnych warunkach oraz lista genów wykazujących istotną ogólną aktywność splicingu jednego lub ...

Dostęp ograniczony. Zaloguj się lub rozpocznij wersję próbną, aby wyświetlić tę treść.

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...

Dostęp ograniczony. Zaloguj się lub rozpocznij wersję próbną, aby wyświetlić tę treść.

Oświadczenia

Autorzy nie mają nic do ujawnienia.

Podziękowania

To badanie było wspierane przez Australian Research Council (ARC) Future Fellowship (FT16010043) oraz ANU Futures Scheme.

Dostęp ograniczony. Zaloguj się lub rozpocznij wersję próbną, aby wyświetlić tę treść.

Materiały

Lista materiałów użytych w tym artykule
Nazwa	Firma	Numer katalogowy	Komentarze
Nie ma znaczenia dla badań

obliczeniowych

Bibliografia

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

Dostęp ograniczony. Zaloguj się lub rozpocznij wersję próbną, aby wyświetlić tę treść.

Przedruki i uprawnienia

Tagi

Analiza poliadenylacji r nicowe wykorzystanie ekson w detekcja zdarze splicingu pakiety Bioconductor analiza DEXSeq Limma diffSplice przep yw pracy rMATS wykres Sashimi

Identyfikacja alternatywnego splicingu i poliadenylacji w danych sekwencyjnych RNA