يسمح هذا البروتوكول بمراقبة الجودة الأولية لتجارب تسلسل الحمض النووي الريبي لعلماء الأحياء في المختبرات الرطبة ذوي الخبرة المحدودة في المعلوماتية الحيوية.
مقالة منهجية
يسمح هذا البروتوكول بمراقبة الجودة الأولية لتجارب تسلسل الحمض النووي الريبي لعلماء الأحياء في المختبرات الرطبة ذوي الخبرة المحدودة في المعلوماتية الحيوية.
غالبا ما تتطلب الأساليب الحديثة في علوم النبات الجزيئي تجارب تسلسل الحمض النووي الريبي الجماعي ، على سبيل المثال ، لتتبع التغييرات العالمية في النسخ عند المعالجات أو لتحديد المكونات الرئيسية للمسارات التنظيمية. وبالتالي ، تعتمد مجالات متنوعة من علوم النبات على بيانات تسلسل الحمض النووي الريبي عالية الجودة والقابلة للتكرار للتقدم العلمي. ومع ذلك ، من خلال تجربتنا ، غالبا ما تكون المعرفة حول وتطبيق تدابير مراقبة الجودة في مجموعات بيانات تسلسل الحمض النووي الريبي. هنا ، نقدم Rup (خط أنابيب تقييم قابلية الاستخدام RNA-seq) لمراقبة جودة بيانات RNA-seq المجمعة ، لتحليلات التعبير الجيني اللاحقة ، والتي تكون قائمة بذاتها وقابلة للتطبيق بسهولة لعلماء الأحياء في المختبرات الرطبة الذين لديهم معرفة أساسية ب R. يساعد Rup على التمييز بين بيانات التسلسل عالية الجودة والمناسبة لتجارب التعبير الجيني النهائي ، وتلك غير المناسبة لمزيد من التحليل العام. يتضمن Rup اختبارات للعديد من المشكلات الشائعة مثل عدم كفاية أرقام القراءة أو رسم الخرائط ، وتحديد التلوثات ، والتقدير الكمي لكسور الحمض النووي الريبي في إجمالي بيانات تسلسل الحمض النووي الريبي ، وتكرار اختبار التشابه ، واستخدام البيانات الحقيقية للعرض التوضيحي وتقديم تصور بديهي. يوفر Rup مجموعة من الأدوات لتحديد أوجه القصور التجريبية قبل تحليل النسخ القياسي ، وبالتالي تحسين جودة البيانات للباحثين الفرديين والميدان. وهذا يعزز الثقة في تحليل بيانات تسلسل الحمض النووي الريبي بالجملة ويوفر أساسا للإرشادات المستقبلية التي تحدد الحد الأدنى من معايير مراقبة الجودة، وبالتالي تحسين موثوقية وشفافية بيانات تسلسل الحمض النووي الريبي المنشورة. تتوفر بيانات Rup والاختبار في https://github.com/oliverrupp/rup.
تستجوب تجارب النسخ (RNA-seq) بشكل شامل توقيعات النسخ التي تشكل الأنماط الظاهرية. أصبح هذا النهج لا يمكن الاستغناء عنه في علم الوراثة الجزيئية النباتية لتحديد الجينات المفردة أو المعبر عنها بشكل مشترك والعمليات البيولوجية المتضمنة ، على سبيل المثال ، في التطوير ، أو تفاعل مسببات الأمراض النباتية ، أو مقاومة الإجهاد اللاأحيائي1 ، 2 ، 3. أدت التطورات الحديثة في تقنية تسلسل الحمض النووي الريبي إلى زيادة الخصوصية ومكنت من اكتشاف الأشكال الإسوية والمتغيرات المختلفة بدقة أحادية القاعدة ، مما يسمح بتحديد اختلافات التسلسل من إندل أكبر إلى تعدد أشكال النوكليوتيدات أحادية النوكليوتيدات (SNPs). تتميز البيانات التي تم الحصول عليها بواسطة تسلسل الحمض النووي الريبي بنطاق ديناميكي واسع ، مما يسهل الكشف عن كل من النصوص المعبر عنها الوفيرة والمنخفضة ، وتتطلب ظروفا تجريبية مناسبة للاتساق. بالإضافة إلى ذلك ، تعد مجموعات بيانات RNA-seq كبيرة وبالتالي فهي مكثفة حسابيا للتحليل والتخزين ، مما يستلزم إدارة فعالة للبيانات وموارد حوسبة واسعة النطاق. يتطلب RNA-seq مدخلات عالية الجودة في كل خطوة من خطوات سير العمل ، حيث يمكن أن تنتشر نقاط الضعف في أي مرحلة النتائج وتعرضها للخطر. سيؤدي ضعف سلامة الحمض النووي الريبي أو المشكلات الفنية أثناء إعداد المكتبة إلى التحيز وتقليل الدقة. يمكن أن تختلف معلمات الحصول على قراءات أولية عالية الجودة بين مرافق تسلسل الجيل التالي (NGS).
وفقا لتجربتنا ، نوصي باستخدام الحمض النووي الريبي فقط برقم تكامل الحمض النووي الريبي (RIN) أعلى من 7 ، وهو ما يدل على بنية mRNA سليمة إلى حد كبير ، كمدخلات لإعداد مكتبة RNA-seq. من أجل التسلسل الناجح ، يلزم ما يقرب من 2 ميكروغرام من إجمالي الحمض النووي الريبي بتركيز 50-200 نانوغرام / ميكرولتر لبروتوكولات إعداد المكتبة القياسية. يجب تأكيد نقاء الحمض النووي الريبي بنسبةOD 260/280 بين 1.8 و 2.1 ونسبة OD260/230 أكبر من 1.5 باستخدام مقياس الطيف الضوئي. قد يؤدي عمق التسلسل غير الكافي أو معدلات الخرائط المنخفضة أو المحاذاة الخاطئة للجينوم المرجعي إلى زيادة تشويه التعبير الجيني وقياس الربط الكمي. علاوة على ذلك ، فإن التصميم التجريبي غير الكافي ، مثل عدم التمييز بين عينات الأنسجة أو العلاجات المختلفة والتباين الكبير بين التكرارات ، يمكن أن يؤدي إلى حدوث ضوضاء وتقليل قابلية التكاثر. في حين أن العديد من المختبرات تحلل النسخ بشكل روتيني ، غالبا ما لا يتم الإبلاغ عن ضوابط الجودة الصارمة لخطوات التحليل الأساسية الأولى أو قد تكون غائبة تماما عن المنشورات. قد يؤدي هذا إلى الإفراط في تفسير النتائج المستمدة من تحليل النسخ ، وبالتالي إلى نتائج غير قابلة للتكرار.
هنا ، نقدم سير عمل لمراقبة الجودة للخطوات الأولية المطلوبة لتحليلات النسخ عالية الجودة لملفات تعريف mRNA لقياس التعديلات النسخية. هدفنا هو تمكين علماء الأحياء في المختبرات الرطبة ذوي المعرفة المحدودة في المعلوماتية الحيوية من تقييم بيانات النسخ الأولية الخاصة بهم. Rup (الشكل 1) متاح للباحثين الذين هم على دراية بالمعرفة الأساسية ل R. سيوفر تنفيذ سير العمل المقدم هنا للباحثين فهما مفصلا لبياناتهم الأولية ، بما في ذلك قيودهم المحتملة للتحليل اللاحق. على حد علمنا ، لا يوجد حتى الآن خط أنابيب عملي لتقييم بيانات النسخ الأولية جنبا إلى جنب مع إرشادات للتمييز بين البيانات عالية الجودة ومنخفضة الجودة.

الشكل 1: سير عمل تسلسل الحمض النووي الريبي في خط أنابيب مراقبة جودة السيليكو . يتم اشتقاق ملفات الإدخال لمراقبة الجودة من بيانات تسلسل الحمض النووي الريبي التي تم إنشاؤها بواسطة تسلسل مواد المصنع ، بالإضافة إلى مجموعات البيانات المتاحة للجمهور. يقوم خط أنابيب R المقدم بتقييم جودة التسلسل من خلال ثلاثة طرق رئيسية: جودة التسلسل ، وجودة رسم الخرائط ، وجودة التكرار. يتم حساب مقاييس الجودة المختلفة ، ويتم تصور النتائج الإحصائية باستخدام حزم R الشائعة مثل ggplot2 و pheatmap (على سبيل المثال ، المخططات الشريطية وخرائط الحرارة الرجاء النقر هنا لعرض نسخة أكبر من هذا الرقم.
تتطلب مسارات التقييم التي تم الإبلاغ عنها سابقا بيانات تمت معالجتها مسبقا (على سبيل المثال، قراءة المحاذاة كملفات .bam)، ولا تغطي جميع المقاييس، أو لم تعد تحتفظ بها4،5،6. تكمن ميزة سير العمل المقدم هنا في شموليته من خلال دمج مشكلات مراقبة الجودة الأكثر شيوعا في خط أنابيب واحد. علاوة على ذلك ، نقدم أمثلة للبيانات عالية الجودة المناسبة لجميع تطبيقات التحليل النهائي ، ولكن أيضا أمثلة للبيانات منخفضة الجودة ، ونناقش قيودها المحددة لمزيد من التحليل. تصف العديد من التقارير المنشورة سابقا الغرض من أدوات تحليل تسلسل الحمض النووي الريبي وتقدم تقييمات مقارنة لأدائها7،8. ومع ذلك ، فإن معايير مراقبة الجودة والإبلاغ عن المنهجية في تسلسل الحمض النووي الريبي ليست موحدة وتؤدي إلى تفاقم قابلية التكاثر والتفسيرات ذات المغزى البيولوجي لتجارب النسخ
يدمج Rup الأدوات القياسية عالية الجودة وتحليل مراقبة الجودة وتصور النتائج. يتطلب Rup قراءات تسلسل أولية وجينوم مشروح كمدخلات ويعمل على أنظمة Mac OS و Linux و "نظام Windows الفرعي لنظام التشغيل Linux" (WSL) على أنظمة Windows. إنه يدمج أدوات لجودة التسلسل ويحدد رسم خرائط القراءة إلى الجينوم ، ويتضمن قياس محتوى الحمض النووي الريبي في العينات ، ويقدم ارتباط عينة لتقييم الارتباط المتماثل. تم الحصول على بيانات الاختبار باستخدام التشريح الدقيق بالليزر لأنسجة النسيج الإنشائي المركزي لمرحلتين مختلفتين من الأنواع النباتية Eschscholzia californica ، وهو بروتوكول إدخال منخفض للغاية لإعداد المكتبة ، وتم تسلسله على Novaseq 6000.
يمكن تشغيل Rup كبرنامج نصي واحد مع الحد الأدنى من ملفات الإدخال المطلوبة. يتطلب خط الأنابيب ، كما هو موضح هنا ، بيانات تسلسل الحمض النووي الريبي المقترنة من Illumina. يتم أيضا إيداع مسار معدل للتسلسل أحادي الطرف بنفس خطوات التحليل في مستودع GitHub. فقط تسلسل الجينوم كملف fasta ، والنموذج الجيني والتعليقات التوضيحية ل rRNA كملفات gtf ، والتسلسل الخام يقرأ كملفات fastq.gz. تأكد من توفير ملفات الجينوم والتعليقات التوضيحية ك genome.fa و annotation.gtf و rRNA.gtf في مجلد محدد في المتغير reference_folder. عندما يتم حفظ جميع ملفات التسلسل كملفات .fq.gz في المتغير read_file_folder ، يمكن تشغيل Rup على النحو التالي:
1. الاستعدادات:
BiocManager::install(c("getopt", "ggplot2", "reshape2", "pheatmap", "fastqcr", "Rfastp", "Rsubread", "Rsamtools"))
# 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)
n_threads <- 8 # the number of available CPU cores
bamSortMemory <- "1024" # maximum memory for bam file sorting
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))
2. تقييم جودة التسلسل
ملاحظة: ستقوم هذه الخطوة بإنشاء مخطط شريطي يوضح عدد القراءات قبل وبعد قص كل عينة.
# 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)
# 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)
}
# 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)
# 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)

الشكل 2: نتائج التسلسل والتشذيب. أعداد القراءة قبل (أحمر) وبعد التشذيب (أخضر). تحتوي s2_r1 العينة بالفعل على عدد قراءة منخفض قبل التشذيب ، بينما أزال التشذيب جزءا كبيرا من s2_r2 القراءات للعينات. تظهر جميع العينات الأخرى خسارة مقبولة للقراءات من التشذيب. الرجاء النقر هنا لعرض نسخة أكبر من هذا الرقم.
3. جودة رسم الخرائط
ملاحظة: تحسب هذه الخطوة المخططات الشريطية للتمييز بين القراءات المفردة المعينة (على سبيل المثال ، النصوص من جينات ترميز البروتين) من القراءات المتعددة المعينة (على سبيل المثال ، قراءات الحمض النووي الريبوزي الريبي) والقراءات غير المعينة (على سبيل المثال ، الملوثات).
# 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)
# 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"))
}
# 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)
# 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)
# 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"
# 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"
# 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)
# 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)

الشكل 3: اقرأ نظرة عامة على رسم الخرائط والإحصائيات. تظهر معظم العينات عددا كبيرا من القراءات المخصصة (الجانب الأيسر ، البني) وعددا منخفضا من قراءات الحمض النووي الريبي الريبي (الجانب الأيمن ، الأزرق). يحتوي s2_r3 العينة على كمية عالية بشكل غير عادي من القراءات متعددة التعيينات (الجانب الأيسر ، الأخضر) المقابلة لرقم قراءة rRNA مرتفع (الجانب الأيمن). تظهر s2_r4 العينة عددا كبيرا من القراءات غير المعينة بالاقتران مع رقم قراءة rRNA المتوقع ، مما يشير إلى التلوث بقراءات من كائن حي آخر. الرجاء النقر هنا لعرض نسخة أكبر من هذا الرقم.
4. تكرار الجودة.
# 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")

الشكل 4: تصنيف عدد قراءة الجينات. يتم تصنيف جميع الجينات في العينة إلى واحدة من خمس فئات بناء على عدد القراءات المخصصة. يظهر باللون الأحمر عدد الجينات دون تعيين أي قراءات. تحتوي العينات التي تحتوي على عدد منخفض من إجمالي القراءات المخصصة (من s2_r1 إلى s1_r4) على عدد أكبر من الجينات مع 10 إلى 100 قراءة مخصصة وعدد أقل من الجينات بأكثر من 1000 قراءة. قد يتم تفويت الجينات ذات التعبير المنخفض في هذه العينات. الرجاء النقر هنا لعرض نسخة أكبر من هذا الرقم.
يتم تنفيذ خط الأنابيب بالكامل كبرنامج نصي R وتم اختباره على أنظمة تشغيل Linux و Mac OS. يمكن لمستخدمي Windows استخدام نظام Windows الفرعي لنظام التشغيل Linux (WSL). تتوفر التعليمات البرمجية وبيانات الاختبار كمستودع GitHub: https://github.com/oliverrupp/rup. تتوفر بيانات التسلسل في إطار مشروع ENA EBI PRJEB96400.
تم إنشاء عشر عينات بشكل مصطنع من عينتين حقيقيتين لتوضيح المشكلات المتنوعة التي قد تتم مواجهتها أثناء مراقبة جودة تسلسل الحمض النووي الريبي السائب باستخدام Rup للتحليل. تم تصميم s2_r1 العينة لإظهار عدد قراءة إجمالي منخفض ، وتحتوي s2_r2 العينة على جزء كبير من القراءات منخفضة الجودة ليتم التخلص منها بواسطة عملية التشذيب. تتضمن s2_r3 العينة جزءا كبيرا من قراءات الحمض النووي الريبوزي الريبوزي وتتضمن s2_r4 العينة قراءات الملوثات التي فشلت في التعيين إلى الجينوم المرجعي. تم تبديل أسماء العينات s1_r5 و s2_r5 لتوضيح ارتباط التكرار المنخفض.
يتيح قسم البروتوكول 2 تحديد العينات ذات الأرقام المنخفضة للقراءة قبل أو بعد التشذيب (الشكل 2). يظهر مخطط الشريط رقم القراءة الأدنى في s1_r1 العينة قبل التشذيب وبعده. هنا ، كان رقم القراءة الأولي منخفضا. يشير عدد القراءة المنخفض للعينة s1_r2 بعد التشذيب إلى كمية كبيرة من تسلسلات المحولات ، وأخطاء التسلسل ، وتسلسلات التمهيدي ، وامتدادات poly-A / T المتسلسلة التي تمت إزالتها أثناء عملية التشذيب. قد يؤدي تدهور الحمض النووي الريبي المدخل أيضا إلى تقليل عدد القراءات عالية الجودة.
يحدد القسم 3 من البروتوكول المشكلات في تعيينات القراءة. يوضح الشكل 3 العدد المرتفع من القراءات متعددة الخرائط في s2_r3 العينة (أخضر) ، بالإضافة إلى العدد الكبير من قراءات الحمض النووي الريباسي. يتضح تلوث s2_r4 العينة من خلال جزء كبير من القراءات التي لم يتم تعيينها على الجينوم المرجعي (الوردي). لا ترتبط هذه القراءات بالحمض النووي الريبوزي الريبي ولكن بتسلسلات من كائن حي غير مستهدف. يوضح الشكل 4 المشكلات العامة المرتبطة بانخفاض عدد القراءات المعينة في النسخ النسخية. في العينات ذات المعدل المنخفض من القراءات التي تم تعيينها بشكل فريد للجينوم المرجعي (s2_r1 إلى r4) ، يحتوي حوالي ثلث الجينات فقط على أكثر من 100 قراءة معينة بينما في العينات الأخرى ، يقع حوالي نصف الجينات في هذه الفئات. قد لا يتم العثور على قراءات للجينات ذات التعبير المنخفض جدا في العينات s2_r1 إلى r4 ، وبالتالي ، فإن تحليل التعبير المقارن مع هذه العينات سيكون غير موثوق به إلى حد كبير ويجب تجنبه.
يمكن استخدام القسم 4 من البروتوكول لتحديد القيم المتطرفة المتماثلة. من المتوقع أن تظهر التكرارات لنفس العينة / الحالة / الأنسجة ارتباطا أعلى بين بعضها البعض مقارنة بالتكرار من عينات / ظروف / أنسجة أخرى. يوضح الشكل 5 خريطة حرارية ارتباطية للعينتين (S1 و S2) مع خمسة مكررات لكل منهما. تظهر مخططات التغصنات الموجودة في الأعلى وعلى جانب المخطط مجموعتين مع خمسة مكررات في كل مجموعة. تحتوي المجموعة اليسرى على أربعة نسخ مكررة من العينة 1 ونسخة مكررة واحدة من العينة 2 (s2_r5) ، وتحتوي المجموعة اليمنى على أربعة نسخ مكررة من العينة 2 ونسخة مكررة واحدة من العينة 1 (s1_r5). في هذه الحالة، عندما يتم تبديل أسماء العينات s1_r5 و s2_r5 مرة أخرى، تحتوي كل مجموعة على جميع النسخ المتماثلة لعينة واحدة، مما قد يشير إلى خطأ في تسمية النسخ المتماثل. قد تكون الأسباب الأخرى لعدم تجمع النسخ المتماثلة معا هي عدم التمايز بين العينات / الظروف / الأنسجة ، أو تجميع التكرارات لأسباب تتعلق بجودة التسلسل فقط. قد يحدث هذا الأخير عندما تشكل جميع النسخ المتماثلة ذات الانخفاض المنخفض بشكل استثنائي أو تلك التي تحتوي على عدد قراءة مرتفع بشكل استثنائي بعد الاقتطاع والتعيين نظام مجموعة.

الشكل 5: عينة خريطة حرارية للارتباط. تستند خريطة التمثيل اللوني لارتباط العينة إلى قيم TPM المحولة في السجل2 وتظهر مجموعتين متميزتين من خمس عينات لكل منهما. من المتوقع أن تظهر التكرارات البيولوجية / التقنية ارتباطا أعلى بين بعضها البعض مقارنة بتلك الموجودة في الأنسجة / العلاجات الأخرى. تحتوي المجموعة اليسرى على أربعة نسخ مكررة من العينة 1 ونسخة مكررة واحدة من العينة 2 (s2_r5) ، وتحتوي المجموعة اليمنى على أربعة نسخ مكررة من العينة 2 ونسخة مكررة واحدة من العينة 1 (s1_r5). يجب التحقق من التكرارات الفردية بحثا عن تبديل العينة المحتملة أو تأثيرات الدفعات لشرح تجميعها في خريطة التمثيل اللوني. الرجاء النقر هنا لعرض نسخة أكبر من هذا الرقم.
الجدول 1: مقارنة خطوط أنابيب تقييم جودة RNA-seq. للحصول على تحليل مفصل ، راجع الملف التكميلي 1. الرجاء النقر هنا لتنزيل هذا الجدول.
الملف التكميلي 1: اختيار المعلمات والمراجع. تظهر التحليلات تأثير اختيار المعلمات والمراجع على النتائج الإجمالية لمراقبة الجودة. الرجاء النقر هنا لتنزيل هذا الملف.
تعتمد جودة تحليلات التعبير الجيني التفاضلي بشكل كبير على عاملين: عدد القراءات المتسلسلة في كل عينة وتكرار16 وعدد التكرارات لكل عينة17,18. هنا ، نقدم خط أنابيب Rup سهل الاستخدام لتحديد عدد القراءات المناسبة لقياس كمية التعبير الجيني في كل تكرار. تسمح المقاييس المختلفة للباحثين بفهم سبب إظهار النسخ المتماثلة لأرقام قراءة منخفضة مخصصة وتحديد المشكلات في ارتباط العينة. على الرغم من أن Rup تم تطويره لتقييم جودة عينات الحمض النووي الريبي النباتي ، إلا أنه مناسب بنفس القدر للكائنات حقيقية النواة الأخرى ، ولا يتطلب مزيدا من التعديلات للعينات غير النباتية (الملف التكميلي 1).
تحدد الخطوة الأولى من Rup إجمالي عدد القراءة قبل وبعد قص القراءة. يحدد عدد القراءات المتسلسلة عدد الجينات التي يمكن اكتشافها ، وإذا كان رقم القراءة منخفضا جدا ، فستظل العديد من الجينات المعبر عنها تفاضليا غير مكتشفة. قد يشير جزء كبير من القراءات التي تم التخلص منها أثناء عملية تشذيب وتصفية الجودة إلى تدهور الحمض النووي الريبي ، مما قد يؤدي إلى التقليل من تقدير التعبير عن الجينات شديدة التحلل. ومع ذلك ، فإن ما إذا كان يمكن استخدام عينة للتطبيقات النهائية يعتمد على الكائن الحي المستهدف والهدف من الدراسة. على سبيل المثال ، لالتقاط التعبير عن معظم الجينات النباتية ، في تجربتنا ، هناك حاجة إلى 30 مليون إلى 50 مليون قراءة ، بينما بالنسبة للفطريات ، قد يكون 10 ملايين فقط كافيا. وبالتالي ، لن يحدد Rup عتبات الجودة لاستبعاد العينات الإشكالية بل سيوفر مقاييس لتحديد المشكلات المتنوعة التي قد تواجهها.
تقيم الخطوة الثانية (جودة رسم الخرائط) دقة وموثوقية محاذاة القراءة للجينوم المرجعي ، مع تحديد مدى ملاءمتها للتحليلات النهائية. لا يمكن استخدام جميع القراءات المتسلسلة لحساب وفرة الجينات. يتم تجاهل القراءات التي تفشل في محاذاة الجينوم أو النسخ أو تقرأ تلك الخريطة إلى مواقع متعددة على الجينوم في معظم الحالات19 ولا تساهم في عدد القراءة الإجمالي. يمكن استخدام نتائج الخطوة الثانية من خط الأنابيب لفهم سبب عدم استخدام القراءات في حساب الوفرة. يمكن أن يشير العدد الكبير من القراءات غير المعينة إلى التلوث أثناء استخراج الحمض النووي الريبي (على سبيل المثال ، بمسببات الأمراض النباتية أو العاشبة) أو جينوم مرجعي غير مكتمل. قد تؤدي النماذج الجينية غير المكتملة للجينوم المرجعي إلى عدد كبير من قراءات "عدم وجود ميزة". يمكن أن يكون سبب جزء كبير من القراءات متعددة التعيينات هو عدد كبير من جينات الرنا الريبوزي في المكتبة مما يشير إلى عدم كفاية إزالة الحمض النووي الريبي أثناء عملية إعداد مكتبة التسلسل (في حالة اشتقاق القراءات المتعددة المتعددة حقا من الحمض النووي الريبي ، يمكن تحليلها بشكل أكبر من خلال توفير ملف تعليق توضيحي للحمض النووي الريبي إلى خط الأنابيب ، والذي يحسب تلقائيا عدد قراءات الحمض النووي الريبي المحتملة).
مقياس الجودة الثالث الذي لا يقل أهمية عن ذلك هو الارتباط بين التكرارات لنفس العينة. بشكل عام ، يجب أن يكون الارتباط بين التكرارات لنفس العينة أعلى من الارتباط بين التكرارات لعينات مختلفة. يمكن أن يشير الارتباط المنخفض بين التكرارات إلى تباين بيولوجي كبير بين التكرارات ، أو التشابه العالي بين العينات / الظروف / الأنسجة ، أو تأثيرات الدفعات الأخرى ، أو حتى تبديل العينات أو وضع العلامات الخاطئة. يحسب Rup الارتباطات الزوجية بين جميع العينات وينتج خريطة حرارية مجمعة لارتباطات العينة. بالإضافة إلى ذلك ، يمكن حساب تحليل المكون الرئيسي (PCA) على العينات لتحديد تأثيرات الدفعات المحتملة التي يجب الاعتراف بها في المزيد من التحليلات النهائية. على الرغم من أنه يمكن تصحيح تأثيرات الدفعات في التحليلات النهائية ، فقد يكون من الأفضل إزالة العينات التي بها مشاكل إذا كان الفرق المقاس مرتفعا جدا.
المواد المدخلة لها تأثير مباشر على جودة التسلسل: يعد أخذ عينات الأنسجة واستخراج الحمض النووي الريبي وإعداد المكتبة خطوات حاسمة لتقليل مشكلات جودة تسلسل الحمض النووي الريبي. يجب الحفاظ على ظروف متسقة ، إن أمكن ، على سبيل المثال ، باستخدام غرف النمو. يجب تنفيذ تدابير مكافحة الآفات في الوقت المناسب ، حيث يجب اختيار الأفراد الأصحاء فقط لأخذ العينات. ينصح أيضا بجمع العينات في نفس اليوم و / أو في نفس الوقت لتقليل تباين النسخ اليومي. تتوفر العديد من مجموعات استخراج الحمض النووي الريبي ، ويمكن أن يؤدي اختيار مجموعة مناسبة للأنواع المستهدفة إلى تحسين جودة mRNA. يجب أن تتضمن بروتوكولات إعداد المكتبة خطوات لإثراء polyA + RNA لتقليل جزء الحمض النووي الريباسي (rRNA).
يمكن استخدام Rup كخطوة أولية لمراقبة الجودة في تحليل التعبير الجيني التفاضلي. مراقبة الجودة هذه مطلوبة لعدة أسباب حاسمة ، لأنها تضمن جودة بيانات الإدخال (تحدد النسخ المتماثلة / العينات منخفضة الجودة ، ويمكنها تعيين عتبات الجودة لعمق القراءة ومعدلات رسم الخرائط) وتحدد المشكلات الفنية مثل أخطاء التسلسل أو تأثيرات الدفعات أو التسمية الخاطئة. إذا كان إدخال الحمض النووي الريبي ذا جودة كافية ، فيمكن أن يساعد Rup عن طريق قص القراءات منخفضة الجودة أو تحديد العينات ذات العلامات الخاطئة. ومع ذلك ، لا يمكن ل Rup تعويض جودة الحمض النووي الريبي للإدخال الضعيف. في حالات الحمض النووي الريبي منخفض الجودة أو بيانات التسلسل الخاطئة ، قد يكون من الضروري إعادة جمع العينات وضبط بروتوكول استخراج الحمض النووي الريبي و / أو تكرار التسلسل. وينطبق الشيء نفسه على العينات التي تحتوي على جزء كبير من القراءات غير المعينة التي قد تكون ناجمة عن التلوث. ومع ذلك ، من تجربتنا ، لا يمكن تكرار جميع عمليات استخراج الحمض النووي الريبي ببساطة بسبب نقص توافر المواد. في هذه الحالة ، فهي غير مناسبة لبعض التطبيقات النهائية: لا ينبغي تحليل العينات التي تحتوي على أعداد منخفضة من القراءات الفردية في تحليلات التعبير الجيني التفاضلي ، لكنها لا تزال تحتفظ بمعلومات لتحليل وجود النص. في هذه الحالة ، لا يمكن النظر في عدم وجود نسخ في الأنسجة / العلاجات / الحالات لتحليل وقياس وفرة النص.
يتضمن Rup بعض القيود. على سبيل المثال ، لا يختبر بشكل مباشر تدهور الحمض النووي الريبي حيث أن القياسات المباشرة لسلامة الحمض النووي الريبي مطلوبة قبل إعداد المكتبة وتسلسلها. ومع ذلك، يتم تضمين برنامج نصي لتحديد تدهور الحمض النووي الريبي باستخدام وحدات RSeQC geneBodyCoverage.py و tin.py في مستودع هذا التنابيب. علاوة على ذلك ، لا يختبر Rup محتوى GC وتحيز طول النص. يقتصر حجم الجينوم المرجعي حاليا على 4 جيجابايت ، ووفقا لتجربتنا ، فإن وحدة رسم الخرائط المقروءة تبالغ في تقدير القراءات متعددة التعيينات في polyploids ، كما هو موضح في الملف التكميلي في مقارنة رسم خرائط القراءة إلى إصدارات الجينوم أحادية الصيغة الصبغية مقابل إصدارات الجينوم ثنائية الصبغيات من E. californica ، وبالتالي فإن الجينوم أحادي الصيغة الصبغية هو المدخل المفضل لخط الأنابيب هذا. تعتمد جودة إخراج Rup إلى حد كبير على جودة الجينوم المرجعي والتعليق التوضيحي ، بحيث يمكن أن يؤدي تسلسل الجينوم غير المكتمل أو المجزأ للغاية إلى المبالغة في تقدير القراءات غير المعينة. علاوة على ذلك ، يؤدي التعليق التوضيحي غير المكتمل لنموذج الجينات إلى المبالغة في تقدير القراءات غير المخصصة. يمكن الاستدلال على اكتمال وازدواجية تسلسل الجينوم والتعليق التوضيحي للجينات باستخدام أدوات مثل BUSCO20.
Rup هي أداة قائمة بذاتها لجميع خطوات مراقبة الجودة الأولية الضرورية في تجارب تسلسل الحمض النووي الريبي. على عكس RSeQC و RNA-SeQC ، فإن المعالجة المسبقة لقراءات التسلسل الخام غير مطلوبة. لا يمكن إجراء تحليل جودة التسلسل باستخدام RNA-SeQC ، ولا يتم حساب خرائط الحرارة لارتباط العينة بواسطة RSeQC و RNA-SeQC (الجدول 1). علاوة على ذلك ، يكون إخراج التطبيع في FPKM في RSeQC ولا يتم تنفيذ تصور الإخراج كمعيار لجميع الوحدات النمطية في RSeQC و RNA-SeQC. وبالتالي ، لم يتم اختبار العديد من مشكلات مراقبة الجودة في الأداتين البديلتين ، والتي تشمل انخفاض عدد القراءة ، وجزء كبير من القراءات المقتطعة وتحديد القيم المتطرفة المكررة. تتوفر مقارنة بين Rup و RSeQC و RNA-SeQC في الملف التكميلي 1. علاوة على ذلك ، يمكن استخدام Rup مكملا ل RNA-SeQC أو RSeQC ، على سبيل المثال ، يمكن استخدام ملفات BAM المصنفة التي ينتجها خط الأنابيب هذا كمدخل لهذه الأدوات.
حاليا، تم تحسين Rup لعدد صغير من العينات، ولكن يمكن حساب الخطوات الأكثر استهلاكا للوقت، مثل قص الجودة ورسم خرائط القراءة، مسبقا، على سبيل المثال، على مجموعة كمبيوتر أو بنية تحتية سحابية. بالنسبة للجينومات التي يزيد حجمها عن 4 جيجابايت ، يعد هذا إلزاميا ، نظرا لأن محاذاة Rsubread تقتصر على الجينومات الأصغر من 4 جيجابايت. يتم إجراء التعليق التوضيحي للحمض النووي الريبي باستخدام barrnap ، والذي يحدد جينات الحمض النووي الريبي المحفوظة للغاية. وفقا لتجربتنا ، لا يتطلب التعليق التوضيحي لجين الرنا الريبوزي الريبوزي شمولية ، لأنه حتى لو كانت بعض جينات الحمض النووي الريبي غير الطبيعية مفقودة ، فإن نظرة عامة إجمالية على وجود الحمض النووي الريبوزي الريبوزي في مجموعة البيانات كافية لتقييم جودة مجموعة البيانات. قد تتحول الإصدارات المستقبلية من Rup إلى طريقة تعيين مختلفة ، مثل أداة المحاذاة الزائفة salmon21 ، لتقليل وقت التشغيل. علاوة على ذلك ، يمكن إضافة المزيد من طرق التطبيع والتصحيح ، مثل تحيز GC أو تصحيح تحيز الطول ، إلى Rup.
باختصار ، يوفر Rup معلومات أساسية لضمان موثوقية وقابلية استنساخ تحليل بيانات تسلسل الحمض النووي الريبي. يقوم خط الأنابيب هذا بالإبلاغ بشكل شامل عن مقاييس جودة تسلسل الحمض النووي الريبي الرئيسية وينتج ملفات الإخراج للاستخدام المباشر في التحليلات النهائية. تم تصميمه كأداة قائمة بذاتها للباحثين الذين لديهم الحد الأدنى من المعرفة بالمعلوماتية الحيوية لتقييم جودة بيانات التسلسل الأولية الخاصة بهم مع تصور بديهي.
ليس لدى أصحاب البلاغ تضارب في المصالح للإعلان.
نحن نقدر المساعدة الفنية من قبل مرفق المعلوماتية الحيوية الأساسي في أستاذية بيولوجيا الأنظمة في JLU Giessen وتوفير موارد الحوسبة والدعم العام من قبل مركز خدمة BiGi (منحة BMFB 031A533) داخل المركز. شبكة NBI. تم تمويل العمل المقدم هنا من قبل منحة مؤسسة الأبحاث الألمانية (DFG) BE2547 / 24-1 إلى AB ، ونحن ممتنون أيضا للدعم المقدم من جامعة Justus Liebig Giessen ، ألمانيا.
| الاسم | الشركة | رقم فهرسي | التعليقات |
|---|---|---|---|
| فاستQCR | R | 0.1.3 | |
| FUJITSU ESPRIMO D958 سطح المكتب | فوجيتسو | تم تطوير واختبار خط الأنابيب على أوبونتو لينكس 24.02 (ذاكرة 32 جيجابايت) | |
| جيتوبت | R | 1.20.4 | |
| ggplot2 | R | 3.5.2 | |
| ماك بوك برو | آبل | تم اختبار خط الأنابيب على maxOS 15.4.1 (ذاكرة عشوائية 16 جيجابايت) | |
| خريطة الخيمة | R | 1.0.13 | |
| R | 4.4.3 | ||
| إعادة الشكل2 | R | 1.4.4 | |
| rfastp | الموصل الحيوي | 1.16.0 | |
| rsamtools | الموصل الحيوي | 2.22.0 | |
| rsubread | الموصل الحيوي | 2.20.0 |
طلب إذن لإعادة استخدام النص أو الأشكال في مقالة JoVE هذه
طلب إذن