Rup (RNA-seq Kullanılabilirlik Değerlendirme Hattı) - Ökaryotlarda Toplu RNA-seq Deneyleri için Kalite Kontrol

Özet

Bu protokol, sınırlı biyoinformatik deneyimine sahip ıslak laboratuvar biyologları için RNA-seq deneyleri için ilk kalite kontrolüne izin verir.

Özet

Moleküler bitki bilimindeki modern yaklaşımlar, örneğin, tedaviler üzerine transkriptomlardaki küresel değişiklikleri izlemek veya düzenleyici yolların temel bileşenlerini belirlemek için genellikle toplu RNA-seq deneyleri gerektirir. Sonuç olarak, bitki bilimlerinin çeşitli alanları, bilimsel ilerleme için yüksek kaliteli ve tekrarlanabilir toplu RNA-seq verilerine güvenmektedir. Bununla birlikte, deneyimlerimize göre, RNA-seq veri kümelerinde kalite kontrol önlemlerinin uygulanması ve uygulanması hakkında bilgi ve uygulama genellikle eksiktir. Burada, Rup'u tanıtıyoruz (RNA-seq kullanılabilirlik değerlendirme boru hattı) toplu RNA-seq verilerinin kalite kontrolü için, sonraki gen ekspresyon analizleri için, R. Rup, yüksek kaliteli, aşağı akış gen ekspresyon deneyleri için uygun olan ve genel ileri analizler için uygun olmayan yüksek kaliteli dizileme verileri arasında ayrım yapmaya yardımcı olur. Rup, yetersiz okuma sayıları veya haritalama, kontaminasyonların tanımlanması, toplam RNA-seq verilerindeki rRNA fraksiyonlarının miktarının belirlenmesi, benzerlik testinin çoğaltılması, gösterim için gerçek verilerin kullanılması ve sezgisel görselleştirme sunulması gibi yaygın olarak karşılaşılan birkaç soruna yönelik testler içerir. Rup, standartlaştırılmış transkriptom analizinden önce deneysel eksiklikleri belirlemek için bir dizi araç sağlar ve böylece bireysel araştırmacılar ve alan için veri kalitesini artırır. Bu, toplu RNA-seq veri analizine olan güveni artırır ve minimum kalite kontrol kriterlerini tanımlayan gelecekteki kılavuzlar için bir temel sağlar, böylece yayınlanan RNA-seq verilerinin güvenilirliğini ve şeffaflığını artırır. Rup ve test verileri https://github.com/oliverrupp/rup'da mevcuttur.

Giriş

Transkriptomik (RNA-seq) deneyleri, fenotipleri şekillendiren transkripsiyonel imzaları kapsamlı bir şekilde sorgular. Bu yaklaşım, örneğin gelişim, bitki patojen etkileşimi veya abiyotik stres direnci^{1, 2 ,3} ile ilgili tek veya birlikte eksprese edilen genleri ve biyolojik süreçleri tanımlamak için bitki moleküler genetiğinde yeri doldurulamaz hale geldi. RNA-seq teknolojisindeki son gelişmeler, özgüllüğü artırmış ve farklı izoformların ve varyantların tek baz çözünürlüğünde saptanmasını sağlayarak, daha büyük indellerden tek nükleotid polimorfizmlerine (SNP'ler) kadar dizi varyasyonlarının tanımlanmasına izin vermiştir. RNA-seq ile elde edilen veriler, hem yüksek oranda bol hem de düşük eksprese edilmiş transkriptlerin saptanmasını kolaylaştıran geniş bir dinamik aralık ile karakterize edilir ve tutarlılık için uygun deneysel koşullar gerektirir. Ek olarak, RNA-seq veri kümeleri büyüktür ve bu nedenle analiz edilmesi ve depolanması hesaplama açısından yoğundur, bu da verimli veri yönetimi ve kapsamlı bilgi işlem kaynakları gerektirir. RNA-seq, iş akışının her adımında yüksek kaliteli girdi gerektirir, çünkü herhangi bir aşamadaki zayıflıklar yayılabilir ve sonuçları tehlikeye atabilir. Kütüphane hazırlığı sırasında zayıf RNA bütünlüğü veya teknik sorunlar, önyargıya ve doğruluğun azalmasına yol açacaktır. Yüksek kaliteli ham okumalar elde etmek için parametreler, yeni nesil dizileme (NGS) tesisleri arasında farklılık gösterebilir.

Deneyimlerimize göre, RNA-seq kitaplığı hazırlığı için girdi olarak yalnızca RNA bütünlük numarası (RIN) 7'nin üzerinde olan RNA'yı kullanmanızı öneririz, bu da büyük ölçüde bozulmamış mRNA yapısının göstergesidir. Başarılı dizileme için, standart kitaplık hazırlama protokolleri için 50-200 ng/μL konsantrasyonda yaklaşık 2 μg toplam RNA gereklidir. RNA saflığı, bir spektrofotometre kullanılarak 1.8 ile_{2.1 arasında bir OD 260/280} oranı ve 1.5'ten büyük bir OD_260/230 oranı ile doğrulanmalıdır. Yetersiz dizileme derinliği, düşük haritalama oranları veya referans genoma yanlış hizalama, gen ekspresyonunu ve ekleme nicelemesini daha da bozabilir. Ayrıca, farklı doku veya tedavi örnekleri arasında ayrım yapılmaması ve kopyalar arasındaki yüksek değişkenlik gibi yetersiz deneysel tasarım, gürültüye neden olabilir ve tekrarlanabilirliği azaltabilir. Birçok laboratuvar transkriptomları rutin olarak analiz ederken, ilk temel analiz adımlarının sıkı kalite kontrolleri genellikle rapor edilmez veya yayınlarda tamamen bulunmayabilir. Bu, transkriptom analizinden elde edilen sonuçların aşırı yorumlanmasına ve sonuç olarak tekrarlanamayan sonuçlara yol açabilir.

Burada, transkripsiyonel değişiklikleri ölçmek için mRNA profillerinin yüksek kaliteli transkriptom analizleri için gereken ilk adımların kalite kontrolü için bir iş akışı sağlıyoruz. Amacımız, biyoinformatik konusunda sınırlı bilgiye sahip ıslak laboratuvar biyologlarının birincil transkriptomik verilerini değerlendirmelerini sağlamaktır. Rup'a (Şekil 1), R'nin temel bilgilerine aşina olan araştırmacılar erişebilir. Burada sunulan iş akışının yürütülmesi, araştırmacılara, sonraki analizler için potansiyel sınırlamaları da dahil olmak üzere, birincil verilerinin ayrıntılı bir şekilde anlaşılmasını sağlayacaktır. Bildiğimiz kadarıyla, yüksek kaliteli verileri düşük kaliteli verilerden ayırt etmek için yönergelerle birlikte pratik bir birincil transkriptom veri değerlendirme hattı şu ana kadar eksiktir.

Şekil 1: RNA-seq'in silico kalite kontrol hattının iş akışı. Kalite kontrol için girdi dosyaları, bitki materyalinin dizilenmesiyle oluşturulan RNA-seq verilerinden ve ayrıca halka açık veri kümelerinden türetilir. Sağlanan R ardışık düzeni, dizi kalitesini üç ana yaklaşımla değerlendirir: sıralama kalitesi, haritalama kalitesi ve çoğaltma kalitesi. Çeşitli kalite ölçümleri hesaplanır ve istatistiksel sonuçlar, ggplot2 ve pheatmap (örneğin, çubuk grafikler ve ısı haritaları) gibi yaygın R paketleri kullanılarak görselleştirilir. Bu rakamın daha büyük bir sürümünü görüntülemek için lütfen buraya tıklayın.

Daha önce bildirilen değerlendirme işlem hatları, önceden işlenmiş veriler gerektiriyordu (örneğin, .bam dosyaları olarak okuma hizalaması), tüm metrikleri kapsamıyor veya artık korunmuyor ^4,5,6. Burada sunulan iş akışının avantajı, en yaygın kalite kontrol sorunlarını tek bir ardışık düzende birleştirerek kapsamlı olmasında yatmaktadır. Ayrıca, tüm aşağı akış analiz uygulamalarına uygun yüksek kaliteli veriler için örneklerin yanı sıra düşük kaliteli veriler için örnekler de sağlıyoruz ve daha fazla analiz için bunların belirli sınırlamalarını tartışıyoruz. Daha önce yayınlanmış birkaç rapor, RNA-seq analiz araçlarının amacını açıklar ve performanslarının karşılaştırmalı değerlendirmelerini sağlar ^7,8. Bununla birlikte, RNA-seq kalite kontrol ve metodoloji raporlama standartları standartlaştırılmamıştır ve transkriptomik deneylerin tekrarlanabilirliğini ve biyolojik olarak anlamlı yorumlarını ağırlaştırmaktadır.

Rup, standart yüksek kaliteli araçları, kalite kontrol analizini ve sonuçların görselleştirilmesini entegre eder. Rup, girdi olarak ham dizileme okumaları ve açıklamalı bir genom gerektirir ve Windows sistemlerinde Mac OS, Linux ve "Linux için Windows Alt Sistemi" (WSL) üzerinde çalışır. Dizilemenin kalitesi için araçları entegre eder ve bir genoma okuma eşlemesini ölçer, örneklerde rRNA içerik ölçümünü içerir ve kopya korelasyonunu değerlendirmek için örnek korelasyonu sunar. Test verileri, kütüphane hazırlığı için ultra düşük girdili bir protokol olan Eschscholzia californica bitki türünün iki farklı aşamasının merkezi meristem dokusunun lazer mikrodiseksiyonu kullanılarak elde edildi ve Novaseq 6000'de dizilendi.

Protokol

Rup, gereken minimum girdi dosyasıyla tek bir komut dosyası olarak çalıştırılabilir. Burada gösterildiği gibi boru hattı, Illumina çift uçlu RNA dizileme verilerini gerektirir. Aynı analiz adımlarıyla tek uçlu sıralama için ayarlanmış bir işlem hattı da GitHub deposunda depolanır. Yalnızca bir fasta dosyası olarak genom dizisi, gtf dosyaları olarak gen modeli ve rRNA açıklamaları ve fastq.gz dosyaları olarak ham dizileme okumaları gereklidir. Genom ve açıklama dosyalarının, reference_folder değişkeninde belirtilen bir klasörde genome.fa, annotation.gtf ve rRNA.gtf olarak sağlandığından emin olun. Tüm sıralama dosyaları .fq.gz dosyalar read_file_folder değişkenine kaydedildiğinde, Rup aşağıdaki gibi çalıştırılabilir:

1. Hazırlıklar:

1. Bioconductor kullanarak gerekli tüm R paketlerini kurun.

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

2. Gerekli tüm dosyaları ayarlayın.
   1. Gerekli tüm dosyaları içerecek bir kaynak klasör oluşturun. Referans genom dizisini kaynak klasöre "reference/genome.fa" adlı fasta dosyasına ekleyin. Gen modeli açıklamasını bir gtf dosyası olarak sağlayın: "reference/annotation.gtf", İsteğe bağlı olarak, rRNA gen açıklamasını gtf dosyası olarak sağlayın: "reference/rRNA.gtf".
   2. Tüm sıralama okumalarını fastq dosyaları olarak "reads" klasörüne ekleyin. Tüm fastq dosya adlarının aynı kalıbı izlediğinden emin olun; İleri ve geri okumalar için _1.fastq.gz ve _2.fastq.gz.

      # 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. Parametreyi bilgi işlem kaynaklarına göre ayarlayın.

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

4. Parametreyi sıralama yöntemine göre ayarlayın.
    

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

   
   
   NOT: Minimum okuma uzunluğu 25 bp'den küçük olmamalıdır, çünkü daha küçük okumalar eşlemenin bozulmasına neden olabilir. Daha büyük değerler, sıralama kalitesine bağlı olarak kırpma sırasında daha benzersiz eşlenmiş okumaların yanı sıra daha fazla atılan okumayla sonuçlanacaktır. Bu parametreler hakkında daha fazla ayrıntı eklerde açıklanmıştır.

2. Sıralama kalite değerlendirmesi

NOT: Bu adım, her bir numuneyi kırpmadan önce ve sonra okuma sayısını gösteren bir çubuk grafiği oluşturacaktır.

1. Yüksek verimli dizi verilerinin genel kalite kontrolü için bir araç olan fastqc^9'u çalıştırın ve toplam sıralı okuma sayısı, okuma uzunluğu, ortalama GC içeriği, dizileme adaptörü kontaminasyonu ve aşırı temsil edilen diziler (örneğin, rRNA okumaları) hakkında bilgi sağlar. Sıralama kalitesini baz ve dizi başına kalite raporlarına göre analiz edin. R paketi fastqcr^10'u kullanarak ham sıralama dosyalarında fastqc'yi çalıştırın ve tüm girdi dosyalarını aynı dizine (read_file_folder) yerleştirin. qc_aggregate işlevini kullanarak fastqc_folder'den gelen çıktı dosyalarını bir R nesnesine özetleyin.
   NOT: Sonuçlar fastqc_folder dizinine kaydedilecektir. Ek istatistikler için fastqc tarafından çıktı klasöründe oluşturulan .html dosyalarına erişilebilir.
    

   # 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. Vektör sample_prefixes depolanan tüm örnek önekleriyle, önekler üzerinde yineleme yapın ve tüm örneklerde kalite kırpma aracı fastp'yi çalıştırın. Adaptör dizilerini, primerleri, düşük kaliteli bazları, poli-A/T dizilerini çıkarmak için kırpın ve R paketi Rfastp^12'de bulunan Fastp¹¹ tarafından yürütülür. Kırpılan fastq dosyaları trimmed_read_folder klasöründe saklanacaktır.
    

   # 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. Kırpılan okumalarda fastqc komutunu yeniden çalıştırın.
    

   # 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. Tüm numuneler için kırpmadan önce ve sonra okuma sayısını toplayarak kırpmayı değerlendirin. ggplot2^13'ü kullanarak okunan sayılarla bir çubuk grafiği oluşturun.
    

   # 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)

Şekil 2: Sıralama ve kırpma sonuçları. Kırpmadan önce (kırmızı) ve sonra (yeşil) sayıları okuyun. Örnek s2_r1, kırpmadan önce zaten düşük bir okuma sayısına sahipken, kırpma, örnek s2_r2 okumalarının büyük bir bölümünü kaldırdı. Diğer tüm örnekler, kırpmadan kaynaklanan kabul edilebilir okuma kaybı gösterir. Bu rakamın daha büyük bir sürümünü görüntülemek için lütfen buraya tıklayın.

3. Haritalama kalitesi

NOT: Bu adım, tek eşlenmiş okumaları (örneğin, protein kodlayan genlerden gelen transkriptler) çoklu eşlenmiş okumalardan (örneğin, rRNA okumaları) ve eşlenmemiş okumalardan (örneğin, kontaminasyonlar) ayırt etmek için çubuk grafiklerini hesaplar.

1. Okumaları referans genoma eşlemek için Rsubread¹⁴ paketini kullanın. İlk olarak, referans genom fasta dosyasının (genome_fasta_file) bir dizinini oluşturun.
   NOT: Bu adım, her referans genom için bir kez gerçekleştirilir.
    

   # 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. Tüm örnekler üzerinde yineleyin ve Rsubread paketinin align() işlevini kullanarak okumaları referans genoma hizalayın. Bu işlev, çıktı klasöründe .bam dosyaları oluşturur.
    

   # 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. Rsubread paketindeki featureCounts() işlevini kullanarak her gen için okumaları sayın. Ek açıklama dosyasının (annotation_file) GTF biçiminde olduğundan emin olun. Sayım, yalnızca genomla tek bir eşleşme ile okunur.
    

   # 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. Numunelerdeki rRNA içeriğini değerlendirmek için rRNA genlerine eşleme okumalarını sayın. Çoklu eşlenmiş okumaların burada sayılmasına izin verin. Bir GTF dosyasında (rrna_file) rRNA gen lokuslarını belirtin.
    

   # 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. featureCounts tarafından oluşturulan okuma ataması istatistiklerini toplayın. İstatistikler, farklı kategorilerdeki (örneğin, Atanmış, Eşlenmemiş, Çoklu Eşlenmiş) her örnek için hizalama sayısını içerir.
    

   # 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. rRNA atamasının istatistiklerini toplayın.
    

   # 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. Okuma eşleme istatistiklerini çubuk grafiği olarak çizin.
    

   # 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. Genleri, kendilerine atanan okuma sayısına göre gruplara ayırın ve sonuçları bir çubuk grafiği olarak çizin.
    

   # 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)

Şekil 3: Haritalamaya genel bakışı ve istatistikleri okuyun. Çoğu örnek, çok sayıda atanmış okuma (sol taraf, kahverengi) ve düşük sayıda rRNA okuması (sağ taraf, mavi) gösterir. Örnek s2_r3, yüksek bir rRNA okuma sayısına (sağ taraf) karşılık gelen alışılmadık derecede yüksek miktarda çoklu eşlenmiş okumaya (sol taraf, yeşil) sahiptir. Örnek s2_r4, beklenen bir rRNA okuma sayısı ile kombinasyon halinde çok sayıda eşlenmemiş okuma gösterir ve bu, başka bir organizmadan gelen okumalarla kontaminasyonu düşündürür. Bu rakamın daha büyük bir sürümünü görüntülemek için lütfen buraya tıklayın.

4. Kaliteyi çoğaltın.

1. Aynı numune veya dokudan alınan kopyaların, farklı örneklerden veya dokulardan alınan kopyalara göre birbirleriyle daha yüksek bir korelasyon gösterdiğinden emin olun. Ham okuma sayılarını veya Milyon Başına Transkript (TPM)¹⁵ normalleştirilmiş sayılarını kullanarak verileri görselleştirmek için yinelenen korelasyonların kümelenmiş bir ısı haritasını çizin.
    

   # 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")

Şekil 4: Gen okuma sayısı sınıflandırması. Bir numunedeki tüm genler, atanan okuma sayısına göre beş kategoriden birinde sınıflandırılır. Kırmızı ile gösterilenler, herhangi bir okuma atanmamış genlerin sayısıdır. Düşük sayıda toplam atanmış okumaya (s2_r1 ila s1_r4) sahip numuneler, 10 ila 100 atanmış okumaya sahip daha fazla sayıda gene ve 1000'den fazla okumaya sahip daha az sayıda gene sahiptir. Bu örneklerde düşük ekspresyona sahip genler gözden kaçabilir. Bu rakamın daha büyük bir sürümünü görüntülemek için lütfen buraya tıklayın.

Sonuçlar

Ardışık düzen tamamen bir R betiği olarak uygulanmıştır ve Linux ve Mac OS işletim sistemlerinde test edilmiştir. Windows kullanıcıları Linux için Windows Alt Sistemini (WSL) kullanabilir. Kod ve test verileri bir GitHub deposu olarak kullanılabilir: https://github.com/oliverrupp/rup. Sıralama verileri ENA EBI proje PRJEB96400 altında mevcuttur.

Analiz için Rup kullanılarak toplu RNA dizisinin kalite kontrolü sırasında karşılaşılabilecek çeşitli sorunları örneklemek için iki gerçek örnekten yapay olarak on örnek oluşturuldu. Örnek s2_r1, düşük bir toplam okuma sayısı sergileyecek şekilde tasarlanmıştır, örnek s2_r2, kırpma işlemi tarafından atılacak düşük kaliteli okumaların büyük bir kısmını içerir. Örnek s2_r3, rRNA okumalarının büyük bir bölümünü içerir ve örnek s2_r4, referans genomla eşlenemeyen kirletici okumaları içerir. s1_r5 ve s2_r5 örneklerinin adları, düşük kopya korelasyonunu göstermek için değiştirildi.

Protokol bölüm 2, kırpmadan önce veya sonra düşük okuma sayılarına sahip numunelerin tanımlanmasını sağlar (Şekil 2). Çubuk grafiği, kırpmadan önce ve sonra örnek s1_r1 daha düşük okuma sayısını gösterir. Burada ilk okuma sayısı düşüktü. Kırpmadan sonra s1_r2 düşük okuma sayısı örneği, kırpma işlemi sırasında çıkarılan çok sayıda adaptör dizisi, sıralama hataları, primer dizileri, poli-A/T dizileri dizisini gösterir. Giriş RNA'sının bozulması, yüksek kaliteli okumaların sayısını da azaltabilir.

Protokol bölüm 3, okuma atamalarındaki sorunları tanımlar. Şekil 3 , numune s2_r3'daki (yeşil) yüksek sayıda çoklu eşlenmiş okumayı ve ayrıca yüksek sayıda rRNA okumasını göstermektedir. Numune s2_r4 kontaminasyonu, referans genoma (pembe) eşlenmeyen okumaların büyük bir kısmı ile belirgindir. Bu okumalar rRNA ile değil, hedef olmayan bir organizmadan gelen dizilerle ilişkilidir. Şekil 4 , transkriptomlarda düşük sayıda atanmış okuma ile ilişkili genel sorunları göstermektedir. Referans genoma (s2_r1 ila r4) benzersiz bir şekilde eşlenen düşük okuma oranına sahip örneklerde, genlerin sadece yaklaşık üçte biri 100'den fazla okumaya sahipken, diğer örneklerde genlerin yaklaşık yarısı bu sınıflara girer. Çok düşük ekspresyona sahip genlerin okumaları, r4'e s2_r1 numunelerde bulunmayabilir ve sonuç olarak, bu numunelerle karşılaştırmalı ekspresyon analizi oldukça güvenilmez olacaktır ve bundan kaçınılmalıdır.

Protokol bölüm 4, yinelenen aykırı değerleri tanımlamak için kullanılabilir. Aynı numunenin/koşulun/dokunun replikalarının, diğer numunelerden/koşullardan/dokulardan alınan replikalardan daha yüksek bir korelasyon göstermesi beklenir. Şekil 5 , her biri beş kopya içeren iki numunenin (S1 ve S2) bir korelasyon ısı haritasını göstermektedir. Grafiğin üstündeki ve yanındaki dendrogramlar, her kümede beş kopya bulunan iki kümeyi gösterir. Sol küme örnek 1'in dört yinelemesini ve örnek 2'nin (s2_r5) bir yinelemesini içerir, sağ küme örnek 2'nin dört yinelemesini ve örnek 1'in (s1_r5) bir yinelemesini içerir. Bu durumda, s1_r5 ve s2_r5 örnek adları yeniden değiştirildiğinde, her küme bir örneğin tüm yinelemelerini içerir ve bu da bir yineleme etiketleme hatasını gösterebilir. Kopyaların birlikte kümelenmemesinin diğer nedenleri, numuneler/koşullar/dokular arasında farklılaşma eksikliği veya kopyaların yalnızca sıralama kalitesi nedenleriyle kümelenmesi olabilir. İkincisi, son derece düşük olan tüm çoğaltmalar veya kırpma ve eşlemeden sonra son derece yüksek okuma sayılarına sahip olanlar bir küme oluşturduğunda ortaya çıkabilir.

Şekil 5: Örnek korelasyon ısı haritası. Örnek korelasyon ısı haritası, günlük₂ dönüştürülmüş TPM değerlerine dayanır ve her biri beş örnekten oluşan iki farklı kümeyi gösterir. Biyolojik/teknik kopyaların birbirleri arasında diğer dokulardan/tedavilerden daha yüksek bir korelasyon göstermesi beklenir. Sol küme örnek 1'in dört yinelemesini ve örnek 2'nin (s2_r5) bir yinelemesini içerir, sağ küme örnek 2'nin dört yinelemesini ve örnek 1'in (s1_r5) bir yinelemesini içerir. Tek kopyalar, ısı haritasındaki kümelenmelerini açıklamak için olası örnek değiştirme veya toplu etkiler açısından kontrol edilmelidir. Bu rakamın daha büyük bir sürümünü görüntülemek için lütfen buraya tıklayın.

Tablo 1: RNA-seq kalite değerlendirme boru hatlarının karşılaştırılması. Ayrıntılı bir analiz için Ek Dosya 1'e bakın. Bu Tabloyu indirmek için lütfen buraya tıklayın.

Ek Dosya 1: Parametre ve referans seçimi. Analizler, parametre ve referans seçiminin genel kalite kontrol sonuçları üzerindeki etkisini göstermektedir. Bu Dosyayı indirmek için lütfen buraya tıklayın.

Tartışma

Diferansiyel gen ekspresyonu analizlerinin kalitesi büyük ölçüde iki faktöre bağlıdır: her numunede sıralanan okuma sayısı ve replikasyon¹⁶ ve numune başına replika sayısı^17,18. Burada, her kopyada gen ekspresyonu ölçümü için uygun okuma sayısını belirlemek için kullanıcı dostu Rup boru hattını sunuyoruz. Farklı ölçümler, araştırmacıların kopyaların neden düşük atanmış okuma sayıları gösterdiğini anlamalarına ve örnek korelasyonundaki sorunları belirlemelerine olanak tanır. Rup, bitki RNA-seq örneklerinin kalite değerlendirmesi için geliştirilmiş olsa da, diğer ökaryotik organizmalar için eşit derecede uygundur ve bitki dışı örnekler için daha fazla ayarlama gerektirmez (Ek Dosya 1).

Rup'un ilk adımı, okuma kırpmadan önce ve sonra toplam okuma sayısını belirler. Sıralı okumaların sayısı, saptanabilir genlerin sayısını belirler ve okuma sayısı çok düşükse, diferansiyel olarak eksprese edilen birçok gen tespit edilmeden kalacaktır. Kalite kırpma ve filtreleme işlemi sırasında atılan okumaların büyük bir kısmı, RNA bozulmasını gösterebilir ve bu da yüksek oranda bozulmuş genler için ekspresyonun hafife alınmasına yol açabilir. Bununla birlikte, bir numunenin aşağı akış uygulamaları için kullanılıp kullanılamayacağı, hedef organizmaya ve çalışmanın amacına bağlıdır. Örneğin, çoğu bitki geninin ekspresyonunu yakalamak için, deneyimlerimize göre, 30 milyon ila 50 milyon okuma gerekirken, mantarlar için sadece 10 milyon yeterli olabilir. Bu nedenle Rup, sorunlu örneklerin hariç tutulması için kalite eşikleri tanımlamayacak, bunun yerine karşılaşılabilecek çeşitli sorunları belirlemek için ölçümler sağlayacaktır.

İkinci adım (haritalama kalitesi), referans genoma okuma hizalamasının doğruluğunu ve güvenilirliğini değerlendirir ve bunların aşağı akış analizleri için uygunluğunu belirtir. Sıralı okumaların tümü gen bolluğu hesaplaması için kullanılamaz; genom veya transkriptomla hizalanamayan okumalar veya genom üzerinde birden çok konumla eşlenen okumalar çoğu durumda göz ardı edilir¹⁹ ve genel okuma sayılarına katkıda bulunmaz. Boru hattının ikinci adımının sonuçları, bolluk hesaplamasında okumaların neden kullanılmadığını anlamak için kullanılabilir. Çok sayıda eşlenmemiş okuma, RNA ekstraksiyonu sırasında kontaminasyonu (örneğin, bitki patojenleri veya otçullarla) veya eksik bir referans genomunu gösterebilir. Referans genomun eksik gen modelleri, çok sayıda "özellik yok" okumasına neden olabilir. Çoklu eşlenmiş okumaların büyük bir kısmı, dizileme kitaplığı hazırlama işlemi sırasında yetersiz rRNA çıkarılmasını gösteren kitaplıktaki çok sayıda rRNA geninden kaynaklanabilir (çoklu eşlenmiş okumaların gerçekten rRNA'lardan türetilmesi durumunda, boru hattına bir rRNA açıklama dosyası sağlanarak daha fazla analiz edilebilirler, bu daha sonra olası rRNA okumalarının sayısını otomatik olarak hesaplar).

Üçüncü, eşit derecede önemli kalite ölçütü, aynı numunenin kopyaları arasındaki korelasyondur. Genel olarak, aynı numunenin kopyaları arasındaki korelasyon, farklı örneklerin kopyaları arasındaki korelasyondan daha yüksek olmalıdır. Kopyalar arasındaki düşük bir korelasyon, kopyalar arasında büyük biyolojik değişkenliği, numuneler/koşullar/dokular arasında yüksek benzerliği, diğer parti etkilerini ve hatta numune değişimini veya yanlış etiketlemeyi gösterebilir. Rup, tüm örnekler arasındaki ikili korelasyonları hesaplar ve örnek korelasyonlarının kümelenmiş bir ısı haritasını üretir. Ek olarak, daha sonraki aşağı akış analizlerinde kabul edilmesi gereken olası parti etkilerini belirlemek için numuneler üzerinde bir temel bileşen analizi (PCA) hesaplanabilir. Parti etkileri aşağı akış analizlerinde düzeltilebilse de, ölçülen fark çok yüksekse sorunlu numuneleri çıkarmak daha iyi olabilir.

Girdi materyalinin dizileme kalitesi üzerinde doğrudan etkisi vardır: doku örneklemesi, RNA ekstraksiyonu ve kitaplık hazırlığı, RNA-seq kalite sorunlarını en aza indirmek için kritik adımlardır. Mümkünse, örneğin büyüme odaları kullanılarak tutarlı koşullar korunmalıdır. Numune alma için yalnızca sağlıklı bireylerin seçilmesi gerektiğinden, haşere kontrol önlemleri zamanında alınmalıdır. Sirkadiyen transkriptom varyasyonunu azaltmak için numunelerin aynı gün ve/veya aynı saatte toplanması da tavsiye edilir. Çok sayıda RNA ekstraksiyon kiti mevcuttur ve hedef türler için uygun bir kit seçmek mRNA kalitesini artırabilir. Kitaplık hazırlama protokolleri, rRNA fraksiyonunu en aza indirmek için poliA+ RNA zenginleştirme adımlarını içermelidir.

Rup, diferansiyel gen ekspresyon analizinde ilk kalite kontrol adımı olarak kullanılabilir. Bu kalite kontrolü, girdi verilerinin kalitesini sağladığı (düşük kaliteli kopyaları/örnekleri tanımlar, okuma derinliği ve eşleme oranları için kalite eşiklerini ayarlayabilir) ve sıralama hataları, toplu etkiler veya yanlış etiketleme gibi teknik sorunları tanımladığı için birkaç kritik nedenden dolayı gereklidir. RNA girişi yeterli kalitedeyse, Rup düşük kaliteli okumaları kırparak veya yanlış etiketlenmiş örnekleri tanımlayarak yardımcı olabilir. Ancak Rup, zayıf giriş RNA kalitesini telafi edemez. Düşük kaliteli RNA veya hatalı dizileme verileri durumunda, numunelerin yeniden toplanması, RNA ekstraksiyon protokolünün ayarlanması ve/veya dizilemenin tekrarlanması gerekebilir. Aynısı, kontaminasyondan kaynaklanmış olabilecek büyük bir eşlenmemiş okuma oranına sahip numuneler için de geçerlidir. Bununla birlikte, deneyimlerimize göre, malzeme mevcudiyeti eksikliği nedeniyle tüm RNA ekstraksiyonları basitçe tekrarlanamaz. Bu durumda, bazı aşağı akış uygulamaları için uygun değildirler: düşük sayıda tek eşlenmiş okumaya sahip numuneler, diferansiyel gen ekspresyonu analizlerinde analiz edilmemelidir, ancak yine de transkript varlığının analizi için bilgi tutarlar. Bu durumda, transkript bolluğunun analizi ve miktarının belirlenmesi için dokularda/tedavilerde/koşullarda transkript bulunmaması dikkate alınamaz.

Rup bazı sınırlamalar içerir. Örneğin, kütüphane hazırlama ve dizilemeden önce RNA bütünlüğünün doğrudan ölçümleri gerektiğinden, RNA bozulmasını doğrudan test etmez. Bununla birlikte, geneBodyCoverage.py ve tin.py RSeQC modüllerini kullanarak RNA bozulmasını tanımlamak için bir komut dosyası, bu boru hattının deposuna dahil edilmiştir. Ayrıca Rup, GC içeriğini ve transkript uzunluğu yanlılığını test etmez. Referans genom boyutu şu anda 4 Gb ile sınırlıdır ve deneyimlerimize göre, okuma haritalama modülü, Ek Dosyada bir haploid ile diploid genom versiyonlarının karşılaştırılmasında örneklendiği gibi, poliploidlerde çoklu eşlenmiş okumaları olduğundan fazla tahmin eder. Rup çıktı kalitesi büyük ölçüde referans genomun ve açıklamanın kalitesine bağlıdır, öyle ki eksik veya oldukça parçalanmış bir genom dizisi, eşlenmemiş okumaların fazla tahmin edilmesine yol açabilir. Ayrıca, eksik gen modeli açıklaması, atanmamış okumaların olduğundan fazla tahmin edilmesine yol açar. Genom dizisinin ve gen açıklamasının eksiksizliği ve kopyası, BUSCO²⁰ gibi araçlarla çıkarılabilir.

Rup, RNA-seq deneylerinde gerekli olan tüm ilk kalite kontrol adımları için bağımsız bir araçtır. RSeQC ve RNA-SeQC'den farklı olarak, ham dizileme okumalarının ön işlenmesi gerekli değildir. RNA-SeQC ile dizileme kalitesi analizi yapılamaz ve örnek korelasyon ısı haritaları RSeQC ve RNA-SeQC tarafından hesaplanmaz (Tablo 1). Ayrıca, normalleştirme çıktısı RSeQC'de FPKM'dedir ve çıktı görselleştirme, RSeQC ve RNA-SeQC'deki tüm modüller için standart olarak uygulanmamaktadır. Bu nedenle, düşük okuma sayısı, kırpılmış okumaların yüksek oranı ve yinelenen aykırı değerlerin tanımlanmasını içeren iki alternatif araçta çeşitli kalite kontrol sorunları test edilmez. Rup, RSeQC ve RNA-SeQC'nin bir karşılaştırması Ek Dosya 1'de mevcuttur. Ayrıca Rup, RNA-SeQC veya RSeQC'ye tamamlayıcı olarak kullanılabilir, örneğin, bu boru hattı tarafından üretilen sıralanmış BAM dosyaları bu araçlar için girdi olarak kullanılabilir.

Şu anda Rup, az sayıda örnek için optimize edilmiştir, ancak kalite kırpma ve okuma eşleme gibi en çok zaman alan adımlar, örneğin bir bilgisayar kümesinde veya bir bulut altyapısında önceden hesaplanabilir. 4 Gb'den büyük genomlar için bu zorunludur, çünkü Rsubread hizalayıcı 4 Gb'den küçük genomlarla sınırlıdır. rRNA açıklaması, yüksek oranda korunmuş rRNA genlerini tanımlayan barrnap ile yapılır. Deneyimlerimize göre, rRNA gen açıklaması kapsamlılık gerektirmez, çünkü bazı anormal rRNA genleri eksik olsa bile, veri kümesindeki rRNA varlığına genel bir bakış, veri kümesi kalite değerlendirmesi için yeterlidir. Rup'un gelecekteki sürümleri, çalışma süresini azaltmak için sözde hizalama aracı salmon²¹ gibi farklı bir eşleme yöntemine geçebilir. Ayrıca, Rup'a GC yanlılığı veya uzunluk yanlılığı düzeltmesi gibi daha fazla normalleştirme ve düzeltme yöntemi eklenebilir.

Özetle Rup, RNA-seq veri analizinin güvenilirliğini ve tekrarlanabilirliğini sağlamak için gerekli bilgileri sağlar. Bu boru hattı, ana RNA-seq kalite ölçümlerini kapsamlı bir şekilde bildirir ve aşağı akış analizlerinde doğrudan kullanım için çıktı dosyaları üretir. Minimum biyoinformatik bilgisine sahip araştırmacıların birincil sıralama veri kalitesini sezgisel görselleştirme ile değerlendirmeleri için bağımsız bir araç olarak tasarlanmıştır.

Malzemeler

Bu makalede kullanılan malzemelerin listesi

Ad	Şirket	Katalog numarası	Yorumlar
fastqcr	R	0.1.3
FUJITSU ESPRIMO D958 Masaüstü	FUJITSU		boru hattı Ubuntu Linux 24.02 (32 GB RAM) üzerinde geliştirildi ve test edildi
getopt	R	1.20.4
ggplot2	R	3.5.2
MacBook Pro	Apple		boru hattı maxOS 15.4.1 (16 GB RAM) üzerinde test edildi
Pheatmap haritası	R	1.0.13
R		4.4.3
Reshape2	R	1.4.4
rfastp	Biyoiletken	1.16.0
rsamtools	Biyoiletken	2.22.0
rsubread	Biyoiletken	2.20.0