Rup (צינור הערכת שימושיות RNA-seq) - בקרת איכות לניסויי RNA-seq בתפזורת באיקריוטים

סיכום

פרוטוקול זה מאפשר בקרת איכות ראשונית לניסויי RNA-seq עבור ביולוגים במעבדה רטובה עם ניסיון מוגבל בביואינפורמטיקה.

תקציר

גישות מודרניות במדעי הצמח המולקולרי דורשות לעתים קרובות ניסויי RNA-seq בכמויות גדולות, למשל, כדי לעקוב אחר שינויים גלובליים בטרנסקריפטומים לאחר טיפולים או כדי לזהות מרכיבי מפתח של מסלולים רגולטוריים. כתוצאה מכך, תחומים מגוונים של מדעי הצמח מסתמכים על נתוני RNA-seq איכותיים וניתנים לשחזור לצורך התקדמות מדעית. עם זאת, מניסיוננו, לעתים קרובות חסר ידע ויישום של אמצעי בקרת איכות במערכי נתונים של RNA-seq. כאן, אנו מציגים את Rup (צינור הערכת שימושיות RNA-seq) לבקרת איכות של נתוני RNA-seq בתפזורת, לניתוחי ביטוי גנים עוקבים, שהוא עצמאי וישים בקלות עבור ביולוגים במעבדה רטובה עם ידע בסיסי ב-R. Rup עוזר להבחין בין נתוני ריצוף באיכות גבוהה, המתאימים לניסויי ביטוי גנים במורד הזרם, לבין אלה שאינם מתאימים לניתוח כללי נוסף. Rup כולל בדיקות למספר בעיות נפוצות כגון מספרי קריאה או מיפוי לא מספיקים, זיהוי זיהומים, כימות שברי rRNA בנתוני ה-RNA-seq הכוללים, בדיקת דמיון שכפול, שימוש בנתונים אמיתיים להדגמה והצעת הדמיה אינטואיטיבית. Rup מספקת חבילת כלים לזיהוי חסרונות ניסיוניים לפני ניתוח טרנסקריפטום סטנדרטי, ובכך משפרת את איכות הנתונים עבור חוקרים בודדים והתחום. זה משפר את הביטחון בניתוח נתוני RNA-seq בתפזורת ומספק בסיס להנחיות עתידיות המגדירות קריטריונים מינימליים לבקרת איכות, ובכך משפר את האמינות והשקיפות של נתוני RNA-seq שפורסמו. נתוני RUP ובדיקות זמינים בכתובת https://github.com/oliverrupp/rup.

מבוא

ניסויי טרנסקריפטומיקה (RNA-seq) חוקרים באופן מקיף חתימות שעתוק המעצבות פנוטיפים. גישה זו הפכה לבלתי ניתנת להחלפה בגנטיקה מולקולרית של צמחים כדי לזהות גנים בודדים או משותפים ותהליכים ביולוגיים המעורבים, למשל, בהתפתחות, אינטראקציה עם פתוגנים בצמחים, או עמידות לעקה אביוטית^{1, 2}, ³. ההתקדמות האחרונה בטכנולוגיית RNA-seq הגדילה את הספציפיות ואפשרה זיהוי של איזופורמים וגרסאות שונות ברזולוציה של בסיס יחיד, מה שמאפשר זיהוי וריאציות רצף מאינדל גדול יותר לפולימורפיזמים של נוקלאוטידים בודדים (SNPs). הנתונים המתקבלים על ידי RNA-seq מאופיינים בטווח דינמי רחב, המאפשר זיהוי של תעתיקים בשפע ובביטוי נמוך, ודורשים תנאי ניסוי מתאימים לעקביות. בנוסף, מערכי הנתונים של RNA-seq הם גדולים ולכן אינטנסיביים מבחינה חישובית לניתוח ואחסון, הכרוכים בניהול נתונים יעיל ומשאבי מחשוב נרחבים. RNA-seq דורש קלט איכותי בכל שלב בזרימת העבודה, מכיוון שחולשות בכל שלב עלולות להתפשט ולפגוע בתוצאות. שלמות RNA ירודה או בעיות טכניות במהלך הכנת הספרייה יובילו להטיות ולהפחתת הדיוק. הפרמטרים להשגת קריאות גולמיות באיכות גבוהה יכולים להשתנות בין מתקני ריצוף מהדור הבא (NGS).

על פי הניסיון שלנו, אנו ממליצים להשתמש רק ב-RNA עם מספר שלמות RNA (RIN) מעל 7, המעיד על מבנה mRNA שלם ברובו, כקלט להכנת ספריית RNA-seq. לריצוף מוצלח, נדרשים כ-2 מיקרוגרם של RNA כולל בריכוז של 50-200 ננוגרם/מיקרוליטר עבור פרוטוקולי הכנת ספרייה סטנדרטיים. יש לאשר את טוהר ה-RNA על ידי יחס OD_260/280 בין 1.8 ל-2.1 ויחס OD_260/230 גדול מ-1.5 באמצעות ספקטרופוטומטר. עומק ריצוף לא מספיק, שיעורי מיפוי נמוכים או חוסר יישור לגנום הייחוס עלולים לעוות עוד יותר את ביטוי הגנים וכימות השחבור. יתר על כן, תכנון ניסוי לקוי, כגון חוסר הבחנה בין דגימות של רקמות או טיפולים שונים ושונות גבוהה בין שכפולים, יכול להכניס רעש ולהפחית את יכולת השחזור. בעוד שמעבדות רבות מנתחות טרנסקריפטומים באופן שגרתי, בקרות איכות מחמירות של שלבי הניתוח החיוניים הראשונים לרוב אינן מדווחות או עשויות להיעדר לחלוטין מהפרסומים. זה עלול להוביל לפרשנות יתר של תוצאות הנגזרות מניתוח טרנסקריפטום וכתוצאה מכך לתוצאות בלתי ניתנות לשחזור.

כאן, אנו מספקים זרימת עבודה לבקרת איכות של השלבים הראשוניים הנדרשים לניתוחי טרנסקריפטום באיכות גבוהה של פרופילי mRNA למדידת שינויים בתעתיק. מטרתנו היא לאפשר לביולוגים במעבדה רטובה עם ידע מוגבל בביואינפורמטיקה להעריך את נתוני הטרנסקריפטומיקה העיקריים שלהם. Rup (איור 1) נגיש לחוקרים שמכירים את הידע הבסיסי של R. ביצוע זרימת העבודה המוצגת כאן יספק לחוקרים הבנה מפורטת של הנתונים הראשוניים שלהם, כולל המגבלות הפוטנציאליות שלהם לניתוח עוקב. למיטב ידיעתנו, חסר עד כה צינור הערכת נתוני טרנסקריפטום ראשוני מעשי בשילוב עם הנחיות להבחנה בין נתונים באיכות גבוהה לנמוכה.

איור 1: זרימת העבודה של ה-RNA-seq בצנרת בקרת איכות סיליקו . קבצי הקלט לבקרת איכות נגזרים מנתוני RNA-seq שנוצרו על ידי ריצוף חומר צמחי, כמו גם מערכי נתונים זמינים לציבור. צינור ה-R המסופק מעריך את איכות הרצף באמצעות שלוש גישות עיקריות: איכות ריצוף, איכות מיפוי ואיכות שכפול. מדדי איכות שונים מחושבים, ותוצאות סטטיסטיות מוצגות באמצעות חבילות R נפוצות כגון ggplot2 ו-pheatmap (למשל, תרשימי עמודות ומפות חום). אנא לחץ כאן לצפייה בגרסה גדולה יותר של איור זה.

צינורות הערכה שדווחו בעבר דרשו נתונים מעובדים מראש (למשל, יישור קריאה כקבצי .bam), אינם מכסים את כל המדדים, או אינם מתוחזקים עוד ^4,5,6. היתרון של זרימת העבודה המוצגת כאן טמון במקיפות שלה על ידי איחוד בעיות בקרת האיכות הנפוצות ביותר לצינור יחיד. יתר על כן, אנו מספקים דוגמאות לנתונים באיכות גבוהה המתאימים לכל יישומי הניתוח במורד הזרם, אך גם דוגמאות לנתונים באיכות נמוכה, ודנים במגבלות הספציפיות שלהם להמשך ניתוח. מספר דוחות שפורסמו בעבר מתארים את המטרה של כלי ניתוח RNA-seq ומספקים הערכות השוואתיות של ביצועיהם ^7,8. עם זאת, תקני בקרת איכות ודיווח מתודולוגיה של RNA-seq אינם סטנדרטיים ומחמירים את יכולת השחזור והפרשנויות המשמעותיות מבחינה ביולוגית של ניסויי טרנסקריפטומיקה.

Rup משלבת כלים סטנדרטיים באיכות גבוהה, ניתוח בקרת איכות והדמיה של התוצאות. Rup דורש קריאות רצף גולמיות וגנום מבואר כקלט ופועל ב-Mac OS, Linux ו-"Windows Subsystem for Linux" (WSL) במערכות Windows. הוא משלב כלים לאיכות הריצוף ומכמת מיפוי קריאה לגנום, כולל מדידת תוכן rRNA בדגימות, ומציע מתאם דגימה להערכת מתאם שכפול. נתוני הבדיקה התקבלו באמצעות מיקרו-דיסקציה בלייזר של רקמת המריסטם המרכזית של שני שלבים שונים של מין הצמח Eschscholzia californica, פרוטוקול קלט נמוך במיוחד להכנת ספרייה, ורצף על Novaseq 6000.

פרוטוקול

ניתן להפעיל את Rup כקובץ Script יחיד עם מינימום קבצי קלט הנדרשים. הצינור, כפי שמוצג כאן, דורש נתוני ריצוף RNA מזווג של Illumina. צינור מותאם לריצוף קצה יחיד עם אותם שלבי ניתוח מופקד גם במאגר GitHub. יש צורך רק ברצף הגנום כקובץ fasta, מודל הגנים והערות rRNA כקבצי gtf, וקריאות הריצוף הגולמיות כקבצי fastq.gz. ודא שקבצי הגנום והביאורים מסופקים כ-genome.fa, annotation.gtf ו-rRNA.gtf בתיקיה המצוינת במשתנה reference_folder. כאשר כל קבצי הרצף כקבצי .fq.gz נשמרים במשתנה read_file_folder, ניתן להפעיל את Rup באופן הבא:

1. הכנות:

1. התקן את כל חבילות ה-R הדרושות באמצעות Bioconductor.

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

2. הגדר את כל הקבצים הדרושים.
   1. צור תיקיית מקור שתכיל את כל הקבצים הדרושים. הוסף את רצף הגנום הייחוס בקובץ fasta שנקרא "reference/genome.fa" לתיקיית המקור. ספק את ביאור מודל הגן כקובץ gtf: "reference/annotation.gtf", לחלופין, ספק את ביאור הגן rRNA כקובץ gtf: "reference/rRNA.gtf".
   2. הוסף את כל קריאות הרצף כקבצי fastq לתיקיית "קריאות". ודא שכל שמות הקבצים של fastq עוקבים אחר אותה תבנית; _1.fastq.gz ו-_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. הגדר את הפרמטר בהתאם למשאבי המחשוב.

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

4. הגדר את הפרמטר בהתאם לשיטת הרצף.
    

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

   
   
   הערה: אורך הקריאה המינימלי לא צריך להיות קטן מ-25 bp, מכיוון שקריאות קטנות יותר עלולות לגרום לשבירת המיפוי. ערכים גדולים יותר יביאו לקריאות ממופות באופן ייחודי יותר, אך גם ליותר קריאות שנמחקו במהלך החיתוך, בהתאם לאיכות הרצף. פרטים נוספים על פרמטרים אלה מתוארים בתוספים.

2. הערכת איכות ריצוף

הערה: שלב זה ייצור עלילת עמודות המציגה את מספר הקריאות לפני ואחרי חיתוך כל דגימה.

1. הפעל את fastqc⁹, שהוא כלי לבקרת איכות כללית של נתוני רצף בתפוקה גבוהה, המספק מידע על המספר הכולל של קריאות רצף, אורך הקריאה, תוכן GC ממוצע, זיהום מתאם רצף ורצפים מיוצגים יתר (למשל, קריאות rRNA). נתח את איכות הרצף לפי דוחות איכות בסיסיים ולפי רצף. הפעל fastqc על קבצי הרצף הגולמיים באמצעות חבילת R fastqcr¹⁰ ומקם את כל קבצי הקלט באותה ספרייה (read_file_folder). סכם את קבצי הפלט מ-fastqc_folder לאובייקט R באמצעות הפונקציה qc_aggregate.
   הערה: התוצאות יישמרו בספרייה fastqc_folder. ניתן לגשת לקבצי .html שנוצרו על ידי fastqc בתיקיית הפלט לקבלת נתונים סטטיסטיים נוספים.
    

   # 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. כאשר כל קידומות הדגימה מאוחסנות sample_prefixes הווקטור, חזור על הקידומת והפעל את כלי החיתוך האיכותי על כל הדגימות. חתוך כדי להסיר רצפי מתאמים, פריימרים, בסיסים באיכות נמוכה, רצפי poly-A/T ומבוצע על ידי Fastp¹¹, הזמין בחבילת R Rfastp¹². קבצי ה- fastq החתוכים יאוחסנו בתיקיה 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. הפעל את fastqc שוב על הקריאות החתוכות.
    

   # 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. הערך את החיתוך על ידי איסוף מספר הקריאות לפני ואחרי החיתוך עבור כל הדגימות. צור תרשים עמודות עם מספרי הקריאה באמצעות 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)

איור 2: תוצאות ריצוף וחיתוך. קרא ספירות לפני (אדום) ואחרי חיתוך (ירוק). ל-Sample s2_r1 כבר יש ספירת קריאה נמוכה לפני החיתוך, בעוד שהחיתוך הסיר חלק גדול מקריאות s2_r2 הדגימה. כל הדוגמאות האחרות מראות אובדן קריאות מקובל מהחיתוך. אנא לחץ כאן לצפייה בגרסה גדולה יותר של איור זה.

3. איכות מיפוי

הערה: שלב זה מחשב עלילות עמודות כדי להבחין בין קריאות ממופות בודדות (למשל, תעתיקים מגנים מקודדי חלבון) מקריאות מרובות ממופות (למשל, קריאות rRNA) וקריאות לא ממופות (למשל, זיהומים).

1. השתמש בחבילת Rsubread¹⁴ כדי למפות את הקריאות לגנום הייחוס. ראשית, בנה אינדקס של קובץ ה-fasta של גנום הייחוס (genome_fasta_file).
   הערה: שלב זה מבוצע פעם אחת עבור כל גנום ייחוס.
    

   # 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. חזור על כל הדגימות והתאם את הקריאות לגנום הייחוס באמצעות פונקציית align() של חבילת Rsubread. פונקציה זו יוצרת קבצי .bam בתיקיית הפלט.
    

   # 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. ספור את הקריאות עבור כל גן באמצעות הפונקציה featureCounts() מחבילת Rsubread. ודאו שקובץ הביאור (annotation_file) הוא בפורמט GTF. הרוזן קורא רק עם התאמה אחת לגנום.
    

   # 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. ספרו את מיפוי הקריאות לגנים rRNA להערכת תכולת ה-rRNA בדגימות. אפשר לספור קריאות מרובות ממופות כאן. ציין את מיקומי הגן rRNA בקובץ 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. אסוף את סטטיסטיקת הקצאת הקריאה שנוצרה על-ידי featureCounts. הסטטיסטיקה כוללת את מספר היישור עבור כל מדגם בקטגוריות השונות (למשל, Assigned, UnMapped, 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. אסוף את הנתונים הסטטיסטיים של הקצאת ה-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. התווה את סטטיסטיקת מיפוי הקריאה כתרשים עמודות.
    

   # 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. סווג גנים לקבוצות לפי מספר הקריאות שהוקצו להם והתווה את התוצאות כעלילת עמודות.
    

   # 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: קרא את סקירת המיפוי והסטטיסטיקה. רוב הדגימות מראות מספר גבוה של קריאות שהוקצו (צד שמאל, חום) ומספר נמוך של קריאות rRNA (צד ימין, כחול). ל-Sample s2_r3 יש כמות גבוהה במיוחד של קריאות מרובות (צד שמאל, ירוק) המתאימות למספר קריאת rRNA גבוה (צד ימין). s2_r4 הדגימה מראה מספר גבוה של קריאות לא ממופות בשילוב עם מספר קריאת rRNA צפוי, מה שמרמז על זיהום בקריאות מאורגניזם אחר. אנא לחץ כאן לצפייה בגרסה גדולה יותר של איור זה.

4. שכפל איכות.

1. ודא כי שכפולים מאותה דגימה או רקמה מראים מתאם גבוה יותר זה עם זה מאשר עם שכפולים מדגימות או רקמות שונות. שרטט מפת חום מקובצת של המתאמים המשוכפלים כדי להמחיש את הנתונים באמצעות ספירת הקריאה הגולמית או הספירה המנורמלת של תמלילים למיליון (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")

איור 4: סיווג ספירת קריאות גנים. כל הגנים בדגימה מסווגים לאחת מחמש קטגוריות על סמך מספר הקריאות שהוקצו. מוצגים באדום מספר הגנים ללא הקצאת קריאות. לדגימות עם מספר נמוך של קריאות מוקצות (s2_r1 עד s1_r4) יש מספר גבוה יותר של גנים עם 10 עד 100 קריאות מוקצות ומספר נמוך יותר של גנים עם יותר מ-1000 קריאות. גנים עם ביטוי נמוך עלולים להחמיץ בדגימות אלה. אנא לחץ כאן לצפייה בגרסה גדולה יותר של איור זה.

תוצאות

הצינור מיושם לחלוטין כסקריפט R ונבדק במערכות הפעלה לינוקס ו-Mac OS. משתמשי Windows יכולים להשתמש במערכת המשנה של Windows עבור Linux (WSL). הקוד ונתוני הבדיקה זמינים כמאגר GitHub: https://github.com/oliverrupp/rup. נתוני הריצוף זמינים תחת פרויקט ENA EBI PRJEB96400.

עשר דגימות נוצרו באופן מלאכותי משתי דגימות אמיתיות כדי להדגים בעיות מגוונות שעלולות להיתקל במהלך בקרת איכות של RNA-seq בתפזורת באמצעות Rup לניתוח. s2_r1 המדגם תוכנן להציג מספר קריאה כולל נמוך, s2_r2 המדגם מכיל חלק גדול של קריאות באיכות נמוכה שיש להשליך בתהליך החיתוך. s2_r3 הדגימה כוללת חלק גדול מקריאות ה-rRNA והדגימה s2_r4 כוללת קריאות מזהמים שלא הצליחו למפות את גנום הייחוס. שמות הדגימות s1_r5 ו-s2_r5 הוחלפו כדי להמחיש מתאם שכפול נמוך.

סעיף 2 של פרוטוקול מאפשר זיהוי של דגימות עם מספרי קריאה נמוכים לפני או אחרי החיתוך (איור 2). תרשים העמודות מציג את מספר הקריאה הנמוך יותר s1_r1 מדגם לפני ואחרי החיתוך. כאן, מספר הקריאה הראשוני היה נמוך. מספר הקריאה הנמוך של s1_r2 הדגימה לאחר החיתוך מרמז על כמות גדולה של רצפי מתאמים, שגיאות רצף, רצפי פריימר, רצפי פולי-A/T שהוסרו במהלך תהליך החיתוך. ירידה ב-RNA הקלט עשויה גם להפחית את מספר הקריאות באיכות גבוהה.

פרוטוקול סעיף 3 מזהה בעיות בהקצאות קריאה. איור 3 מתאר את המספר המוגבר של קריאות מרובות מיפוי במדגם s2_r3 (ירוק), כמו גם את המספר הגבוה של קריאות rRNA. זיהום של s2_r4 הדגימה ניכר על ידי החלק הגדול של הקריאות שלא מופו לגנום הייחוס (ורוד). קריאות אלה אינן מתואמות ל-rRNA אלא לרצפים מאורגניזם שאינו מטרה. איור 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 שעברו טרנספורמציה של log₂ ומציגה שני אשכולות נפרדים של חמש דגימות כל אחד. שכפולים ביולוגיים/טכניים צפויים להראות מתאם גבוה יותר זה לזה מאשר אלה מרקמות/טיפולים אחרים. האשכול השמאלי מכיל ארבעה שכפולים של דגימה 1 ושכפול אחד של דגימה 2 (s2_r5), האשכול הימני מכיל ארבעה שכפולים של דגימה 2 ושכפול אחד של דגימה 1 (s1_r5). יש לבדוק את השכפולים הבודדים עבור החלפת דגימה אפשרית או השפעות אצווה כדי להסביר את האשכולות שלהם במפת החום. אנא לחץ כאן לצפייה בגרסה גדולה יותר של איור זה.

טבלה 1: השוואה בין צינורות הערכת איכות RNA-seq. לניתוח מפורט, ראה קובץ משלים 1. אנא לחץ כאן להורדת טבלה זו.

קובץ משלים 1: בחירת פרמטרים והפניות. הניתוחים מראים את ההשפעה של בחירת פרמטר והתייחסות על תוצאות ה-QC הכוללות. אנא לחץ כאן להורדת קובץ זה.

דיון

איכות ניתוחי ביטוי הגנים הדיפרנציאלי תלויה מאוד בשני גורמים: מספר הקריאות שרצפו בכל דגימה ושכפול¹⁶ ומספר השכפולים לדגימה^17,18. כאן, אנו מציגים את צינור ה-Rup הידידותי למשתמש כדי לקבוע את מספר הקריאות המתאימות לכימות ביטוי גנים בכל שכפול. מדדים שונים מאפשרים לחוקרים להבין מדוע שכפולים מראים מספרי קריאה נמוכים ולזהות בעיות במתאם הדגימה. למרות ש-Rup פותח להערכת איכות של דגימות RNA-seq של צמחים, הוא מתאים באותה מידה לאורגניזמים איקריוטיים אחרים, ואינו דורש התאמות נוספות עבור דגימות שאינן צמחיות (קובץ משלים 1).

השלב הראשון של Rup קובע את ספירת הקריאה הכוללת לפני ואחרי חיתוך הקריאה. מספר הקריאות המרוצפות קובע את מספר הגנים הניתנים לזיהוי, ואם מספר הקריאה נמוך מדי, גנים רבים המתבטאים באופן דיפרנציאלי יישארו בלתי מזוהים. חלק גדול מהקריאות שנזרקות במהלך תהליך החיתוך והסינון האיכותי עשויות להצביע על פירוק RNA, מה שעלול להוביל להערכת חסר של ביטוי לגנים פגומים מאוד. עם זאת, האם ניתן להשתמש בדגימה ליישומים במורד הזרם תלויה באורגניזם היעד ובמטרת המחקר. לדוגמה, כדי ללכוד את הביטוי של רוב הגנים הצמחיים, מניסיוננו, נדרשים 30 עד 50 מיליון קריאות, בעוד שעבור פטריות, רק 10 מיליון עשויים להספיק. לפיכך, Rup לא יגדיר סף איכות להחרגת מדגמים בעייתיים אלא יספק מדדים לזיהוי בעיות מגוונות שעלולות להיתקל בהן.

השלב השני (איכות המיפוי) מעריך את הדיוק והאמינות של יישור הקריאה לגנום הייחוס, ומפרט את התאמתם לניתוחים במורד הזרם. לא כל הקריאות הרצופות יכולות לשמש לחישוב שפע גנים; מתעלמים ברוב המקרים מקריאות שאינן מתיישבות עם הגנום או הטרנסקריפטום או קוראות שממפה למספר מיקומים בגנום¹⁹ ואינן תורמות לספירת הקריאה הכוללת. ניתן להשתמש בתוצאות השלב השני של הצינור כדי להבין מדוע לא נעשה שימוש בקריאות בחישוב השפע. מספר גבוה של קריאות לא ממופות יכול להצביע על זיהום במהלך מיצוי RNA (למשל, עם פתוגנים צמחיים או אוכלי עשב) או גנום ייחוס לא שלם. מודלים גנטיים לא שלמים של גנום הייחוס עלולים לגרום למספר גבוה של קריאות "ללא תכונה". חלק גדול מהקריאות הממופות המרובות יכול להיגרם על ידי מספר גבוה של גנים rRNA בספרייה המצביע על הסרת rRNA לא מספקת במהלך תהליך הכנת ספריית הריצוף (במקרה שהקריאות המרובות נגזרות באמת מ-rRNAs, ניתן לנתח אותן עוד יותר על ידי מתן קובץ הערות rRNA לצינור, אשר לאחר מכן מחשב אוטומטית את מספר קריאות ה-rRNA האפשריות).

מדד האיכות השלישי, החשוב לא פחות, הוא המתאם בין שכפולים של אותו מדגם. באופן כללי, המתאם בין שכפולים של אותה דגימה צריך להיות גבוה יותר מהמתאם בין שכפולים של דגימות שונות. מתאם נמוך בין שכפולים יכול להצביע על שונות ביולוגית גדולה בין השכפולים, דמיון גבוה בין הדגימות/תנאים/רקמות, השפעות אצווה אחרות, או אפילו החלפת דגימות או תיוג שגוי. Rup מחשב מתאמים זוגיים בין כל הדגימות ומייצר מפת חום מקובצת של מתאמי הדגימה. בנוסף, ניתן לחשב ניתוח רכיבים עיקריים (PCA) על הדגימות כדי לזהות השפעות אצווה אפשריות שיש להכיר בהן בניתוחים נוספים בהמשך הזרם. למרות שניתן לתקן את השפעות האצווה בניתוחים במורד הזרם, ייתכן שעדיף להסיר דגימות בעייתיות אם ההבדל הנמדד גבוה מדי.

לחומר קלט יש השפעה ישירה על איכות הריצוף: דגימת רקמות, מיצוי RNA והכנת ספרייה הם צעדים קריטיים למזעור בעיות איכות RNA-seq. יש לשמור על תנאים עקביים, במידת האפשר, למשל, באמצעות תאי גידול. יש לבצע אמצעי הדברה בזמן, שכן יש לבחור רק אנשים בריאים לדגימה. כמו כן, מומלץ לאסוף דגימות באותו יום ו/או באותו זמן כדי להפחית את השונות בשעון הביולוגי. קיימות ערכות רבות למיצוי RNA, ובחירת ערכה מתאימה למין היעד יכולה לשפר את איכות ה-mRNA. פרוטוקולי הכנת הספרייה צריכים לכלול שלבים להעשרת polyA+ RNA כדי למזער את מקטע ה-rRNA.

Rup יכול לשמש כשלב בקרת איכות ראשוני בניתוח ביטוי גנים דיפרנציאלי. בקרת איכות זו נדרשת ממספר סיבות קריטיות, מכיוון שהיא מבטיחה את איכות נתוני הקלט (מזהה שכפולים/דגימות באיכות נמוכה, יכולה להגדיר ספי איכות עבור עומק קריאה ושיעורי מיפוי) ומזהה בעיות טכניות כגון שגיאות רצף, השפעות אצווה או תיוג שגוי. אם כניסת ה-RNA איכותית מספיק, Rup יכולה לסייע על ידי חיתוך קריאות באיכות נמוכה או זיהוי דגימות עם תווית שגויה. עם זאת, Rup אינו יכול לפצות על איכות RNA קלט ירודה. במקרים של RNA באיכות נמוכה או נתוני ריצוף שגויים, ייתכן שיהיה צורך לאסוף מחדש דגימות, להתאים את פרוטוקול מיצוי ה-RNA ו/או לחזור על ריצוף. כך גם לגבי דגימות עם חלק גדול מהקריאות הלא ממופות שייתכן שנגרמו מזיהום. עם זאת, מניסיוננו, לא ניתן פשוט לחזור על כל מיצוי ה-RNA בגלל חוסר זמינות חומר. במקרה זה, הם אינם מתאימים ליישומים מסוימים במורד הזרם: אין לנתח דגימות עם מספר נמוך של קריאות ממופות בודדות בניתוחי ביטוי גנים דיפרנציאליים, אך הן עדיין מחזיקות מידע לניתוח נוכחות תעתיק. במקרה זה, היעדר תעתיקים ברקמות/טיפולים/מצבים לא יכול להילקח בחשבון לצורך ניתוח וכימות של שפע התעתיקים.

Rup כולל כמה מגבלות. לדוגמה, הוא אינו בודק ישירות את פירוק ה-RNA מכיוון שנדרשות מדידות ישירות של שלמות ה-RNA לפני הכנת הספרייה והרצף. עם זאת, סקריפט לזיהוי ירידה ב-RNA באמצעות מודולי RSeQC geneBodyCoverage.py ו-tin.py כלול במאגר של צינור זה. יתר על כן, Rup אינו בודק הטיית תוכן GC ואורך תמלול. גודל גנום הייחוס מוגבל כיום ל-4 ג'יגה-בייט, ועל פי הניסיון שלנו, מודול מיפוי הקריאה מעריך יתר על המידה קריאות מרובות ממופות בפוליפלואידים, כפי שמודגם בקובץ המשלים בהשוואה של מיפוי קריאה לגרסאות גנום הפלואיד לעומת גרסאות הגנום הדיפלואידי של E. californica, גנום הפלואידי הוא אפוא הקלט המועדף לצינור זה. איכות הפלט של Rup תלויה במידה רבה באיכות גנום הייחוס וההערות, כך שרצף גנום לא שלם או מקוטע מאוד עלול להוביל להערכת יתר של קריאות לא ממופות. יתר על כן, ביאור מודל גנים לא שלם מוביל להערכת יתר של קריאות לא מוקציות. ניתן להסיק על שלמות ושכפול של רצף הגנום והערת הגנים באמצעות כלים כמו BUSCO²⁰.

Rup הוא כלי עצמאי לכל שלבי בקרת האיכות הראשוניים החיוניים בניסויי RNA-seq. בניגוד ל-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 Gb, זה חובה, מכיוון ש-Rsubread aligner מוגבל לגנומים קטנים מ-4 Gb. הערת rRNA נעשית עם barrnap, המזהה את גני ה-rRNA השמורים ביותר. על פי הניסיון שלנו, ביאור גן rRNA אינו דורש מקיף, מכיוון שגם אם חסרים כמה גנים חריגים של rRNA, סקירה גסה של נוכחות rRNA במערך הנתונים מספיקה להערכת איכות מערך הנתונים. גרסאות עתידיות של Rup עשויות לעבור לשיטת מיפוי אחרת, כמו כלי היישור המדומה סלמון²¹, כדי להפחית את זמן הריצה. יתר על כן, ניתן להוסיף ל-Rup שיטות נורמליזציה ותיקון נוספות, כמו הטיית GC או תיקון הטיית אורך.

לסיכום, Rup מספקת מידע חיוני כדי להבטיח את האמינות והשחזור של ניתוח נתוני RNA-seq. צינור זה מדווח באופן מקיף על מדדי האיכות העיקריים של RNA-seq ומייצר קבצי פלט לשימוש ישיר בניתוחים במורד הזרם. הוא תוכנן ככלי עצמאי לחוקרים עם ידע מינימלי בביואינפורמטיקה כדי להעריך את איכות נתוני הריצוף העיקריים שלהם באמצעות הדמיה אינטואיטיבית.

חומרים

רשימת החומרים שנעשה בהם שימוש במאמר זה

שם	חברה	מספר קטלוג	הערות
fastQCR	R	0.1.3
FUJITSU ESPRIMO D958 Desktop	פוג'יטסו		הצינור פותח ונבדק על Ubuntu Linux 24.02 (32 Gb RAM)
getopt	R	1.20.4
ggplot2	R	3.5.2
MacBook Pro	אפל		הצינור נבדק על MaxOS 15.4.1 (זיכרון 16 גיגה-בייט)
פיטמאפ	R	1.0.13
R		4.4.3
reshape2	R	1.4.4
rfastp	ביו-קונדוקטור	1.16.0
rsamtools	ביו-קונדוקטור	2.22.0
rsubread	ביו-קונדוקטור	2.20.0