Транскриптомический анализ на основе объёмных данных РНК-секвенции

Summary

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

Abstract

Безалкогольная жировая печёнка (NAFL) обычно считается доброкачественным заболеванием; однако при прогрессировании до неалкогольного стеатогепатита (NASH) пациенты сталкиваются с значительно повышенным риском развития терминальной стадии заболевания печени. Многие исследования пытаются прояснить молекулярный механизм, лежащий в основе перехода от NAFL к NASH. Высокопроизводительные технологии секвенирования (такие как объёмная РНК-секвенция) позволили исследователям глубже понять, изучая транскриптом, выявляя экспрессию молекул, активацию сигнальных путей и другие факторы, связанные с прогрессированием заболевания. Существует множество открытых данных, доступных для исследователей, которые можно проанализировать с целью выявления потенциальных целей для лечения заболеваний. Однако связанные исследования ограничены отсутствием эффективного и надёжного процесса предварительного анализа транскриптома. Здесь предоставляется высоковоспроизводимый и удобный для пользователя анализ в верхней части и последующий конвейер дифференциального анализа генов для обеспечения стандартизированной обработки и глубокого разбора частных или публичных данных. Конвейер делится на четыре этапа: (1) контроль качества данных; (2) картирование генов; (3) дифференциальный генный анализ; и (4) функциональный анализ. Этот процесс направлен на выявление молекулярных механизмов трансформации заболеваний и помощь исследователям в скрининге потенциальных мишеней для лекарств и терапевтических подходов посредством анализа данных Bulk RNA-seq.

Introduction

Неалкогольная жировая болезнь печени (НАЖББ) является самым распространённым хроническим заболеванием печени в мире, затрагивающим более четверти населения. Заболеваемость резко выросла за последние^{десятилетия 1,2,3}. Растущее бремя болезней, особенно её более развитая форма — безалкогольный стеатогепатит (NASH), представляет собой серьёзную глобальную проблему для здоровья и тяжёлую экономическую^{нагрузку 4}. Первая стадия НАЖБП — это безалкогольная жировая печёнка (НАФЛ), сопровождающаяся воспалениями и фиброзом, которые могут прогрессировать в НАСГ. Последний значительно увеличивает риск прогрессирования в терминальную стадию заболевания печени, включая цирроз и гепатоцеллюлярную карциному (ГЦК)^5,6,7. Заболеваемость и смертность от ГЦК связаны с увеличением^{NASH 8,9}, и ожидается, что к 2030 году NAFLD/NASH станет ведущим показателем трансплантации печени^. Однако клиническое прогрессирование НАЖБП крайне^{гетерогенно 11}, что серьёзно затрудняет разработку соответствующих^{препаратов 12}, поэтому особенно важно точно изучать молекулярные механизмы, участвующие в процессе.

Получение информации о клеточном составе на основе массовой РНК-секвенции может значительно прояснить патогенез различных заболеваний. В последние десятилетия было проведено множество исследований с объёмом РНК-секвенации на модельных организмах и людях для выяснения различий в экспрессии генов в прогрессии^NASH 13,14,15 и выявления новых терапевтических целей для вмешательства. Основываясь на массовом анализе РНК-секвенсирующей, Сюн и соавторы обнаружили, что непаренхимальные клетки (NPC) печени участвуют в процессах, таких как формирование внеклеточного матрикса и клеточная адгезия, что способствует прогрессированию^{NASH 16}. Ли и соавторы показали, что белок, ассоциирующий опухоль Вильмса (WTAP) печени (WTAP), регулирует накопление и воспаление эктопических липид, тем самым способствуя формированию^{NASH 17}. Хотя массовый анализ РНК-секвенции является мощным инструментом для разъяснения механизмов NASH, его результаты крайне чувствительны к качеству исходящих данных. Гетерогенность экспериментальных операций и процессов анализа на этапе может серьёзно снизить достоверность данных, тем самым скрывая истинную биологическую информацию и мешая точности последующих анализов. Поэтому важно разработать набор стандартизированных процедур анализа.

По сравнению с секвенированием одноклеточной РНК (scRNA-seq), объёмная РНК-секвенация обладает рядом явных преимуществ как в экспериментальном проектировании, так и в практических применениях. Хотя scRNA-seq позволяет идентифицировать клеточную гетерогенность на уровне одной клетки и даёт точный анализ специфических к типу клеток транскрипционных признаков, он связан с высокой стоимостью, сложными требованиями к обработке данных и ограниченной чувствительностью к обнаружению транскриптов с низким^{количеством 18}. В отличие от этого, объемный РНК-секвенир обеспечивает большую глубину секвенирования, меньшую стоимость и большую пропускную способность выборки, что делает его особенно подходящим для анализа дифференциальной экспрессии генов на уровне популяции и изучения молекулярных^{механизмов 19}. Таким образом, при стандартизированном аналитическом процессе объёмная РНК-секвенция остаётся эффективным, экономически эффективным и надёжным подходом для изучения молекулярной основы сложных заболеваний.

Этот протокол разработан специально для массовых наборов данных по РНК-секвенциям, полученных из человеческих тканей с высокой целостностью РНК (RIN ≥ 7.0) и достаточной колькостью входной РНК (≥ 500 нг на образец). Для обеспечения надёжного выполнения шагов выравнивания и количественной оценки рекомендуется локальная рабочая станция с не менее 10-ядерным процессором, 32 ГБ оперативной памяти и минимум 200 ГБ свободного дискового пространства. Основываясь на этих требованиях, протокол обеспечивает эффективный и удобный аналитический рабочий процесс, включая подробные операционные инструкции и стандартизированные конфигурации параметров, чтобы удовлетворить потребности исследователей, анализирующих крупномасштабные транскриптомические данные.

Protocol

Для демонстрационных целей публично доступный набор данных, PRJNA1023502 созданный Лань Бай и др., был использован для иллюстрации каждого шага как в верхнем, так и в нижнем^{анализе 20}. Поскольку этот набор данных основан на базе данных NCBI SRA с открытым доступом, дополнительные разрешения или этические одобрения не требуются. См. таблицу материалов для проверки всех необходимых программ и версий R-пакета. Общедоступный набор данных PRJNA1023502 включает 6 образцов РНК-секвенации печени без NASH, 6 NAFL и 6 NASH. В этом протоколе набор данных использовался для демонстрации всех этапов рабочего процесса массовой РНК-секвенации, включая извлечение данных из базы данных SRA, контроль качества (fastp), выравнивание (HISAT2), количественную оценку (featureCounts), а также анализ дифференциальной экспрессии и функционального обогащения.

1. Установка набора инструментов SRA

1. Посетите официальный сайт SRA Toolkit и скачайте версию 3.2.1.

2. Публичная загрузка данных

1. Получите номер SRA.
   1. Получите доступ к нему в дополнении^{статьи 20}, в разделе доступности данных или ища ключевые слова в базе данных NCBI SRA.
   2. Введите предвыборку <номер SRA> в терминале, чтобы скачать его.

3. Генерация матрици по количеству генов

1. Определите пути и настройки, как описано ниже:
   REFERENCE=~/reference/human/GRCh38/GRCh38.primary_assembly.genome.fa
GTF=~/reference/human/GRCh38/gencode.v44.annotation.gtf
INDEX=~/reference/human/GRCh38/GRCh38_index
FASTQ_DIR=~/SRA_tutorial/fastq
OUT_FASTP=~/RNAseq/fastp
OUT_HISAT2=~/RNAseq/hisat2
OUT_COUNTS=~/RNAseq/counts
​mkdir -p $FASTQ_DIR $OUT_FASTP $OUT_HISAT2 $OUT_COUNTS
   1. Скачайте человеческий эталонный геном (первичная сборка GRCh38) и соответствующий файл аннотации гена (gencode.v44, комплексный набор аннотаций генов) из официальной базы данных GENCODE (https://www.gencodegenes.org/human/).
      for f in SRR*; do [[ ! $f =~ \.sra$ ]] && mv "$f" "$f.sra"; done
​for f in SRR*; do [[ ! $f =~ \.sra$ ]] && mv "$f" "$f.sra"; done
2. Переименуйте файлы SRA.
   1. Убедитесь, что все SRA-файлы используют расширение .sra для точного распознавания и обработки с помощью последующих инструментов.
   2. Выполните следующие команды:
      for f in SRR*; do [[ ! $f =~ \.sra$ ]] && mv "$f" "$f.sra"; done
3. Выполните следующую команду для преобразования SRA-файлов в формат FASTQ:
   for f in *.sra; do fasterq-dump "$f" --split-files -O $FASTQ_DIR - e 20; done
4. Постройте индекс HISAT2 на основе эталонного генома.
   hisat2-build $REFERENCE $INDEX
5. Обрабатывайте файлы FASTQ с визуальным мониторингом прогресса.
   for fq in $FASTQ_DIR/*.fastq; do
​sample=$(basename "$fq" .fastq)
   1. Для данных секвенирования парных концов выполните следующую команду
      for fq1 in $FASTQ_DIR/*_1.fastq; do
sample=$(basename "$fq1" _1.fastq)
​fq2=$FASTQ_DIR/${sample}_2.fastq
   2. Выполните контроль качества и фильтрацию с помощью fastp и выводите результаты в каталог $OUT_FASTP. Для одноконцевых данных секвенирования выполните следующую команду:
      fastp \
-i "${fq}" \
-o $OUT_FASTP/${sample}.clean.fastq \
-h $OUT_FASTP/${sample}.html \
-j $OUT_FASTP/${sample}.json \
-w 20
      Для данных секвенирования парных концов выполните следующую команду:
      fastp \
-i "${fq}" \ -I "$fq2" \
-o $OUT_FASTP/${sample}_1.clean.fastq \
-O $OUT_FASTP/${sample}_2.clean.fastq \
-h $OUT_FASTP/${sample}.html \
-j $OUT_FASTP/${sample}.json \
​-w 20
   3. Для каждого образца проверьте HTML-отчёт для визуального контроля качества и обратитесь к отчету JSON для получения структурированных метрик, включая качество чтения, распределение GC, избыточность чтения и загрязнение адаптеров. Убедитесь, что все образцы соответствуют следующим порогам качества: Q30 ≥ 85%, уровень загрязнения адаптера < 5%, и отсутствие аномальных паттернов GC. Храните все выходные файлы в директории $OUT_FASTP.
   4. Выполните выравнивание последовательностей с помощью HISAT2, генерируйте SAM-файлы и выводите их в каталог $OUT_HISAT2. Для одноконцевых данных секвенирования выполните следующую команду:
      hisat2 -p 20 \ -x $INDEX \-U $OUT_FASTP/${sample}.clean.fastq \
-S $OUT_HISAT2/${sample}.sam
      Для данных секвенирования парных концов выполните следующую команду:
      hisat2 -p 20 \-x $INDEX \-1 $OUT_FASTP/${sample}_1.clean.fastq \
-2 $OUT_FASTP/${sample}_2.clean.fastq \
​-S $OUT_HISAT2/${sample}.sam
   5. Выполните выравнивание HISAT2 считываний с фильтром качества с человеческим эталонным геномом, используя индекс, построенный из основного FASTA-файла GRCh38. Предоставить примеры команд для данных секвенирования с одним концом и парным концом.
   6. Конвертируйте SAM в BAM, сортируйте и индексируйте.
      samtools view -@ 20 -bS $OUT_HISAT2/${sample}.sam \
| samtools sort -@ 20 -o $OUT_HISAT2/${sample}.sorted.bam
samtools index $OUT_HISAT2/${sample}.sorted.bam
​done
   7. Сгенерируйте SAM-файл для каждого образца и конвертируйте его в отсортированный и индексированный BAM-файл с помощью SAMtools. Для высококачественных человеческих транскриптомных образцов (число целостности РНК, RIN ≥ 8,0) обеспечите, чтобы общая скорость выравнивания превышала 85%; для типичных объёмных образцов РНК-секвенции (RIN ≥ 7.0) считайте приемлемыми скорости выравнивания ≥ 70%.
   8. Выполните количественную оценку генов с помощью featureCount.
      featureCounts -T 20 -p -s 0 \
-a $GTF \
-o $OUT_COUNTS /${sample}.counts.txt \
$OUT_HISAT2/${sample}.sorted.bam
Done
   9. Проверьте выходной файл с разделением по табуляциям (*.counts.txt) и его сводный отчет (*.counts.txt.summary), сгенерированный featureCounts для каждого выборки. Убедиться, что частота распределения чтения соответствует типичному порогу ≥70% для человеческой массовой РНК-секвенции; Значительно более низкая частота может указывать на несовпадающую ориентацию нитей, проблемы с аннотированием или плохое качество выравнивания. Используйте параметр -s 0 для этого не-цепочно-специфического набора данных РНК-секвенации. Для библиотек, специфичных для ниточки, замените -s 0 на -s 1 или -s 2 в команде.
6. Сгенерируйте матрицу подсчетов генов, как описано ниже.
   1. Инициализируйте матрицу подсчетов с помощью идентификаторов генов и подсчётов из первого образца.
      cut -f1 $(ls $OUT_COUNTS/*.counts.txt | head -1) > all_counts.txt
   2. Итеративно добавляйте счёты из каждой выборки в одну матрицу.
      for f in $OUT_COUNTS/*.counts.txt; do
cut -f7 "$f" | paste all_counts.txt - > tmp && mv tmp
all_counts.txt
​done
   3. Добавьте строку заголовка с идентификаторами образцов в матрицу подсчёта.
      samples=$(ls *.counts.txt | sed 's/.counts.txt//' | paste -sd "\t")
echo -e "Geneid\t$samples" | cat - all_counts.txt > counts_matrix.txt
   4. Извлекать длины генов из файла GTF (сумма длин экзонов на ген).
      awk '$3=="exon"{match($0,/gene_id "([^"]+)"/,a); if(a[1]!=""){len=$5-$4+1; gene_len[a[1]]+=len}} END{print "GENE_ID\tLENGTH"; for(g in gene_len) print g"\t"gene_len[g]}' \$GTF > gene_length.txt

4. Обработка матриц исходного счёта и аннотация генов

1. Инициализация среды R и загрузка соответствующего пакета R.
   1. Используйте функцию install.packages() для установки tidyverse, ggplot2, ggrepel, RColorBrewer, ggridges, FactoMineR. Пакет tidyverse используется для обработки данных и построения графиков; Пакет ggplot2 используется для визуализации; Пакет ggrepel используется для неперекрывающихся текстовых меток; пакет RColorBrewer используется для цветовых палитр; пакет Grridges применяется для участков гребней; пакет FactoMineR используется для PCA и многомерного анализа.
   2. Используйте функцию BiocManager::install() для установки biomaRt, DESeq2, clusterProfiler. Пакет biomaRt используется для аннотации генов от Ensembl; пакет DESeq2 используется для анализа дифференциальных выражений; пакет clusterProfiler используется для функционального анализа обогащения.
   3. Прочитайте исходную матрицу подсчёта, сгенерированную на шаге 2, используйте biomaRt для сопоставления Ensembl ID с именем гена HGNC для последующего анализа, удалите дублирующиеся столбцы Geneid и сгенерируйте очищенную матрицу счёта. Используйте исходную матрицу подсчёта (counts_matrix.csv) в качестве входных данных, с идентификаторами генов Ensembl в виде строк и образцами в виде столбцов. Сохраните очищенную матрицу подсчёта как выход (clean_counts_SRA.csv), с символами генов HGNC в виде строк и образцами в виде столбцов.
      mart <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
id_map <- getBM(attributes = c("ensembl_gene_id", "hgnc_symbol"),
filters = "ensembl_gene_id",
values = exprSet$GeneID,
mart = mart)
exprSet <- exprSet %>%
left_join(id_map, by = c("GeneID" = "ensembl_gene_id")) %>%
filter(!is.na(hgnc_symbol), hgnc_symbol != "") %>%
distinct(hgnc_symbol, .keep_all = TRUE) %>%
column_to_rownames("hgnc_symbol")

5. Количественная оценка экспрессии генов

ПРИМЕЧАНИЕ: См. дополнительный файл 1 для подробного сценария.

1. Выполните следующую команду для расчёта фрагментов на килобазу транскрипта на миллион отображённых чтений (FPKM) и транскриптов на миллион (TPM).
   counts <- read.csv("output/clean_counts_SRA.csv", header=TRUE, row.names=1)
gene_len <- read.delim("data/gene_length.txt", header=FALSE, col.names=c("gene_symbol","length"))
gene_len <- gene_len %>% distinct(gene_symbol, .keep_all=TRUE)
rownames(gene_len) <- gene_len$gene_symbol
gene_len <- gene_len[match(rownames(counts), gene_len$gene_symbol),]
length_bp <- gene_len$length
fpkm <- (counts / length_bp) * 1e9 / colSums(counts)
write.csv(fpkm, "output/clean_fpkm_SRA.csv")
tpm <- (counts / length_bp) / colSums(counts / length_bp) * 1e6
write.csv(tpm, "output/clean_tpm_SRA.csv")

6. Кластеризация выборок и визуализация различий

1. Оценивайте качество данных, изучая кластеризацию выборок с помощью PCA. См. дополнительный файл 1 для подробного сценария.
2. Для визуализации PCA используйте код ниже.
   gene.pca <- PCA(exprSet, ncp = 2, scale.unit = TRUE, graph = FALSE)
ggplot(pca_sample, aes(x = Dim.1, y = Dim.2)) +
geom_point(aes(color = group)) +
labs(x = paste('PC1:', pca_eig1, '%'),
y = paste('PC2:', pca_eig2, '%'))

7. Анализ дифференциальных выражений и визуализация результатов

ПРИМЕЧАНИЕ: См. дополнительный файл 1 для подробного сценария.

1. Выполните следующие команды для создания набора данных DESeq2 и установите порог. До моделирования сохранялись только гены с общим количеством > 1 во всех образцах. Идентифицировать дифференцированно экспрессированные гены (DEGs) с использованием порогов значимости скорректированных p-значений < 0,05 и |log2FC| > 0,5.
   dds <- DESeq(DESeqDataSetFromMatrix(countData = exprSet, colData = colData, design = ~group)); sizeFactors(dds); res <- results(dds); dds <- dds[rowSums(counts(dds)) > 1,]
dd1 <- results(dds, contrast = contrast, alpha = 0.05)
dd2 <- lfcShrink(dds, contrast = contrast, res = dd1, type = "ashr")
2. Визуализируйте дифференциально экспрессируемые гены с помощью вулканического графика.
   ggplot(data = data, aes(x = log2FoldChange, y = -log10(padj))) +
geom_point(aes(color = group), alpha = 1, size = 1.2) +
geom_hline(yintercept = -log10(0.05), lty = 4) +
geom_vline(xintercept = c(-0.5, 0.5), lty = 4) +
geom_text_repel(data = subset(data, abs(log2FoldChange) >= 1.5 & padj < 0.05),
aes(label = gene_id))

8. Проведение анализа функционального обогащения и визуализации

ПРИМЕЧАНИЕ: См. дополнительный файл 1 для подробного сценария.

1. Визуализируйте результаты обогащения KEGG.
   EGG <- enrichKEGG(gene = gene$ENTREZID, organism = 'hsa',
pvalueCutoff = 0.05, qvalueCutoff = 0.05)
ggplot(symboldata, aes(richFactor, Description)) +
geom_point(aes(color = p.adjust, size = Count))
2. Визуализируйте результаты обогащения GO. Настройте визуализацию, изменяя разрешение дуги (n=500), цветовые отображения и расположение фасетей в соответствии с экспериментальными спецификациями.
   ego <- enrichGO(gene = gene$ENTREZID, OrgDb = "org.Hs.eg.db", ont = "ALL",
pvalueCutoff = 0.05, qvalueCutoff = 0.05, pAdjustMethod = "BH")
ggplot(df) +
ggforce::geom_link(aes(x = 0, y = Description, xend = -log10(p.adjust),
yend = Description, color = ONTOLOGY), n = 500, show.legend = FALSE) +
facet_wrap(~ONTOLOGY, scales = "free", ncol = 1)
3. Визуализируйте результаты анализа обогащения генных наборов (GSEA). Отрегулировать количество путей, масштабирование плотных гребней и прозрачность слоёв, чтобы учесть разные размеры эффектов.
   genelist <- sort(res$log2FoldChange, decreasing = TRUE)
names(genelist) <- rownames(res)
hallmarks <- read.gmt('resource/h.all.v2023.2.Hs.symbols.gmt')
y <- GSEA(genelist, TERM2GENE = hallmarks, pvalueCutoff = 0.05)
gsearesult <- yd %>% arrange(desc(NES)) %>% slice_head(n = 10)
ggplot(gsearesult, aes(x = logFC, y = Description, fill = -log10(pvalue))) +
geom_density_ridges(alpha = 0.8, scale = 0.8) +
geom_point(aes(size = abs(NES), x = -0.4, color = NES)) +
scale_fill_distiller(palette = 'Spectral') +
scale_color_distiller(palette = 'Reds') +
scale_size_continuous(range = c(2, 6))

Results

Рабочий процесс восходящего анализа для массового РНК-секвенации иллюстрирован на рисунке 1A. Этот рабочий процесс последовательно выполняет следующие ключевые шаги на платформе Linux: во-первых, строгий контроль качества исходных данных секвенирования осуществляется с помощью fastp для удаления низкокачественных чтений и последовательностей адаптеров; впоследствии HISAT2 выравнивает высококачественные чтения с эталонным геномом, при этом Samtools конвертирует и сортирует файлы выравнивания;...

Discussion

Анализ объёмных данных РНК-секвенации характеризуется как междисциплинарная задача, объединяющая геномику, биоинформатику, статистику и информатику. Полный аналитический рабочий процесс включает несколько этапов на начальном и следующем этапе, включая предварительную обработку исходных данных, контроль качества, выравнивание последовательностей, количественную оценку на уровне генов, нормализацию данных, анализ дифференциальной экспрессии и биологическую интерпретацию. Среди этих этапов ...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
biomaRt	Биопроводник	2.64.0	Аннотация гена из Ensembl
clusterProfiler	Биопроводник	4.16.0	Анализ функционального обогащения
DESeq2	Биопроводник	1.48.1	Анализ дифференциальных выражений
FactoMineR	AgroParisTech	2.11.0	PCA и многомерный анализ
FASTP	OpenGene	1.0.1	Контроль качества и фильтрация данных FASTQ
FeatureCounts	Отдел биоинформатики, Институт медицинских исследований Уолтера и Элизы Холл	2.0.0	Посчитайте количество чтений, отображённых на каждый ген для количественной оценки экспрессии генов
ggplot2	Положить	3.5.2	Визуализация данных
ggrepel	Камиль Словиковский	0.9.6	Неперекрывающиеся текстовые метки
Гриджи	Клаус О. Уилке	0.5.6	Создание графиков гребней
HISAT2	Университет Джонса Хопкинса	2.2.1	Выровнять отфильтрованные высококачественные чтения с эталонным геномом
R	R Core Team	4.5.0	Среда для вычислений, анализа и визуализации данных
RColorBrewer	Эрих Нойвирт	1.1.3	Цветовые палитры для построения графиков
samtools	Крупномасштабный рабочий поток по геномике	1.22.0	Конвертировать и обрабатывать SAM-файлы для эффективного поиска и доступа
Набор инструментов SRA	Национальный центр биотехнологической информации	3.2.1	Получение и предварительная обработка исходных данных секвенирования из базы данных NCBI SRA