Rup (RNA-seq Benchmark Assessment Pipeline) – Kontrola jakości eksperymentów masowych RNA-seq u eukariontów

Podsumowanie

Protokół ten umożliwia początkową kontrolę jakości eksperymentów RNA-seq dla biologów laboratoryjnych z ograniczonym doświadczeniem w bioinformatyce.

Streszczenie

Nowoczesne podejścia w molekularnej nauce o roślinach często wymagają eksperymentów RNA-seq w całości, na przykład do śledzenia globalnych zmian w transkryptomach po zabiegach lub identyfikacji kluczowych elementów szlaków regulacyjnych. W konsekwencji różnorodne dziedziny nauk o roślinach opierają się na wysokiej jakości i powtarzalnych danych RNA-seq w zakresie rozwoju nauki. Jednak z naszego doświadczenia wynika, że wiedza o zakresie i stosowanie miar kontroli jakości w zbiorach danych RNA-seq często są niewystarczające. Wprowadzamy tutaj Rup (RNA-seq assessment assessment pipeline) do kontroli jakości danych RNA-seq objętości, do późniejszych analiz ekspresji genów, które są samodzielne i łatwo przydatne dla biologów laboratoryjnych z podstawową wiedzą o R. Rup pomaga rozróżnić dane sekwencjonowania wysokiej jakości, odpowiednie do dalszych eksperymentów ekspresji genów, od tych, które nie nadają się do dalszej analizy. Rup obejmuje testy na kilka często spotykanych problemów, takich jak niewystarczające liczby odczytów lub mapowania, identyfikacja zanieczyszczeń, ilościowa identyfikacja frakcji rRNA w całkowitych danych RNA-seq, testy podobieństwa replikacyjnego, wykorzystanie rzeczywistych danych do demonstracji oraz oferowanie intuicyjnej wizualizacji. Rup oferuje zestaw narzędzi do identyfikacji niedociągnięć eksperymentalnych przed standaryzowaną analizą transkryptomową, co poprawia jakość danych zarówno dla poszczególnych badaczy, jak i dla całej dziedziny. Zwiększa to zaufanie do masowej analizy danych RNA-seq i stanowi podstawę dla przyszłych wytycznych definiujących minimalne kryteria kontroli jakości, poprawiając tym samym wiarygodność i przejrzystość opublikowanych danych RNA-seq. Dane z rozbicia i testów są dostępne w https://github.com/oliverrupp/rup.

Wprowadzenie

Eksperymenty transkryptomiki (RNA-seq) kompleksowo badają sygnatury transkrypcyjne kształtujące fenotypy. To podejście stało się niezastąpione w genetyce molekularnej roślin, aby identyfikować pojedyncze lub współekspresyjne geny oraz procesy biologiczne zaangażowane na przykład w rozwoju, interakcjach patogenów roślin czy odporności na stres abiotyczny¹, ², ³. Najnowsze osiągnięcia w technologii RNA-seq zwiększyły specyficzność i umożliwiły wykrywanie różnych izoform i wariantów w rozdzielczości pojedynczej zasady, umożliwiając identyfikację wariacji sekwencji od większych indel do polimorfizmów pojedynczych nukleotydów (SNP). Dane uzyskane przez RNA-seq charakteryzują się szerokim zakresem dynamicznym, ułatwiającym wykrywanie zarówno bardzo obfitych, jak i słabo eksprimowanych transkryptów, i wymagają odpowiednich warunków eksperymentalnych dla spójności. Ponadto zbiory danych RNA-seq są duże, a więc wymagające analizy i przechowywania obliczeniowo, co wymaga efektywnego zarządzania danymi i rozległych zasobów obliczeniowych. RNA-seq wymaga wysokiej jakości danych wejściowych na każdym etapie pracy, ponieważ słabości na każdym etapie mogą się propagować i zaburzać wyniki. Słaba integralność RNA lub problemy techniczne podczas przygotowań bibliotecznych prowadzą do stronniczości i obniżenia dokładności. Parametry uzyskiwania wysokiej jakości surowych odczytów mogą się różnić w zależności od obiektów sekwencjonowania nowej generacji (NGS).

Według naszego doświadczenia zalecamy używanie wyłącznie RNA o liczbie integralności RNA (RIN) powyżej 7, co wskazuje na w dużej mierze nienaruszoną strukturę mRNA, jako wejścia do przygotowania biblioteki RNA-seq. Do pomyślnego sekwencjonowania wymagane jest około 2 μg całkowitego RNA przy stężeniu 50–200 ng/μL dla standardowych protokołów przygotowania biblioteki. Czystość RNA powinna być potwierdzona przez stosunek OD_260/280 między 1,8 a 2,1 oraz OD_260/230 większy niż 1,5 za pomocą spektrofotometru. Niewystarczająca głębokość sekwencjonowania, niskie tempo mapowania lub niedopasowanie do genomu referencyjnego mogą dodatkowo zakłócać ekspresję genów i kwantyfikację splicingu. Ponadto niewystarczające projektowanie eksperymentalne, takie jak brak rozróżnienia między próbkami różnych tkanek lub zabiegów oraz duża zmienność między replikacjami, może wprowadzać szum i zmniejszać powtarzalność. Chociaż wiele laboratoriów rutynowo analizuje transkryptomy, rygorystyczne kontrole jakości pierwszych kluczowych kroków analiz często nie są podawane lub mogą być całkowicie nieobecne w publikacjach. Może to prowadzić do nadinterpretacji wyników uzyskanych na podstawie analizy transkryptomu, a w konsekwencji do wyników niepowtarzalnych.

Przedstawiamy tutaj workflow do kontroli jakości początkowych kroków wymaganych do wysokiej jakości analiz transkryptomu profili mRNA w celu pomiaru zmian transkrypcyjnych. Naszym celem jest umożliwienie biologom laboratoryjnym z ograniczoną wiedzą z zakresu bioinformatyki oceny ich podstawowych danych transkryptomicznych. Rup (Rysunek 1) jest dostępny dla badaczy zaznajomionych z podstawową wiedzą o R. Wykonanie przedstawionego tutaj workflow pozwoli badaczom uzyskać szczegółowe zrozumienie ich danych pierwotnych, w tym potencjalnych ograniczeń do dalszej analizy. Z naszej wiedzy brakuje praktycznego procesu oceny danych pierwotnych transkryptomów w połączeniu z wytycznymi pozwalającymi rozróżnić dane wysokiej od niskiej jakości.

Rysunek 1: Przepływ pracy w potoku kontroli jakości RNA-seq in silico . Pliki wejściowe do kontroli jakości pochodzą z danych RNA-seq wygenerowanych podczas sekwencjonowania materiałów zakładowych, a także z publicznie dostępnych zbiorów danych. Dostarczony potok R ocenia jakość sekwencji za pomocą trzech głównych podejść: jakości sekwencjonowania, jakości mapowania oraz jakości replikacji. Obliczane są różne metryki jakości, a wyniki statystyczne wizualizowane są za pomocą popularnych pakietów R, takich jak ggplot2 i pheatmap (np. wykresy słupkowe i mapy ciepła). Proszę kliknąć tutaj, aby zobaczyć większą wersję tej figurki.

Wcześniej zgłaszane potoki ewaluacyjne wymagały wstępnie przetworzonych danych (np. odczyt wyrównania jako pliki .bam), nie obejmują wszystkich wskaźników lub nie są już utrzymywane^{w 4,5,6}. Zaletą przedstawionego tutaj przepływu pracy jest jego kompleksowość poprzez konsolidację najczęstszych problemów kontroli jakości w jednym procesie. Ponadto podajemy przykłady wysokiej jakości danych odpowiednich do wszystkich zastosowań analizy downstream, ale także przykłady danych niskiej jakości, oraz omawiamy ich konkretne ograniczenia w dalszej analizie. Kilka wcześniej opublikowanych raportów opisuje cel narzędzi analizy RNA-seq oraz przedstawia porównawcze oceny ich wydajności ^7,8. Jednak standardy kontroli jakości RNA-seq i raportowania metodologicznego nie są standaryzowane i pogarszają powtarzalność oraz biologicznie istotne interpretacje eksperymentów transkryptomicznych.

Rup integruje standardowe, wysokiej jakości narzędzia, analizę kontroli jakości oraz wizualizację wyników. Rup wymaga surowych odczytów sekwencjonowania oraz annotowanego genomu jako wejścia i działa na Mac OS, Linux oraz "Windows Subsystem for Linux" (WSL) na systemach Windows. Integruje narzędzia do oceny jakości sekwencjonowania i ilościowo określa odkładanie odczytu do genomu, uwzględnia pomiar zawartości rRNA w próbkach oraz oferuje korelację próbek do oceny korelacji replikacyjnych. Dane testowe zostały uzyskane za pomocą laserowej mikrosekcji centralnej tkanki merysty dwóch różnych stadiów gatunku rośliny Eschscholzia californica, protokołu ultra-niskowejściowego do przygotowania biblioteki, i zsekwencjonowanego na Novaseq 6000.

Protokół

Rup może być uruchamiany jako pojedynczy skrypt z minimalnym wymaganiem plików wejściowych. Potok, jak pokazano tutaj, wymaga danych sekwencjonowania RNA na parowanych końcach Illumina. Zmodyfikowany potok do sekwencjonowania jednokierunkowego z tymi samymi krokami analizy jest również deponowany w repozytorium GitHub. Potrzebne są jedynie sekwencje genomu jako plik fasta, model genu i adnotacje rRNA jako pliki gtf oraz surowe odczyty sekwencjonowania jako fastq.gz plików. Upewnij się, że pliki genomu i annotacji są dostępne jako genome.fa, annotation.gtf oraz rRNA.gtf w folderze określonym w zmiennej reference_folder. Gdy wszystkie pliki sekwencjonowania jako .fq.gz pliki są zapisane w zmiennej read_file_folder, Rup można uruchomić w następujący sposób:

1. Przygotowania:

1. Zainstaluj wszystkie niezbędne pakiety R za pomocą Bioconductor.

   BiocManager::install(c("getopt", "ggplot2", "reshape2", "pheatmap", "fastqcr", "Rfastp", "Rsubread", "Rsamtools"))

2. Przygotuj wszystkie potrzebne pliki.
   1. Stwórz folder źródłowy, który będzie zawierał wszystkie potrzebne pliki. Dodaj sekwencję genomu referencyjnego w pliku fasta o nazwie "reference/genome.fa" do folderu źródłowego. Dostarcz adnotację modelu genu jako plik gtf: "reference/annotation.gtf", opcjonalnie podaj adnotację genu rRNA jako plik gtf: "reference/rRNA.gtf".
   2. Dodaj wszystkie odczyty sekwencjonowania jako pliki fastq do folderu "reads". Upewnij się, że wszystkie nazwy plików fastq podążają za tym samym wzorcem; _1.fastq.gz oraz _2.fastq.gz dla odczytów do przodu i do tyłu.

      # define source folder and all derived subfolders for input and output files
source_folder <- "data"
reference_folder <- file.path(source_folder, "reference")
read_file_folder <- file.path(source_folder, "reads")
results_folder <- file.path(source_folder, "results")

# define input reference files
genome_fasta_file <- file.path(reference_folder, "genome.fa")
annotation_file <- file.path(reference_folder, "annotation.gtf")
rrna_file <- file.path(reference_folder, "rRNA.gtf")

# define result subfolders
fastqc_folder <- file.path(results_folder, "fastqc")
trimmed_fastqc_folder <- file.path(results_folder, "trimmed_fastqc")
trimmed_read_folder <- file.path(results_folder, "trimmed")
bam_folder <- file.path(results_folder, "bam")
sorted_bam_folder <- file.path(results_folder, "sorted_bam")
counts_folder <- file.path(results_folder, "counts")

# create results folders
for(folder in c(results_folder, fastqc_folder, trimmed_fastqc_folder, trimmed_read_folder, bam_folder, sorted_bam_folder, counts_folder)) {
 if(!dir.exists(folder)) {
  dir.create(folder)
 }
}

# get sample prefixes from input fastq files
fastq_files <- list.files(read_file_folder, pattern = "*_1.f(ast)?q.gz")
sample_prefixes <- gsub("_1.f(ast)?q.gz", "", fastq_files)

3. Ustaw parametr zgodnie z zasobami obliczeniowymi.

   n_threads <- 8 # the number of available CPU cores
bamSortMemory <- "1024" # maximum memory for bam file sorting

4. Ustaw parametr zgodnie z metodą sekwencjonowania.
    

   MinReadLength <- 25 # minimum read length, should not be less than 25
minFragLength <- 0 # minimum fragment length distribution
maxFragLength <- 300 # maximum fragment length distribution
orientation <- "fr" # read orientation for read mapping ("fr", "rf", "ff")
stranded <- 0 # stranded sequencing
         # (0 (unstranded), 1 (stranded) and 2 (reversely stranded))

   
   
   UWAGA: Minimalna długość odczytu nie powinna być mniejsza niż 25 bp, ponieważ mniejsze odczyty mogą powodować przerwanie mapowania. Większe wartości skutkują bardziej unikalnie mapowanymi odczytami, ale także większą liczbą odrzuconych odczytów podczas przycinania, w zależności od jakości sekwencjonowania. Więcej szczegółów na temat tych parametrów opisano w suplementach.

2. Sekwencjonowanie oceny jakości

UWAGA: Ten krok utworzy wykres słupkowy pokazujący liczbę odczytów przed i po przycięciu każdej próbki.

1. Uruchom fastqc⁹, narzędzie do ogólnej kontroli jakości danych sekwencji o wysokiej przepustowości, dostarczające informacji o łącznej liczbie sekwencjonowanych odczytów, długości odczytu, średniej zawartości GC, zanieczyszczeniu adapterów sekwencjonujących oraz sekwencjach nadreprezentowanych (np. odczyty rRNA). Analizuj jakość sekwencjonowania według raportów bazowych i sekwencyjnych. Uruchom fastqc na surowych plikach sekwencjonowania za pomocą pakietu R fastqcr¹⁰ i umieść wszystkie pliki wejściowe w tym samym katalogu (read_file_folder). Podsumowanie plików wyjściowych z fastqc_folder do obiektu R za pomocą funkcji qc_aggregate.
   UWAGA: Wyniki zostaną zapisane w katalogu fastqc_folder. Pliki .html utworzone przez fastqc w folderze wyjściowym są dostępne dla dodatkowych statystyk.
    

   # load the fastqc library
library(fastqcr)
# run fastqc on all raw fastq files
fastqc(fq.dir=read_file_folder, qc.dir=fastqc_folder, threads=n_threads)

# aggregate the fastqc statistics
qc <- qc_aggregate(fastqc_folder, progress=F)
qc$tot.seq <- as.numeric(qc$tot.seq)
# remove suffixes from sample names
qc$sample = gsub(".f(ast)?q.gz", "" , qc$sample)

2. Gdy wszystkie prefiksy próbek są zapisane w sample_prefixes wektorowym, przejdź przez prefiksy i uruchom narzędzie do przycinania jakości fastp na wszystkich próbkach. Cut, aby usunąć sekwencje adapterów, primery, niskiej jakości podstawy, sekwencje poly-A/T, jest wykonywany przez Fastp¹¹, dostępny w pakiecie R Rfastp¹². Przycięte pliki fastq będą przechowywane w folderze trimmed_read_folder.
    

   # load the rfastp library
library(Rfastp)

# iterate over sample prefixes
for(prefix in sample_prefixes) {
 # create output prefix
 # !!! Rfastp automatically adds _R1.fastq.gz and _R2.fastq.gz to the prefix
outputPrefix <- file.path(trimmed_read_folder, prefix)

# get input reads based on sample prefix
read1 = file.path(read_file_folder, paste(prefix, "_1.fastq.gz", sep=""))
read2 = file.path(read_file_folder, paste(prefix, "_2.fastq.gz", sep=""))
if(!file.exists(read1)) { read1 = file.path(read_file_folder, paste(prefix, "_1.fq.gz", sep="")) }
if(!file.exists(read2)) { read2 = file.path(read_file_folder, paste(prefix, "_2.fq.gz", sep="")) }

# run fastp trimmig with minimum read length
fastp_stats <- rfastp(read1 = read1,
     read2 = read2,
     minReadLength = minReadLength,
     outputFastq = outputPrefix,
     thread = n_threads)
}

3. Ponownie uruchom fastqc na przyciętych odczytach.
    

   # run fastqc on the trimmed reads
fastqc(fq.dir=trimmed_read_folder, qc.dir=trimmed_fastqc_folder, threads=n_threads)

# aggreagate the fastqc statistics
qc_trimmed <- qc_aggregate(trimmed_fastqc_folder, progress=F)
qc_trimmed$tot.seq <- as.numeric(qc_trimmed$tot.seq)
# correct sample names (change the fastp "R1","R2" suffix to "1","2")
qc_trimmed$sample = gsub("_R([12])$", "_\\1" , qc_trimmed$sample)

4. Oceń przycinanie, zbierając liczbę odczytów przed i po przycinaniu dla wszystkich próbek. Stwórz wykres słupkowy z liczbami odczytu, używając ggplot2¹³.
    

   # load the ggplot2 library for plotting
library(ggplot2)

# add the trimming status to the fastqc results
qc$Trimming = "raw"
qc_trimmed$Trimming = "trimmed"

# combine the results into one vector
qc_all = rbind(qc, qc_trimmed)

# plot the read number before and after trimming as barplot
read_number_plot <- ggplot(qc_all, aes(x=sample, y=tot.seq, fill=Trimming)) +
geom_bar(stat="identity", position = "dodge") +
xlab("Samples") + ylab("Number of Reads") +
theme(axis.text.x = element_text(angle = 90, hjust = 0)) +
theme(text = element_text(size = 18)) +
ggtitle("Number of reads before and afer trimming")
print(read_number_plot)

Rysunek 2: Wyniki sekwencjonowania i przycinania. Odczyt liczy się przed (czerwony) i po przycięciu (zielony). Próbka s2_r1 już przed przycięciem ma niską liczbę odczytów, podczas gdy przycinanie usuwało dużą część odczytów próbki s2_r2. Wszystkie pozostałe próbki wykazują akceptowalną utratę odczytów spowodowaną przycinaniem. Proszę kliknąć tutaj, aby zobaczyć większą wersję tej figurki.

3. Jakość mapowania

UWAGA: Ten krok oblicza wykresy słupkowe, aby rozróżnić odczyty pojedynczo mapowane (np. transkrypcje genów kodujących białka) od wielokrotnie mapowanych odczytów (np. odczytów rRNA) oraz niemapowanych odczytów (np. zanieczyszczenia).

1. Użyj pakietu Rsubread¹⁴ , aby odwzorować odczyty na genom referencyjny. Najpierw zbuduj indeks referencyjnego pliku genomu fasta (genome_fasta_file).
   UWAGA: Ten krok wykonuje się raz dla każdego genomu referencyjnego.
    

   # load the Rsubread library for read mapping and read counting
library(Rsubread)

# create a subread index from the reference genome sequence
buildindex(file.path(reference_folder,"subread.index"), genome_fasta_file)

2. Iteruj wszystkie próbki i wyrównuj odczyty z genomem referencyjnym, korzystając z funkcji align() pakietu Rsubread. Ta funkcja tworzy pliki .bam w folderze wyjściowym.
    

   # load the Rsamtools library for bam file sorting and indexing
library(Rsamtools)

# iterate over sample names
for(prefix in sample_prefixes) {
# create the bam output file name
output_bam = file.path(bam_folder, paste(prefix, ".bam", sep=""))

# align the reads to the reference genome index
mapping_stats <- align(index=file.path(reference_folder,"subread.index"),
     # create the forward and reverse read file names based on the prefix
     # "_R1.fastq.gz" and "_R2.fastq.gz" are forced by fastp
     readfile1=file.path(trimmed_read_folder, paste(prefix, "_R1.fastq.gz", sep="")),
     readfile2=file.path(trimmed_read_folder, paste(prefix, "_R2.fastq.gz", sep="")),
     output_file=output_bam, # set the output file name
     type=0, # the reads are RNA-seq
     minFragLength=minFragLength, # the minimal allowed fragment length
     maxFragLength=maxFragLength, # the minimal allowed fragment length
     PE_orientation=orientation, # read orientation
     nthreads=n_threads,
     # use a GTF annotation file to support the mapping
     useAnnotation=TRUE,
     annot.ext=annotation_file,
     isGTF=TRUE,
     nBestLocations=2) # to distinguish between unique and multi-mapping reads

# create the sorted output file name (the .bam suffix will be added automatically)
output_sorted_bam = file.path(sorted_bam_folder, prefix)

# sort and index the bam file
sortBam(output_bam, output_sorted_bam, maxMemory=bamSortMemory, nThreads=n_threads)
indexBam(paste0(output_sorted_bam, ".bam"))
}

3. Policz odczyty dla każdego genu za pomocą funkcji featureCounts() z pakietu Rsubread. Upewnij się, że plik adnotacji (annotation_file) jest w formacie GTF. Count odczytuje tylko pojedyncze dopasowanie do genomu.
    

   # collect all bam files
bam_files <- list.files(bam_folder, pattern = "*.bam$")

# count the number of reads for each gene
     gene_feature_counts = featureCounts(file.path(bam_folder, bam_files),
     annot.ext = annotation_file, # the gene models
     isGTFAnnotationFile = TRUE,
     countMultiMappingReads = FALSE, # only count unique reads
     strandSpecific = stranded, # for strand specific data
     isPairedEnd = TRUE, # paired-end sequencing
# the following paramter allow to control for false-positive read assignments
     requireBothEndsMapped = TRUE,
     checkFragLength = TRUE,
     minFragLength = minFragLength,
     maxFragLength = maxFragLength,
     nthreads = n_threads)

4. Policz odczyty mapujące na geny rRNA w celu oceny zawartości rRNA w próbkach. Pozwól na liczenie wielomapowych odczytów tutaj. Określ lokusy genu rRNA w pliku GTF (rrna_file).
    

   # count the reads mapping to rRNA genes
rrna_feature_counts = featureCounts(file.path(bam_folder, bam_files),
     annot.ext = rrna_file, # the rRNA gene model file
     isGTFAnnotationFile = TRUE,
     countMultiMappingReads = TRUE, # rRNA reads usually map to multiple loci
     fraction = TRUE,
     strandSpecific = stranded,
     isPairedEnd = TRUE,
     nthreads = n_threads)

5. Zbieraj statystyki przypisania odczytu utworzone przez featureCounts. Statystyki obejmują liczbę wyrównań dla każdej próbki w różnych kategoriach (np. Przypisane, Niemapowane, MultiMapped).
    

   # load the reshape2 library for data restructuring
library(reshape2)

# get the assignement stats
gene_count_stats <- gene_feature_counts$stat
# set the sample names as column names
colnames(gene_feature_counts$counts) = gsub(".bam", "", colnames(gene_feature_counts$counts))

# set the assignment type as row names
rownames(gene_count_stats) <- gene_count_stats$Status
# select the columns with the assignment statistics
gene_count_stats <- gene_count_stats[,seq(2, ncol(gene_count_stats))]
# set the sample names as column names
colnames(gene_count_stats) <- gsub(".bam", "", colnames(gene_count_stats))

# transform the dataframe for plotting
transformed_stats <- melt(t(gene_count_stats))

# set the column names
colnames(transformed_stats) <- c("Sample", "Group", "Alignments")

# adjust the groupa and sample ordering for plotting
transformed_stats$Group <- factor(transformed_stats$Group,
          levels = rev(levels(transformed_stats$Group)[order(levels(transformed_stats$Group))]))
transformed_stats$Sample <- factor(transformed_stats$Sample,
          levels = rev(levels(transformed_stats$Sample)[order(as.character(transformed_stats$Sample))]))

# remove assignment classes without any alignments
transformed_stats <- transformed_stats[transformed_stats$Alignments > 0,]
# the statistics refer to proteins coding genes
transformed_stats$Reference = "Genes"

6. Zbieraj statystyki przypisania rRNA.
    

   # collect the rrna assignment statistics
rrna_count_stats <- rrna_feature_counts$stat
# set the sample names as column names
colnames(rrna_count_stats) <- gsub(".bam", "", colnames(rrna_count_stats))

# only the number of assigned reads are important in this case
rrna_counts = as.data.frame(t(rrna_count_stats[rrna_count_stats$Status == "Assigned",2:ncol(rrna_count_stats)]))
colnames(rrna_counts) = "Alignments"

# setup names and categories for plotting
rrna_counts$Sample = rownames(rrna_counts)
rrna_counts$Group = "rRNA"
rrna_counts$Reference = "rRNA"

7. Statystyki odczytu odwzorowania jako wykres słupkowy.
    

   # combine the proteind coding and rRNA gene statistics
mapping_stats = rbind(transformed_stats, rrna_counts)

# plot the assignment statistics
mapping_stats_plot = ggplot(data = mapping_stats, aes(x = Sample, y = Alignments)) +
 geom_col(aes(fill = Group), width = 0.7) +
 theme_bw() + facet_wrap(~Reference) +
 ylab("Number of Alignments") + xlab("Samples") +
 ggtitle("Read Mapping Numbers") +
 theme(text = element_text(size = 18)) +
 theme(axis.text.x = element_text(angle = 90, hjust = 0))
print(mapping_stats_plot)

8. Klasyfikuj geny w grupy według liczby przypisanych im odczytów i nakreśl wyniki jako wykres słupkowy.
    

   # get the reads per gene counts
gene_count_matrix = gene_feature_counts$counts
# set the sample names as column names
colnames(gene_count_matrix) = gsub(".bam", "", colnames(gene_count_matrix))

# group and count genes in classes of specific read counts
read_count_classes = data.frame("no_reads"=colSums(gene_count_matrix == 0),
     "at_least_1_read"=colSums(gene_count_matrix >= 1 & gene_count_matrix < 10),
     "at_least_10_reads"=colSums(gene_count_matrix >= 10 & gene_count_matrix < 100),
     "at_least_100_reads"=colSums(gene_count_matrix >= 100 & gene_count_matrix < 1000),
     "at_least_1000_reads"=colSums(gene_count_matrix >= 1000))

# setup data.frame for plotting
read_count_classes$Sample = rownames(read_count_classes)
melt_rcc = melt(read_count_classes)
melt_rcc$variable = as.character(melt_rcc$variable)

# sort the gene groups
melt_rcc$variable = factor(melt_rcc$variable, levels=c("no_reads",
         "at_least_1_read",
         "at_least_10_reads",
         "at_least_100_reads",
         "at_least_1000_reads"))

# plot the data as bar plot
gene_coverage_plot <- ggplot(melt_rcc, aes(x=Sample, y=value, fill=variable)) +
geom_bar(stat="identity") +
ylab("Number of Genes") + xlab("Samples") +
ggtitle("Number of Reads per Gene") +
guides(fill=guide_legend(title="Number of assigned reads")) +
theme(text = element_text(size = 18)) +
theme(axis.text.x = element_text(angle = 90, hjust = 0))
print(gene_coverage_plot)

Rysunek 3: Przeczytaj przegląd mapowania i statystyki. Większość próbek wykazuje dużą liczbę przypisanych odczytów (lewa strona, brązowa) i niską liczbę odczytów rRNA (prawa strona, niebieska). Próbka s2_r3 ma wyjątkowo dużą liczbę odczytów wielomapowanych (lewa strona, zielona), odpowiadających wysokiej liczbie odczytów rRNA (prawa strona). s2_r4 próbki wykazuje dużą liczbę niemapowanych odczytów w połączeniu z oczekiwaną liczbą odczytów rRNA, co sugeruje skażenie odczytami od innego organizmu. Proszę kliknąć tutaj, aby zobaczyć większą wersję tej figurki.

4. Odtworzenie jakości.

1. Upewnij się, że replikacje z tej samej próbki lub tkanki wykazują wyższą korelację między sobą niż replikacje z różnych próbek lub tkanek. Wykonaj zgrupowaną mapę ciepła korelacji replikacyjnych, aby zwizualizować dane, używając surowych liczb odczytów lub normalizowanych liczb transkrypcji na milion (TPM)¹⁵ .
    

   # load the pheatmap library
library(pheatmap)

# compute TPM values
Length_kb = gene_feature_counts$annotation$Length / 1000 # get gene lengths in kb
RPK = gene_feature_counts$counts / Length_kb # normalize read counts by gene length
scaling_factors = colSums(RPK) / 1e6 # get scaling factor based on the sum normalized read counts
TPM = RPK / scaling_factors # scale the normalized read counts to 1e6

# plot the sample/replicate correlation heatmap
# the pearson correlation is computed on the log2 transformed TPM values
pheatmap(cor(log2(TPM+1)), fontsize = 18, main="Sample Correlation Heatmap")

Rysunek 4: Klasyfikacja liczby odczytów genów. Wszystkie geny w próbce klasyfikuje się w jedną z pięciu kategorii na podstawie liczby przypisanych odczytów. Na czerwono pokazano liczbę genów bez przypisanych odczytów. Próbki z niską liczbą przypisanych odczytów (s2_r1 do s1_r4) mają większą liczbę genów z 10 do 100 przypisanych odczytów oraz mniejszą liczbę genów z ponad 1000 odczytów. Geny o niskiej ekspresji mogą zostać przeoczone w tych próbkach. Proszę kliknąć tutaj, aby zobaczyć większą wersję tej figurki.

Wyniki

Potok jest w pełni zaimplementowany jako skrypt R i był testowany na systemach Linux i Mac OS. Użytkownicy Windows mogą korzystać z Windows Subsystem for Linux (WSL). Kod i dane testowe są dostępne jako repozytorium GitHub: https://github.com/oliverrupp/rup. Dane sekwencjonowania są dostępne w ramach projektu ENA EBI PRJEB96400.

Dziesięć próbek zostało sztucznie utworzonych z dwóch rzeczywistych próbek, aby zilustrować różnorodne problemy, które mogą wystąpić podczas kontroli jakości masowego RNA-seq z wykorzystaniem Rup do analizy. Próbka s2_r1 została zaprojektowana tak, aby wykazywać niską całkowitą liczbę odczytów, próbka s2_r2 zawiera dużą część odczytów niskiej jakości do odrzucenia w procesie przycinania. s2_r3 próbki obejmuje dużą część odczytów rRNA, a s2_r4 próbki zawierają odczyty zanieczyszczeń, które nie dopasowały się do genomu referencyjnego. Nazwy próbek s1_r5 i s2_r5 zostały zamienione, aby zilustrować niską korelację replikacji.

Sekcja protokołu 2 umożliwia identyfikację próbek o niskich liczbach odczytów przed lub po przycięciu (Rysunek 2). Wykres słupkowy pokazuje niższą wartość odczytu w próbkowej s1_r1 zarówno przed, jak i po przycięciu. Tutaj początkowa liczba odczytów była niska. Niska liczba odczytów próbek s1_r2 po przycięciu sugeruje dużą liczbę sekwencji adapterów, błędów sekwencjonowania, sekwencji primerów, rozciągnięć poli-A/T, które zostały usunięte podczas procesu przycinania. Degradacja wejściowego RNA może również zmniejszyć liczbę wysokiej jakości odczytów.

Sekcja protokołu 3 identyfikuje problemy w przypisaniach odczytu. Rysunek 3 przedstawia podwyższoną liczbę odczytów wielomapowanych w próbie s2_r3 (zielono), a także wysoką liczbę odczytów rRNA. Zanieczyszczenie s2_r4 próbki widoczne jest po dużej części odczytów nieprzypisanych do genomu referencyjnego (różowy). Odczyty te nie korelują z rRNA, lecz z sekwencjami organizmu niebędącego celem. Rysunek 4 przedstawia ogólne problemy związane z niską liczbą przypisanych odczytów w transkryptomach. W próbkach o niskim wskaźniku odczytów przypisanych unikalnie do genomu referencyjnego (s2_r1 do r4) tylko około jednej trzeciej genów przypisano więcej niż 100 odczytów, podczas gdy w pozostałych próbkach około połowa genów należy do tych klas. Odczyty genów o bardzo niskiej ekspresji mogą nie być znalezione w próbkach s2_r1 r4, dlatego analiza porównawcza ekspresji z użyciem tych próbek będzie bardzo niewiarygodna i należy jej unikać.

Sekcja protokołu 4 może być użyta do identyfikacji replikowanych wartości odstających. Replikacje tej samej próbki/warunku/tkanki powinny wykazywać wyższą korelację między sobą niż replikacje z innych próbek/warunków/tkanek. Rysunek 5 przedstawia korelację heatmapę dwóch próbek (S1 i S2) z pięcioma replikacjami każda. Dendrogramy na górze i z boku wykresu pokazują dwa skupiska z pięcioma replikami w każdej gromadzie. Lewy klaster zawiera cztery repliki próbki 1 i jedną replikę próbki 2 (s2_r5), prawy klaster zawiera cztery repliki próbki 2 oraz jedną replikę próbki 1 (s1_r5). W takim przypadku, gdy nazwy próbek s1_r5 i s2_r5 są zamieniane z powrotem, każdy klaster zawiera wszystkie replikacje jednej próbki, co może wskazywać na błąd w etykietowaniu replikacji. Inne powody, dla których replikacje nie grupują się razem, mogą być brak rozróżnienia między próbkami/stanami/tkankami lub skupianie replikacji wyłącznie ze względów jakościowych. To drugie może wystąpić, gdy wszystkie replikacje o wyjątkowo niskiej liczbie odczytów lub te o wyjątkowo wysokiej liczbie odczytów po przycięciu i mapowaniu tworzą klaster.

Rysunek 5: Przykładowa mapa ciepła korelacji. Mapa korelacji próbek opiera się na logariatucie₂ przekształconych wartości TPM i pokazuje dwa odrębne klastry po pięć próbek każda. Oczekuje się, że repliki biologiczne/techniczne będą wykazywać wyższą korelację między sobą niż te z innych tkanek/terapii. Lewy klaster zawiera cztery repliki próbki 1 i jedną replikę próbki 2 (s2_r5), prawy klaster zawiera cztery repliki próbki 2 oraz jedną replikę próbki 1 (s1_r5). Pojedyncze repliki należy sprawdzić pod kątem możliwej wymiany próbek lub efektów wsadowych, aby wyjaśnić ich klasteryzację na mapie ciepła. Proszę kliknąć tutaj, aby zobaczyć większą wersję tej figurki.

Tabela 1: Porównanie pipeline'ów oceny jakości RNA-seq. Szczegółową analizę można znaleźć w Pliku Uzupełniającym 1. Prosimy kliknąć tutaj, aby pobrać tę tabelę.

Plik uzupełniający 1: Wybór parametrów i odniesień. Analizy pokazują wpływ wyboru parametrów i odniesień na ogólne wyniki QC. Proszę kliknąć tutaj, aby pobrać ten plik.

Dyskusja

Jakość analiz różnicowej ekspresji genów zależy w dużej mierze od dwóch czynników: liczby odczytów zsekwencjonowanych w każdej próbce i replikacji¹⁶ oraz liczby replikacji na próbkę^17,18. Tutaj przedstawiamy przyjazny dla użytkownika pipeline Rup, aby określić liczbę odczytów odpowiednich do kwantyfikacji ekspresji genów w każdym replikacie. Różne metryki pozwalają badaczom zrozumieć, dlaczego repliki wykazują niskie przypisane liczby odczytów oraz identyfikować problemy w korelacji próbek. Chociaż Rup został opracowany do oceny jakości próbek RNA-seq roślin, jest równie odpowiedni dla innych organizmów eukariotycznych, nie wymagając dalszych korekt dla próbek nieroślinnych (Plik uzupełniający 1).

Pierwszy krok Rup określa całkowitą liczbę odczytów przed i po przycięciu odczytów. Liczba zsekwencjonowanych odczytów decyduje o liczbie wykrywalnych genów, a jeśli liczba odczytów jest zbyt niska, wiele genów różnicowo wyrażanych pozostanie niewykrytych. Znaczna część odczytów odrzuconych podczas procesu przycinania jakości i filtrowania może wskazywać na degradację RNA, co może prowadzić do niedoszacowania ekspresji genów o wysokiej degradacji. Jednak to, czy próbka może być użyta do zastosowań w dalszych etapach, zależy od organizmu docelowego i celu badania. Na przykład, aby uchwycić ekspresję większości genów roślinnych, z naszego doświadczenia potrzebne jest 30 do 50 milionów odczytów, podczas gdy dla grzybów wystarcza tylko 10 milionów. W związku z tym Rup nie będzie definiował progów jakości wykluczających problematyczne próby, lecz dostarczał miar do identyfikacji różnorodnych problemów, które mogą się pojawić.

Drugi etap (jakość mapowania) ocenia dokładność i wiarygodność wyrównania odczytu do genomu referencyjnego, określając ich przydatność do dalszych analiz. Nie wszystkie odczyty sekwencjonowane mogą być używane do obliczeń obfitości genów; odczyty, które nie są zgodne z genomem lub transkryptomem, albo odwzorowują się na wiele miejsc w genomie, są w^{większości przypadków} ignorowane i nie przyczyniają się do ogólnej liczby odczytów. Wyniki drugiego etapu potoku można wykorzystać do zrozumienia, dlaczego odczyty nie są używane w obliczeniach obfitości. Duża liczba niezmapowanych odczytów może wskazywać na skażenie podczas ekstrakcji RNA (np. przez patogeny roślinne lub roślinożerne) lub niepełny genom referencyjny. Niekompletne modele genów genomu referencyjnego mogą skutkować dużą liczbą odczytów "brak cech". Znaczna część wielokrotnie mapowanych odczytów może być spowodowana dużą liczbą genów rRNA w bibliotece, co wskazuje na niewystarczające usunięcie rRNA podczas procesu przygotowania biblioteki sekwencjonowania (jeśli wielokrotnie mapowane odczyty rzeczywiście pochodzą z rRNA, można je dalej analizować, dostarczając plik adnotacji rRNA do potoku, który automatycznie oblicza liczbę możliwych odczytów rRNA).

Trzecim, równie ważnym wskaźnikiem jakości jest korelacja między replikacjami tej samej próby. Zazwyczaj korelacja między replikacjami tej samej próby powinna być wyższa niż korelacja między replikacjami różnych próbek. Niska korelacja między replikacjami może wskazywać na dużą zmienność biologiczną między replikacjami, wysokie podobieństwo między próbkami/warunkami/tkankami, inne efekty partii, a nawet zamianę próbek lub błędne oznakowanie. Rup oblicza korelacje parowe między wszystkimi próbkami i tworzy grupową mapę cieplną korelacji próbek. Dodatkowo na próbkach można wykonać analizę głównych składników (PCA), aby zidentyfikować możliwe efekty wsadowe, które należy uwzględnić w dalszych analizach. Chociaż efekty partii można skorygować w analizach dalszych, może być lepiej usunąć problematyczne próbki, jeśli zmierzona różnica jest zbyt duża.

Materiał wejściowy ma bezpośredni wpływ na jakość sekwencjonowania: pobieranie próbek tkanek, ekstrakcja RNA oraz przygotowanie biblioteki to kluczowe kroki minimalizujące problemy z jakością RNA-seq. Utrzymuje się stabilne warunki, jeśli to możliwe, na przykład poprzez stosowanie komór wzrostowych. Środki zwalczania szkodników powinny być przeprowadzane terminowo, ponieważ do pobierania próbek należy wybierać wyłącznie zdrowe osoby. Zaleca się również pobieranie próbek tego samego dnia i/lub o tej samej porze, aby zmniejszyć zmienność transkryptomu okołodobowego. Dostępnych jest wiele zestawów do ekstrakcji RNA, a wybór odpowiedniego zestawu dla gatunku docelowego może poprawić jakość mRNA. Protokoły przygotowania biblioteki powinny zawierać kroki wzbogacania RNA poliA+ w celu zminimalizowania frakcji rRNA.

Rup może być stosowany jako początkowy krok kontroli jakości w różnicowej analizie ekspresji genów. Ta kontrola jakości jest wymagana z kilku kluczowych powodów, ponieważ zapewnia jakość danych wejściowych (identyfikuje repliki/próbki niskiej jakości, może ustalać progi jakości dla głębi odczytu i szybkości mapowania) oraz identyfikuje problemy techniczne, takie jak błędy sekwencjonowania, efekty wsadowe czy błędne oznakowanie. Jeśli wejście RNA jest odpowiedniej jakości, Rup może pomóc poprzez przycinanie odczytów niskiej jakości lub identyfikację błędnie oznakowanych próbek. Jednak Rup nie jest w stanie zrekompensować niskiej jakości wejściowego RNA. W przypadku niskiej jakości danych RNA lub błędnych sekwencjonowania może być konieczne ponowne pobranie próbek, dostosowanie protokołu ekstrakcji RNA i/lub powtórzenie sekwencjonowania. To samo dotyczy próbek z dużą częścią niezmapowanych odczytów, które mogły być spowodowane skażeniem. Jednak z naszego doświadczenia wynika, że nie wszystkie ekstrakcje RNA można po prostu powtórzyć ze względu na brak dostępności materiału. W tym przypadku nie nadają się do niektórych dalszych zastosowań: próbki z niską liczbą odczytów pojedynczych mapowanych nie powinny być analizowane w analizach różnicowej ekspresji genów, ale nadal zawierają informacje do analizy obecności transkrypcji. W takim przypadku brak transkrypcji w tkankach/leczeniu/warunkach nie może być brany pod uwagę przy analizie i ilościowym określeniu obfitości transkryptów.

Rup wiąże się z pewnymi ograniczeniami. Na przykład nie testuje bezpośrednio degradacji RNA, ponieważ przed przygotowaniem i sekwencjonowaniem biblioteki wymagane są bezpośrednie pomiary integralności RNA. Jednak w repozytorium tego potoku znajduje się skrypt do identyfikacji degradacji RNA z wykorzystaniem modułów RSeQC geneBodyCoverage.py i tin.py. Ponadto Rup nie testuje treści GC ani długości transkrypcji. Rozmiar genomu referencyjnego jest obecnie ograniczony do 4 Gb, a według naszego doświadczenia moduł mapowania odczytów zawyża odczyty wielokrotnie mapowane w poliploidach, co ilustruje Plik Uzupełniający w porównaniu mapowania odczytu na haploida z wersjami diploidalnymi E. californica, dlatego preferowanym wejściem dla tego pipeline'u jest genom haploidalny. Jakość wyjścia Rup w dużej mierze zależy od jakości genomu referencyjnego i adnotacji, tak że niepełna lub silnie fragmentowana sekwencja genomu może prowadzić do przeszacowania liczby odczytów niemapowanych. Co więcej, niekompletne adnotacje modeli genu prowadzą do przeszacowania liczby nieprzypisanych odczytów. Kompletność i duplikację sekwencji genomu oraz adnotacji genu można wywnioskować za pomocą narzędzi takich jak BUSCO²⁰.

Rup jest samodzielnym narzędziem dla wszystkich początkowych kroków kontroli jakości, niezbędnych w eksperymentach RNA-seq. W przeciwieństwie do RSeQC i RNA-SeQC, wstępne przetwarzanie surowych odczytów sekwencjonowania nie jest wymagane. Analiza jakości sekwencjonowania nie może być przeprowadzona za pomocą RNA-SeQC, a heatmapy korelacji próbek nie są obliczane przez RSeQC i RNA-SeQC (Tabela 1). Ponadto wyjście normalizacyjne jest w FPKM w RSeQC, a wizualizacja wyników nie jest standardem wszystkich modułów w RSeQC i RNA-SeQC. W związku z tym w dwóch alternatywnych narzędziach nie testuje się kilku problemów kontroli jakości, takich jak niska liczba odczytów, wysoki udział przyciętych odczytów oraz identyfikacja replikowanych wartości odstawających. Porównanie Rup, RSeQC i RNA-SeQC jest dostępne w Pliku Uzupełniającym 1. Ponadto Rup może być używany jako uzupełnienie RNA-SeQC lub RSeQC, na przykład posortowane pliki BAM wytworzone przez ten pipeline mogą służyć jako dane wejściowe dla tych narzędzi.

Obecnie Rup jest zoptymalizowany pod niewielką liczbę próbek, ale najbardziej czasochłonne kroki, takie jak przycinanie jakości i mapowanie odczytu, można obliczyć precom, na przykład w klastrze komputerowym lub infrastrukturze chmurowej. Dla genomów większych niż 4 GB jest to obowiązkowe, ponieważ nakładnik Rsubread jest ograniczony do genomów mniejszych niż 4 Gb. Adnotacja rRNA wykonywana jest za pomocą barrnap, który identyfikuje wysoce zachowane geny rRNA. Według naszego doświadczenia, adnotacja genu rRNA nie wymaga kompleksowości, ponieważ nawet jeśli brakuje niektórych nieprawidłowych genów rRNA, ogólny przegląd obecności rRNA w zbiorze danych jest wystarczający do oceny jakości zbioru. Przyszłe wersje Rup mogą przejść na inną metodę mapowania, jak narzędzie pseudo-wyrównania salmon²¹, aby skrócić czas działania. Ponadto do Rup można dodać więcej metod normalizacji i korekty, takich jak błędy GC lub długości.

Podsumowując, Rup dostarcza niezbędnych informacji, aby zapewnić wiarygodność i powtarzalność analizy danych RNA-seq. Ten pipeline kompleksowo raportuje główne metryki jakości RNA-seq i generuje pliki wyjściowe do bezpośredniego wykorzystania w dalszych analizach. Został zaprojektowany jako samodzielne narzędzie dla badaczy z minimalną wiedzą bioinformatyczną do oceny jakości danych z zakresu sekwencjonowania podstawowego za pomocą intuicyjnej wizualizacji.

Materiały

Lista materiałów użytych w tym artykule

Nazwa	Firma	Numer katalogowy	Komentarze
fastqcr	R	0.1.3
FUJITSU ESPRIMO D958 Desktop	FUJITSU		pipeline został opracowany i przetestowany na Ubuntu Linux 24.02 (32 Gb RAM)
getopt	R	1.20.4
ggplot2	R	3.5.2
MacBook Pro	Apple		pipeline był testowany na maxOS 15.4.1 (16 GB RAM)
mapa pheatmap	R	1.0.13
R		4.4.3
Reshape2	R	1.4.4
rfastp	Bioconductor	1.16.0
rsamtools	Bioconductor	2.22.0
rsubread	Bioconductor	2.20.0