Rup (RNA-seq Usability Assessment Pipeline) - Контроль качества для экспериментов по объемному RNA-seq у эукариот

Краткое содержание

Этот протокол позволяет проводить начальный контроль качества экспериментов с РНК-секвенированием для биологов, работающих в мокрых лабораториях с ограниченным опытом в области биоинформатики.

Аннотация

Современные подходы в молекулярной науке о растениях часто требуют экспериментов с объемным RNA-seq, например, для отслеживания глобальных изменений в транскриптомах при обработке или для определения ключевых компонентов регуляторных путей. Следовательно, различные области науки о растениях полагаются на высококачественные и воспроизводимые объемные данные РНК-секвенирования для научного прогресса. Однако, по нашему опыту, часто отсутствуют знания и применение мер контроля качества в наборах данных RNA-seq. Здесь мы представляем Rup (конвейер оценки юзабилити РНК-секвенирования) для контроля качества объемных данных РНК-секвенирования для последующего анализа экспрессии генов, который является автономным и легко применимым для биологов, работающих в мокрых лабораториях с базовыми знаниями R. Rup помогает различать данные секвенирования высокого качества, пригодные для последующих экспериментов по экспрессии генов, и те, которые не подходят для общего дальнейшего анализа. Rup включает в себя тесты для нескольких часто встречающихся проблем, таких как недостаточное количество считываемых или картографирование, идентификация загрязнений, количественная оценка фракций рРНК в общих данных RNA-seq, тестирование сходства реплик, использование реальных данных для демонстрации и интуитивная визуализация. Rup предоставляет набор инструментов для выявления экспериментальных недостатков перед стандартизированным транскриптомным анализом, тем самым повышая качество данных для отдельных исследователей и в этой области. Это повышает доверие к объемному анализу данных РНК-секвенирования и обеспечивает основу для будущих руководящих принципов, определяющих минимальные критерии контроля качества, тем самым повышая надежность и прозрачность публикуемых данных РНК-секвенирования. Данные о разрыве и испытаниях доступны по адресу https://github.com/oliverrupp/rup.

Введение

Эксперименты по транскриптомике (РНК-секвенирование) всесторонне исследуют транскрипционные сигнатуры, формирующие фенотипы. Этот подход стал незаменимым в молекулярной генетике растений для идентификации одиночных или совместно экспрессируемых генов и биологических процессов, участвующих в развитии, взаимодействии растительных патогенов или устойчивости к абиотическому стрессу^{1, 2 ,3}. Последние достижения в технологии RNA-seq повысили специфичность и позволили обнаруживать различные изоформы и варианты с разрешением одного основания, что позволяет идентифицировать вариации последовательностей от более крупных инделов до однонуклеотидных полиморфизмов (SNP). Данные, полученные с помощью РНК-секвенирования, характеризуются широким динамическим диапазоном, облегчающим обнаружение как высокораспространенных, так и низкоэкспрессированных транскриптов, и требуют соответствующих экспериментальных условий для согласованности. Кроме того, наборы данных RNA-seq имеют большие размеры и, следовательно, требуют больших вычислительных ресурсов для анализа и хранения, что влечет за собой эффективное управление данными и обширные вычислительные ресурсы. РНК-секвенирование требует высококачественного ввода данных на каждом этапе рабочего процесса, так как слабые места на любом этапе могут распространяться и ставить под угрозу результаты. Плохая целостность РНК или технические проблемы при подготовке библиотеки приведут к систематической ошибке и снижению точности. Параметры для получения высококачественных необработанных прочтений могут различаться в зависимости от средств секвенирования нового поколения (NGS).

По нашему опыту, мы рекомендуем использовать только РНК с числом целостности РНК (RIN) выше 7, что свидетельствует о практически неповрежденной структуре мРНК, в качестве исходных данных для подготовки библиотеки РНК-секвенирования. Для успешного секвенирования требуется примерно 2 мкг общей РНК в концентрации 50–200 нг/мкл для стандартных протоколов подготовки библиотек. Чистота РНК должна быть подтверждена соотношением OD_260/280 от 1,8 до 2,1 и отношением OD_260/230 больше 1,5 с помощью спектрофотометра. Недостаточная глубина секвенирования, низкая скорость картирования или смещение с референсным геномом могут еще больше исказить экспрессию генов и количественную оценку сплайсинга. Кроме того, неадекватный дизайн эксперимента, такой как отсутствие дискриминации между образцами различных тканей или методов лечения и высокая вариабельность между репликациями, может привести к появлению шума и снижению воспроизводимости. В то время как многие лаборатории регулярно анализируют транскриптомы, строгий контроль качества первых основных этапов анализа часто не сообщается или может полностью отсутствовать в публикациях. Это может привести к чрезмерной интерпретации результатов, полученных в результате транскриптомного анализа, и, как следствие, к невоспроизводимым результатам.

Здесь мы обеспечиваем рабочий процесс для контроля качества начальных этапов, необходимых для высококачественного транскриптомного анализа профилей мРНК для измерения транскрипционных изменений. Наша цель состоит в том, чтобы дать возможность биологам с ограниченными знаниями в области биоинформатики оценить свои первичные данные транскриптомики. Rup (рис. 1) доступен для исследователей, знакомых с базовыми знаниями R. Выполнение представленного здесь рабочего процесса предоставит исследователям подробное понимание их первичных данных, включая их потенциальные ограничения для последующего анализа. Насколько нам известно, до сих пор отсутствует практический конвейер оценки первичных транскриптомных данных в сочетании с рекомендациями по различению данных высокого и низкого качества.

Рисунок 1: Рабочий процесс конвейера контроля качества RNA-seq in silico . Входные файлы для контроля качества извлекаются из данных РНК-секвенирования, полученных в результате секвенирования растительного материала, а также из общедоступных наборов данных. Предоставленный конвейер R оценивает качество последовательности с помощью трех основных подходов: качество секвенирования, качество картирования и качество репликации. Вычисляются различные метрики качества, а статистические результаты визуализируются с помощью общих пакетов R, таких как ggplot2 и pheatmap (например, гистограммы и тепловые карты). Пожалуйста, нажмите здесь, чтобы просмотреть увеличенную версию этой цифры.

Ранее заявленные конвейеры оценки требовали предварительно обработанных данных (например, выравнивание чтения в виде файлов .bam), не охватывают все метрики или больше не поддерживаются ^4,5,6. Преимущество представленного здесь рабочего процесса заключается в его комплексности за счет объединения наиболее распространенных вопросов контроля качества в единый конвейер. Кроме того, мы приводим примеры высококачественных данных, подходящих для всех последующих приложений анализа, а также примеры данных низкого качества, а также обсуждаем их конкретные ограничения для дальнейшего анализа. В нескольких ранее опубликованных отчетах описывается назначение инструментов анализа РНК-секвенирования и даются сравнительные оценки их эффективности ^7,8. Тем не менее, стандарты контроля качества РНК-секвенирования и методологических отчетов не стандартизированы и ухудшают воспроизводимость и биологически значимые интерпретации транскриптомных экспериментов.

Rup интегрирует стандартные высококачественные инструменты, анализ контроля качества и визуализацию результатов. Rup требует необработанных прочтений секвенирования и аннотированного генома в качестве входных данных и работает в Mac OS, Linux и «Windows Subsystem for Linux» (WSL) в системах Windows. Он интегрирует инструменты для обеспечения качества секвенирования и количественной оценки сопоставления прочтений с геномом, включает измерение содержания рРНК в образцах и предлагает корреляцию образцов для оценки корреляции репликации. Тестовые данные были получены с помощью лазерной микродиссекции центральной меристемной ткани двух разных стадий растения вида Eschscholzia californica, протокола со сверхнизким уровнем ввода для подготовки библиотеки, и секвенированы на Novaseq 6000.

Протокол

Rup может быть запущен как один скрипт с минимальным количеством входных файлов. Для конвейера, как показано здесь, требуются данные секвенирования парных концов РНК Illumina. Настроенный конвейер для одностороннего секвенирования с теми же шагами анализа также помещается в репозиторий GitHub. Требуется только последовательность генома в виде файла fasta, модель гена и аннотации рРНК в виде файлов 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 в папку "reads". Убедитесь, что все имена файлов 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.н., так как меньшие объемы чтения могут привести к нарушению отображения. Более высокие значения приведут к большему количеству уникально сопоставленных операций чтения, а также к большему количеству отброшенных операций чтения во время обрезки, в зависимости от качества секвенирования. Более подробно об этих параметрах рассказывается в приложениях.

2. Оценка качества секвенирования

ПРИМЕЧАНИЕ: На этом шаге будет создана столбчатая диаграмма, показывающая количество прочтений до и после обрезки каждого образца.

1. Запустите fastqc⁹, который представляет собой инструмент для общего контроля качества данных последовательностей с высокой пропускной способностью, предоставляющий информацию об общем количестве секвенированных прочтений, длине чтения, среднем содержании GC, загрязнении адаптера секвенирования и чрезмерно представленных последовательностях (например, чтениях рРНК). Анализируйте качество секвенирования в соответствии с отчетами о качестве базы и последовательности. Запустите 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, выполните перебор префиксов и запустите инструмент обрезки качества fastp для всех образцов. Обрезка для удаления переходных последовательностей, праймеров, некачественных баз, поли-А/Т последовательностей и выполняется с помощью 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: Результаты секвенирования и обрезки. Счетчики чтения до (красный) и после обрезки (зеленый). Образец s2_r1 уже имеет низкое количество прочтений до обрезки, в то время как обрезка удалила большую часть s2_r2 чтения образца. Все остальные образцы показывают приемлемые потери чтения от обрезки. Пожалуйста, нажмите здесь, чтобы просмотреть увеличенную версию этой цифры.

3. Качество картографирования

ПРИМЕЧАНИЕ: На этом шаге рассчитываются столбчатые диаграммы для различения одиночных картированных прочтений (например, транскриптов генов, кодирующих белки) от многокартированных прочтений (например, рРНК) и некартированных прочтений (например, контаминаций).

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. Count считывает только одно совпадение с геномом.
    

   # 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. Подсчитайте сопоставление прочтений с генами рРНК для оценки содержания рРНК в образцах. Позвольте здесь подсчитывать многосопоставленные чтения. Укажите локусы гена рРНК в файле 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. Статистика включает количество выравниваний для каждой выборки в различных категориях (например, Назначено, Несопоставлено, Мультисопоставлено).
    

   # 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. Сбор статистики присвоения рРНК.
    

   # 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: Чтение обзора карт и статистики. Большинство образцов показывают большое количество назначенных прочтений (левая сторона, коричневый) и низкое количество прочтений рРНК (правая сторона, синий). Образец s2_r3 имеет необычно большое количество мультикартированных прочтений (левая сторона, зеленый), что соответствует высокому числу прочтений рРНК (правая сторона). Образец s2_r4 показывает большое количество некартированных прочтений в сочетании с ожидаемым числом прочтений рРНК, что указывает на загрязнение чтениями из другого организма. Пожалуйста, нажмите здесь, чтобы просмотреть увеличенную версию этой цифры.

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 и был протестирован на операционных системах Linux и Mac OS. Пользователи Windows могут использовать подсистему Windows для Linux (WSL). Код и тестовые данные доступны в виде репозитория GitHub: https://github.com/oliverrupp/rup. Данные секвенирования доступны в рамках проекта ENA EBI PRJEB96400.

Десять образцов были созданы искусственно из двух реальных образцов, чтобы проиллюстрировать различные проблемы, с которыми можно столкнуться во время контроля качества объемного RNA-seq с использованием Rup для анализа. Образец s2_r1 был разработан для демонстрации низкого общего числа прочтений, образец s2_r2 содержит большую долю некачественных прочтений, которые должны быть отброшены в процессе обрезки. Образец s2_r3 включает большую часть прочтений рРНК, а образец s2_r4 включает загрязняющие чтения, которые не смогли сопоставить с референсным геномом. Названия образцов s1_r5 и s2_r5 были поменяны местами, чтобы проиллюстрировать низкую корреляцию реплик.

Раздел 2 протокола позволяет идентифицировать образцы с низким числом считываемых до или после обрезки (рисунок 2). Линейчатая диаграмма показывает меньшее число считываний в выборке s1_r1 как до, так и после обрезки. Здесь первоначальное число прочтений было низким. Низкое число считываемых образцов s1_r2 после обрезки свидетельствует о большом количестве последовательностей адаптеров, ошибок секвенирования, последовательностей праймеров, секвенированных поли-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 log₂ и показывает два отдельных кластера по пять выборок в каждом. Ожидается, что биологические/технические реплики продемонстрируют более высокую корреляцию друг с другом, чем реплики из других тканей/методов лечения. Левый кластер содержит четыре реплики образца 1 и одну реплику образца 2 (s2_r5), правый кластер содержит четыре реплики образца 2 и одну реплику образца 1 (s1_r5). Отдельные реплики должны быть проверены на возможную подмену образцов или эффекты пакетной обработки, чтобы объяснить их кластеризацию на тепловой карте. Пожалуйста, нажмите здесь, чтобы просмотреть увеличенную версию этой цифры.

Таблица 1: Сравнение конвейеров оценки качества RNA-seq. Подробный анализ см. в Дополнительном файле 1. Пожалуйста, нажмите здесь, чтобы скачать эту таблицу.

Дополнительный файл 1: Выбор параметров и эталонов. Анализ показывает влияние выбора параметров и референса на общие результаты контроля качества. Пожалуйста, нажмите здесь, чтобы загрузить этот файл.

Обсуждение

Качество анализа дифференциальной экспрессии генов в значительной степени зависит от двух факторов: количества прочтений, секвенированных в каждом образце и реплицитных¹⁶, и количества повторений на образце ^17,18. В этой статье мы представляем удобный для пользователя конвейер Rup для определения количества прочтений, пригодных для количественной оценки экспрессии генов в каждой реплике. Различные метрики позволяют исследователям понять, почему реплики показывают низкие присвоенные номера прочтений, и выявить проблемы в корреляции выборки. Несмотря на то, что Rup был разработан для оценки качества образцов РНК-секвенирования растений, он в равной степени подходит и для других эукариотических организмов, не требуя дальнейших корректировок для образцов, не относящихся к растениям (Дополнительный файл 1).

На первом шаге Rup определяется общее количество прочтений до и после обрезки чтения. Количество секвенированных прочтений определяет количество обнаруживаемых генов, и если число считываемых слишком мало, многие дифференциально экспрессируемые гены останутся необнаруженными. Большая часть прочтений, отброшенных в процессе качественной обрезки и фильтрации, может указывать на деградацию РНК, что может привести к недооценке экспрессии сильно деградированных генов. Однако возможность использования образца для последующих применений зависит от целевого организма и цели исследования. Например, чтобы зафиксировать экспрессию большинства генов растений, по нашему опыту, требуется от 30 до 50 миллионов прочтений, в то время как для грибов может быть достаточно только 10 миллионов. Таким образом, Rup не будет определять пороговые значения качества для исключения проблемных выборок, а скорее предоставит метрики для выявления различных проблем, с которыми можно столкнуться.

На втором этапе (качество картирования) оценивается точность и надежность выравнивания прочтений по референсному геному, а также указывается их пригодность для последующих анализов. Не все секвенированные чтения могут быть использованы для вычисления численности генов; Прочтения, которые не совпадают с геномом или транскриптомом, или чтения, которые сопоставляются с несколькими местами в геноме, в^{большинстве случаев} игнорируются и не влияют на общее количество прочтений. Результаты второго шага конвейера могут быть использованы для понимания того, почему чтения не используются при вычислении численности. Большое количество некартированных прочтений может указывать на загрязнение во время экстракции РНК (например, растительными патогенами или травоядными) или на неполный референсный геном. Неполные модели генов референсного генома могут привести к большому количеству прочтений «без признаков». Большая часть мультикартированных прочтений может быть вызвана большим количеством генов рРНК в библиотеке, что указывает на недостаточное удаление рРНК в процессе подготовки библиотеки секвенирования (в случае, если мультикартированные чтения действительно получены из рРНК, их можно дополнительно проанализировать путем предоставления файла аннотации рРНК конвейеру, который затем автоматически вычисляет количество возможных прочтений рРНК).

Третья, не менее важная метрика качества — корреляция между тиражами одной и той же выборки. Как правило, корреляция между репликами одной и той же выборки должна быть выше, чем корреляция между репликами разных образцов. Низкая корреляция между репликатами может указывать на большую биологическую изменчивость между репликатами, высокую схожесть между образцами/условиями/тканями, другие эффекты партии или даже на замену образцов или неправильную маркировку. Rup вычисляет попарные корреляции между всеми выборками и создает кластеризованную тепловую карту корреляций выборки. Кроме того, на образцах может быть проведен анализ главных компонент (PCA) для выявления возможных эффектов партии, которые должны быть учтены в дальнейших анализах по технологической цепочке. Несмотря на то, что эффекты партии могут быть скорректированы при последующих анализах, возможно, лучше удалить проблемные образцы, если измеренная разница слишком велика.

Входной материал оказывает непосредственное влияние на качество секвенирования: забор образцов ткани, экстракция РНК и подготовка библиотеки являются критически важными этапами для минимизации проблем с качеством РНК-секвенирования. Следует поддерживать постоянные условия, по возможности, например, с помощью камер выращивания. Мероприятия по борьбе с вредителями должны проводиться своевременно, так как для отбора проб следует отбирать только здоровых особей. Кроме того, рекомендуется собирать образцы в тот же день и/или в одно и то же время, чтобы уменьшить вариацию циркадного транскриптома. Существует множество наборов для экстракции РНК, и выбор подходящего набора для целевого вида может улучшить качество мРНК. Протоколы подготовки библиотеки должны включать этапы обогащения polyA+ РНК для минимизации фракции рРНК.

Rup может быть использован в качестве начального этапа контроля качества при дифференциальном анализе экспрессии генов. Этот контроль качества необходим по нескольким важным причинам, поскольку он обеспечивает качество входных данных (выявляет реплики/образцы низкого качества, может устанавливать пороговые значения качества для глубины чтения и скорости отображения) и выявляет технические проблемы, такие как ошибки секвенирования, эффекты пакетной обработки или неправильная маркировка. Если вход РНК имеет достаточное качество, Rup может помочь, обрезая некачественные чтения или выявляя неправильно меченые образцы. Однако Rup не может компенсировать плохое качество входной РНК. В случаях низкого качества РНК или ошибочных данных секвенирования может потребоваться повторный сбор образцов, корректировка протокола экстракции РНК и/или повторное секвенирование. То же самое относится и к образцам с большой долей некартированных прочтений, которые могли быть вызваны загрязнением. Однако, по нашему опыту, не все экстракции РНК можно просто повторить из-за нехватки материала. В этом случае они не подходят для некоторых последующих применений: образцы с небольшим числом одиночных картированных прочтений не должны анализироваться при дифференциальном анализе экспрессии генов, но они все еще содержат информацию для анализа наличия транскриптов. В этом случае отсутствие транскриптов в тканях/методах лечения/условиях не может быть рассмотрено для анализа и количественной оценки распространенности транскриптов.

Rup включает в себя некоторые ограничения. Например, он не проверяет деградацию РНК напрямую, поскольку перед подготовкой библиотеки и секвенированием требуются прямые измерения целостности РНК. Тем не менее, скрипт для идентификации деградации РНК с помощью модулей RSeQC geneBodyCoverage.py и tin.py включен в репозиторий этого конвейера. Кроме того, Rup не проверяет на погрешность по содержанию GC и длине транскрипта. Размер референсного генома в настоящее время ограничен 4 Гб, и, согласно нашему опыту, модуль картирования чтения переоценивает мультикартированные чтения в полиплоидах, как показано в дополнительном файле при сравнении картирования чтения с гаплоидными и диплоидными версиями генома E. californica, таким образом, гаплоидный геном является предпочтительным входом для этого конвейера. Качество вывода Rup во многом зависит от качества референсного генома и аннотации, так что неполная или сильно фрагментированная последовательность генома может привести к переоценке некартированных прочтений. Кроме того, неполная аннотация генной модели приводит к завышению количества нераспределенных прочтений. Полнота и дупликация последовательности генома и аннотации генов могут быть выведены с помощью таких инструментов, как BUSCO²⁰.

Rup является самостоятельным инструментом для всех начальных этапов контроля качества, необходимых в экспериментах с РНК-секвенированием. В отличие от RSeQC и RNA-SeQC, предварительная обработка необработанных прочитанных данных секвенирования не требуется. Анализ качества секвенирования не может быть выполнен с помощью RNA-SeQC, а тепловые карты корреляции образцов не рассчитываются с помощью RSeQC и RNA-SeQC (Таблица 1). Кроме того, выход нормализации в RSeQC указан в FPKM, а визуализация вывода не реализована в качестве стандарта для всех модулей в RSeQC и RNA-SeQC. Таким образом, в двух альтернативных инструментах не проверяются некоторые проблемы контроля качества, которые включают в себя низкое число прочтений, высокую долю обрезанных прочтений и идентификацию повторяющихся выбросов. Сравнение Rup, RSeQC и RNA-SeQC доступно в Дополнительном файле 1. Более того, Rup может быть использован в дополнение к RNA-SeQC или RSeQC, например, отсортированные BAM-файлы, полученные этим конвейером, могут быть использованы в качестве входных данных для этих инструментов.

В настоящее время Rup оптимизирован для небольшого количества выборок, но наиболее трудоемкие шаги, такие как качественная обрезка и сопоставление чтения, могут быть предварительно вычислены, например, на компьютерном кластере или облачной инфраструктуре. Для геномов размером более 4 Гб это обязательно, так как выравниватель Rsubread ограничен геномами размером менее 4 Гб. Аннотация рРНК выполняется с помощью barrnap, который идентифицирует высококонсервативные гены рРНК. По нашему опыту, аннотация генов рРНК не требует полноты, так как даже при отсутствии некоторых аномальных генов рРНК для оценки качества датасета достаточно общего обзора наличия рРНК в датасете. Будущие версии Rup могут переключиться на другой метод отображения, такой как инструмент псевдовыравнивания salmon²¹, чтобы сократить время выполнения. Более того, в Rup можно было бы добавить больше методов нормализации и коррекции, таких как смещение GC или коррекция смещения длины.

Таким образом, Rup предоставляет важную информацию для обеспечения надежности и воспроизводимости анализа данных RNA-seq. Этот конвейер всесторонне сообщает об основных показателях качества РНК-секвенирования и создает выходные файлы для непосредственного использования в последующих анализах. Он был разработан как автономный инструмент для исследователей с минимальными знаниями в области биоинформатики для оценки качества данных первичного секвенирования с помощью интуитивно понятной визуализации.

Материалы

Список материалов, использованных в этой статье

Имя	Компания	Каталожный номер	Комментарии
fastqcr	R	0.1.3
FUJITSU ESPRIMO D958 Настольный компьютер	FUJITSU		конвейер был разработан и протестирован на Ubuntu Linux 24.02 (32 Гб оперативной памяти)
getopt	R	1.20.4
ggplot2	R	3.5.2
MacBook Pro	Apple		конвейер был протестирован на maxOS 15.4.1 (16 Гб оперативной памяти)
pheatmap	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