Rup (RNA-seq Usability Assessment Pipeline) - Controle de qualidade para experimentos de RNA-seq em massa em eucariotos

Resumo

Este protocolo permite o controle de qualidade inicial para experimentos de RNA-seq para biólogos de laboratório úmido com experiência limitada em bioinformática.

Resumo

As abordagens modernas na ciência molecular de plantas geralmente exigem experimentos de RNA-seq em massa, por exemplo, para rastrear mudanças globais nos transcriptomas após os tratamentos ou para identificar os principais componentes das vias regulatórias. Consequentemente, diversas áreas das ciências vegetais dependem de dados de RNA-seq em massa reprodutíveis e de alta qualidade para o avanço científico. No entanto, a partir de nossa experiência, muitas vezes falta conhecimento e aplicação de medidas de controle de qualidade em conjuntos de dados de RNA-seq. Aqui, apresentamos o Rup (pipeline de avaliação de usabilidade de RNA-seq) para o controle de qualidade de dados de RNA-seq em massa, para análises de expressão gênica subsequentes, que é autônomo e prontamente aplicável para biólogos de laboratório úmido com conhecimento básico de R. Rup ajuda a discriminar entre dados de sequenciamento de alta qualidade, adequados para experimentos de expressão gênica a jusante, e aqueles inadequados para análises gerais posteriores. O Rup inclui testes para vários problemas comumente encontrados, como números de leitura ou mapeamento insuficientes, identificação de contaminações, quantificação de frações de rRNA nos dados totais de RNA-seq, teste de similaridade replicada, usando dados reais para demonstração e oferecendo visualização intuitiva. O Rup fornece um conjunto de ferramentas para identificar deficiências experimentais antes da análise padronizada do transcriptoma, melhorando assim a qualidade dos dados para pesquisadores individuais e para o campo. Isso aumenta a confiança na análise de dados de RNA-seq em massa e fornece uma base para futuras diretrizes que definem critérios mínimos de controle de qualidade, melhorando assim a confiabilidade e a transparência dos dados de RNA-seq publicados. Os dados de rup e teste estão disponíveis em https://github.com/oliverrupp/rup.

Introdução

Experimentos de transcriptômica (RNA-seq) interrogam de forma abrangente as assinaturas transcricionais que moldam os fenótipos. Essa abordagem tornou-se insubstituível na genética molecular de plantas para identificar genes únicos ou co-expressos e processos biológicos envolvidos, por exemplo, no desenvolvimento, interação com patógenos de plantas ou resistência ao estresse abiótico^{1, 2,3}. Avanços recentes na tecnologia de RNA-seq aumentaram a especificidade e permitiram a detecção de diferentes isoformas e variantes em resolução de base única, permitindo a identificação de variações de sequência de indels maiores a polimorfismos de nucleotídeo único (SNPs). Os dados obtidos por RNA-seq são caracterizados por uma ampla faixa dinâmica, facilitando a detecção de transcritos altamente abundantes e pouco expressos, e requerem condições experimentais apropriadas para consistência. Além disso, os conjuntos de dados de RNA-seq são grandes e, portanto, computacionalmente intensivos para analisar e armazenar, envolvendo gerenciamento eficiente de dados e extensos recursos de computação. O RNA-seq requer entrada de alta qualidade em todas as etapas do fluxo de trabalho, pois os pontos fracos em qualquer estágio podem se propagar e comprometer os resultados. A má integridade do RNA ou problemas técnicos durante a preparação da biblioteca levarão a vieses e redução da precisão. Os parâmetros para obter leituras brutas de alta qualidade podem variar entre as instalações de sequenciamento de próxima geração (NGS).

De acordo com nossa experiência, recomendamos o uso apenas de RNA com um número de integridade de RNA (RIN) acima de 7, o que é indicativo de estrutura de mRNA amplamente intacta, como entrada para a preparação da biblioteca de RNA-seq. Para um sequenciamento bem-sucedido, aproximadamente 2 μg de RNA total a uma concentração de 50–200 ng/μL são necessários para protocolos de preparação de biblioteca padrão. A pureza do RNA deve ser confirmada por uma razão OD_260/280 entre 1,8 e 2,1 e uma razão OD_260/230 maior que 1,5 usando um espectrofotômetro. Profundidade de sequenciamento insuficiente, baixas taxas de mapeamento ou desalinhamento com o genoma de referência podem distorcer ainda mais a expressão gênica e a quantificação do splicing. Além disso, o desenho experimental inadequado, como a falta de discriminação entre amostras de diferentes tecidos ou tratamentos e a alta variabilidade entre as repetições, pode introduzir ruído e diminuir a reprodutibilidade. Embora muitos laboratórios analisem transcriptomas rotineiramente, controles de qualidade rigorosos das primeiras etapas essenciais da análise geralmente não são relatados ou podem estar totalmente ausentes das publicações. Isso pode levar à interpretação excessiva dos resultados derivados da análise do transcriptoma e, consequentemente, a resultados irreprodutíveis.

Aqui, fornecemos um fluxo de trabalho para controle de qualidade das etapas iniciais necessárias para análises de transcriptoma de alta qualidade de perfis de mRNA para medir alterações transcricionais. Nosso objetivo é permitir que biólogos de laboratório úmido com conhecimento limitado em bioinformática avaliem seus dados transcriptômicos primários. O Rup (Figura 1) é acessível a pesquisadores familiarizados com o conhecimento básico de R. A execução do fluxo de trabalho apresentado aqui fornecerá aos pesquisadores uma compreensão detalhada de seus dados primários, incluindo suas possíveis limitações para análise subsequente. Até onde sabemos, falta um pipeline prático de avaliação de dados de transcriptoma primário em combinação com diretrizes para distinguir dados de alta qualidade de dados de baixa qualidade.

Figura 1: Fluxo de trabalho do pipeline de controle de qualidade RNA-seq in silico . Os arquivos de entrada para controle de qualidade são derivados de dados de RNA-seq gerados pelo sequenciamento de material vegetal, bem como conjuntos de dados disponíveis publicamente. O pipeline R fornecido avalia a qualidade da sequência por meio de três abordagens principais: qualidade de sequenciamento, qualidade de mapeamento e qualidade de replicação. Várias métricas de qualidade são calculadas e os resultados estatísticos são visualizados usando pacotes R comuns, como ggplot2 e pheatmap (por exemplo, gráficos de barras e mapas de calor). Clique aqui para ver uma versão maior desta figura.

Os pipelines de avaliação relatados anteriormente exigiam dados pré-processados (por exemplo, alinhamento de leitura como arquivos .bam), não cobrem todas as métricas ou não são mais mantidos ^4,5,6. A vantagem do fluxo de trabalho apresentado aqui reside em sua abrangência, consolidando os problemas de controle de qualidade mais comuns em um único pipeline. Além disso, fornecemos exemplos de dados de alta qualidade adequados para todas as aplicações de análise downstream, mas também exemplos de dados de baixa qualidade e discutimos suas limitações específicas para análise posterior. Vários relatórios publicados anteriormente descrevem o propósito das ferramentas de análise de RNA-seq e fornecem avaliações comparativas de seu desempenho ^7,8. No entanto, o controle de qualidade de RNA-seq e os padrões de relatórios de metodologia não são padronizados e agravam a reprodutibilidade e as interpretações biologicamente significativas dos experimentos transcriptômicos.

O Rup integra ferramentas padrão de alta qualidade, análise de controle de qualidade e visualização dos resultados. O Rup requer leituras de sequenciamento bruto e um genoma anotado como entrada e é executado no Mac OS, Linux e "Windows Subsystem for Linux" (WSL) em sistemas Windows. Ele integra ferramentas para a qualidade do sequenciamento e quantifica o mapeamento de leitura para um genoma, inclui a medição do conteúdo de rRNA nas amostras e oferece correlação de amostra para avaliar a correlação replicada. Os dados do teste foram obtidos por meio de microdissecção a laser do tecido do meristema central de dois estágios diferentes da espécie vegetal Eschscholzia californica, um protocolo de entrada ultrabaixa para preparo de bibliotecas, e sequenciados em Novaseq 6000.

Protocolo

O Rup pode ser executado como um único script com o mínimo de arquivos de entrada necessários. O pipeline, como mostrado aqui, requer dados de sequenciamento de RNA de extremidade emparelhada da Illumina. Um pipeline ajustado para sequenciamento de extremidade única com as mesmas etapas de análise também é depositado no repositório GitHub. Apenas a sequência do genoma como um arquivo fasta, o modelo de gene e as anotações de rRNA como arquivos gtf e o sequenciamento bruto são lidos como arquivos fastq.gz são necessários. Certifique-se de que os arquivos de genoma e anotação sejam fornecidos como genome.fa, annotation.gtf e rRNA.gtf em uma pasta especificada na variável reference_folder. Quando todos os arquivos de sequenciamento como arquivos .fq.gz são salvos na variável read_file_folder, o Rup pode ser executado da seguinte maneira:

1. Preparações:

1. Instale todos os pacotes R necessários usando o Bioconductor.

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

2. Configure todos os arquivos necessários.
   1. Crie uma pasta de origem que conterá todos os arquivos necessários. Adicione a sequência do genoma de referência no arquivo fasta chamado "reference/genome.fa" à pasta de origem. Forneça a anotação do modelo de gene como um arquivo gtf: "reference/annotation.gtf", Opcionalmente, forneça a anotação do gene rRNA como arquivo gtf: "reference/rRNA.gtf".
   2. Adicione todas as leituras de sequenciamento como arquivos fastq à pasta "reads". Certifique-se de que todos os nomes de arquivo fastq sigam o mesmo padrão; _1.fastq.gz e _2.fastq.gz para as leituras direta e reversa.

      # 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. Defina o parâmetro de acordo com os recursos de computação.

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

4. Defina o parâmetro de acordo com o método de sequenciamento.
    

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

   
   
   NOTA: O comprimento mínimo de leitura não deve ser menor que 25 bp, pois leituras menores podem causar a quebra do mapeamento. Valores maiores resultarão em leituras mapeadas de forma mais exclusiva, mas também em mais leituras descartadas durante o corte, dependendo da qualidade do sequenciamento. Mais detalhes sobre esses parâmetros são descritos nos suplementos.

2. Avaliação da qualidade do sequenciamento

NOTA: Esta etapa criará um gráfico de barras mostrando o número de leituras antes e depois de cortar cada amostra.

1. Execute o fastqc⁹, que é uma ferramenta para controle de qualidade geral de dados de sequência de alto rendimento, fornecendo informações sobre o número total de leituras sequenciadas, o comprimento de leitura, o conteúdo médio de GC, contaminação do adaptador de sequenciamento e sequências super-representadas (por exemplo, leituras de rRNA). Analise a qualidade do sequenciamento de acordo com os relatórios de qualidade da base e da sequência. Execute fastqc nos arquivos de sequenciamento brutos usando o pacote R fastqcr¹⁰ e coloque todos os arquivos de entrada no mesmo diretório (read_file_folder). Resuma os arquivos de saída do fastqc_folder em um objeto R usando a função qc_aggregate.
   NOTA: Os resultados serão salvos no diretório fastqc_folder. Os arquivos .html criados pelo fastqc na pasta de saída podem ser acessados para obter estatísticas adicionais.
    

   # 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. Com todos os prefixos de amostra armazenados no sample_prefixes vetorial, itere sobre os prefixos e execute a ferramenta de corte de qualidade fastp em todas as amostras. Trim para remover sequências de adaptadores, primers, bases de baixa qualidade, sequências poli-A/T e é executado pelo Fastp¹¹, que está disponível no pacote R Rfastp¹². Os arquivos fastq cortados serão armazenados na pasta 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. Execute fastqc novamente nas leituras cortadas.
    

   # 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. Avalie o corte coletando o número de leituras antes e depois do corte para todas as amostras. Crie um gráfico de barras com os números lidos usando 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)

Figura 2: Resultados de sequenciamento e corte. Contagens de leitura antes (vermelho) e depois do corte (verde). O s2_r1 de amostra já tem uma contagem de leitura baixa antes do corte, enquanto o corte removeu uma grande fração das leituras de s2_r2 de amostra. Todas as outras amostras mostram perda aceitável de leituras de corte. Clique aqui para ver uma versão maior desta figura.

3. Qualidade do mapeamento

NOTA: Esta etapa calcula gráficos de barras para discriminar leituras mapeadas únicas (por exemplo, transcritos de genes codificadores de proteínas) de leituras multimapeadas (por exemplo, leituras de rRNA) e leituras não mapeadas (por exemplo, contaminações).

1. Use o pacote Rsubread¹⁴ para mapear as leituras para o genoma de referência. Primeiro, crie um índice do arquivo fasta do genoma de referência (genome_fasta_file).
   NOTA: Esta etapa é realizada uma vez para cada genoma de referência.
    

   # 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. Itere em todas as amostras e alinhe as leituras ao genoma de referência usando a função align() do pacote Rsubread. Essa função cria arquivos .bam na pasta de saída.
    

   # 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. Conte as leituras de cada gene usando a função featureCounts() do pacote Rsubread. Certifique-se de que o arquivo de anotação (annotation_file) esteja no formato GTF. A contagem só lê com uma única correspondência com o genoma.
    

   # 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. Conte o mapeamento de leituras para genes de rRNA para avaliar o conteúdo de rRNA nas amostras. Permita que leituras multimapeadas sejam contadas aqui. Especifique os loci do gene rRNA em um arquivo 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. Colete as estatísticas de atribuição de leitura criadas por featureCounts. As estatísticas incluem o número de alinhamentos para cada amostra nas diferentes categorias (por exemplo, Atribuído, Não Mapeado, MultiMapeado).
    

   # 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. Colete as estatísticas da atribuição de 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. Plote as estatísticas de mapeamento de leitura como um gráfico de barras.
    

   # 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. Classifique os genes em grupos de acordo com o número de leituras atribuídas a eles e plote os resultados como um gráfico de barras.
    

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

Figura 3: Leia a visão geral e as estatísticas do mapeamento. A maioria das amostras mostra um alto número de leituras atribuídas (lado esquerdo, marrom) e um baixo número de leituras de rRNA (lado direito, azul). A amostra s2_r3 tem uma quantidade excepcionalmente alta de leituras multimapeadas (lado esquerdo, verde) correspondendo a um alto número de leitura de rRNA (lado direito). A amostra s2_r4 mostra um alto número de leituras não mapeadas em combinação com um número esperado de leitura de rRNA, sugestivo de contaminação com leituras de outro organismo. Clique aqui para ver uma versão maior desta figura.

4. Replique a qualidade.

1. Certifique-se de que os replicados da mesma amostra ou tecido apresentam uma correlação mais alta entre si do que com replicados de amostras ou tecidos diferentes. Plote um mapa de calor agrupado das correlações replicadas para visualizar os dados usando as contagens de leitura brutas ou as contagens normalizadas de transcrições por milhão (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")

Figura 4: Classificação da contagem de leitura de genes. Todos os genes em uma amostra são classificados em uma das cinco categorias com base no número de leituras atribuídas. Em vermelho está o número de genes sem nenhuma leitura atribuída. Amostras com um baixo número total de leituras atribuídas (s2_r1 a s1_r4) têm um número maior de genes com 10 a 100 leituras atribuídas e um número menor de genes com mais de 1000 leituras. Genes com baixa expressão podem ser perdidos nessas amostras. Clique aqui para ver uma versão maior desta figura.

Resultados

O pipeline é completamente implementado como um script R e foi testado nos sistemas operacionais Linux e Mac OS. Os usuários do Windows podem usar o WSL (Subsistema do Windows para Linux). O código e os dados de teste estão disponíveis como um repositório GitHub: https://github.com/oliverrupp/rup. Os dados de sequenciamento estão disponíveis no projeto ENA EBI PRJEB96400.

Dez amostras foram criadas artificialmente a partir de duas amostras reais para exemplificar diversos problemas que podem ser encontrados durante o controle de qualidade do RNA-seq em massa usando Rup para análise. O s2_r1 de amostra foi projetado para exibir um baixo número total de leituras, s2_r2 de amostra contém uma grande fração de leituras de baixa qualidade a serem descartadas pelo processo de corte. A amostra s2_r3 inclui uma grande fração de leituras de rRNA e a s2_r4 de amostra inclui leituras de contaminantes que não foram mapeadas para o genoma de referência. Os nomes das amostras s1_r5 e s2_r5 foram trocados para ilustrar a baixa correlação de replicações.

A seção 2 do protocolo permite a identificação de amostras com baixo número de leitura antes ou depois do corte (Figura 2). O gráfico de barras mostra o número de leitura mais baixo no s1_r1 de amostra antes e depois do corte. Aqui, o número inicial de leitura foi baixo. O baixo número de leitura de s1_r2 de amostra após o corte sugere uma grande quantidade de sequências adaptadoras, erros de sequenciamento, sequências de primers, trechos poli-A/T sequenciados que foram removidos durante o processo de corte. A degradação do RNA de entrada também pode reduzir o número de leituras de alta qualidade.

A seção 3 do protocolo identifica problemas nas atribuições de leitura. A Figura 3 mostra o número elevado de leituras multimapeadas no s2_r3 de amostra (verde), bem como o alto número de leituras de rRNA. A contaminação do s2_r4 da amostra é aparente pela grande fração de leituras não mapeadas para o genoma de referência (rosa). Essas leituras não se correlacionam com o rRNA, mas com sequências de um organismo não-alvo. A Figura 4 mostra problemas gerais associados ao baixo número de leituras atribuídas em transcriptomas. Em amostras com baixa taxa de leituras mapeadas exclusivamente para o genoma de referência (s2_r1 a r4), apenas cerca de um terço dos genes têm mais de 100 leituras atribuídas, enquanto nas outras amostras, cerca de metade dos genes se enquadram nessas classes. Leituras de genes com expressão muito baixa podem não ser encontradas nas amostras s2_r1 a r4 e, consequentemente, a análise de expressão comparativa com essas amostras será altamente não confiável e deve ser evitada.

A seção 4 do protocolo pode ser usada para identificar outliers replicados. Espera-se que as réplicas da mesma amostra/condição/tecido apresentem uma correlação mais alta entre si do que as réplicas de outras amostras/condições/tecidos. A Figura 5 mostra um mapa de calor de correlação das duas amostras (S1 e S2) com cinco repetições cada. Os dendrogramas no topo e ao lado do gráfico mostram dois clusters com cinco repetições em cada cluster. O cluster esquerdo contém quatro réplicas da amostra 1 e uma réplica da amostra 2 (s2_r5), o cluster direito contém quatro réplicas da amostra 2 e uma réplica da amostra 1 (s1_r5). Nesse caso, quando os nomes de amostra s1_r5 e s2_r5 são trocados novamente, cada cluster contém todas as réplicas de uma amostra, o que pode indicar um erro de rotulagem de replicação. Outras razões para os replicados não agruparem podem ser a falta de diferenciação entre as amostras/condições/tecidos, ou o agrupamento de replicados apenas por razões de qualidade de sequenciamento. Este último pode ocorrer quando todas as réplicas com contagens de leitura excepcionalmente baixas ou aquelas com contagens de leitura excepcionalmente altas após o corte e o mapeamento formam um cluster.

Figura 5: Mapa de calor de correlação de amostra. O mapa de calor de correlação de amostra é baseado em valores de TPM transformados em log₂ e mostra dois clusters distintos de cinco amostras cada. Espera-se que as réplicas biológicas/técnicas mostrem uma correlação mais alta entre si do que as de outros tecidos/tratamentos. O cluster esquerdo contém quatro réplicas da amostra 1 e uma réplica da amostra 2 (s2_r5), o cluster direito contém quatro réplicas da amostra 2 e uma réplica da amostra 1 (s1_r5). As réplicas individuais devem ser verificadas quanto a possíveis efeitos de troca de amostras ou lotes para explicar seu agrupamento no mapa de calor. Clique aqui para ver uma versão maior desta figura.

Tabela 1: Comparação de pipelines de avaliação de qualidade de RNA-seq. Para uma análise detalhada, consulte o Arquivo Suplementar 1. Clique aqui para baixar esta tabela.

Arquivo Suplementar 1: Seleção de parâmetro e referência. As análises mostram a influência da seleção de parâmetros e referências nos resultados gerais do CQ. Clique aqui para baixar este arquivo.

Discussão

A qualidade das análises de expressão gênica diferencial depende muito de dois fatores: o número de leituras sequenciadas em cada amostra e réplica¹⁶ e o número de réplicas por amostra^17,18. Aqui, apresentamos o pipeline Rup amigável para determinar o número de leituras adequadas para quantificação da expressão gênica em cada replicação. Diferentes métricas permitem que os pesquisadores entendam por que as réplicas mostram baixos números de leitura atribuídos e identifiquem problemas na correlação da amostra. Embora o Rup tenha sido desenvolvido para avaliação da qualidade de amostras de RNA-seq de plantas, ele é igualmente adequado para outros organismos eucarióticos, não exigindo ajustes adicionais para amostras não vegetais (Arquivo Suplementar 1).

A primeira etapa do Rup determina a contagem total de leitura antes e depois do corte de leitura. O número de leituras sequenciadas determina o número de genes detectáveis e, se o número de leituras for muito baixo, muitos genes diferencialmente expressos permanecerão não detectados. Uma grande fração de leituras descartadas durante o processo de corte e filtragem de qualidade pode indicar degradação do RNA, o que pode levar a uma subestimação da expressão de genes altamente degradados. No entanto, se uma amostra pode ser usada para aplicações a jusante depende do organismo alvo e do objetivo do estudo. Por exemplo, para capturar a expressão da maioria dos genes de plantas, em nossa experiência, são necessárias 30 milhões a 50 milhões de leituras, enquanto para fungos, apenas 10 milhões podem ser suficientes. Assim, o Rup não definirá limites de qualidade para exclusão de amostras problemáticas, mas fornecerá métricas para identificar diversos problemas que possam ser encontrados.

A segunda etapa (qualidade do mapeamento) avalia a precisão e a confiabilidade do alinhamento de leitura com o genoma de referência, especificando sua adequação para análises posteriores. Nem todas as leituras sequenciadas podem ser usadas para computação de abundância de genes; leituras que não se alinham ao genoma ou transcriptoma ou leituras que mapeiam para vários locais no genoma são ignoradas na maioria dos casos¹⁹ e não contribuem para a contagem geral de leituras. Os resultados da segunda etapa do pipeline podem ser usados para entender por que as leituras não são usadas na computação de abundância. Um alto número de leituras não mapeadas pode indicar contaminação durante a extração de RNA (por exemplo, com patógenos de plantas ou herbívoros) ou um genoma de referência incompleto. Modelos genéticos incompletos do genoma de referência podem resultar em um grande número de leituras "sem recurso". Uma grande fração de leituras multimapeadas pode ser causada por um alto número de genes de rRNA na biblioteca, indicando remoção insuficiente de rRNA durante o processo de preparação da biblioteca de sequenciamento (no caso de as leituras multimapeadas serem realmente derivadas de rRNAs, elas podem ser analisadas posteriormente fornecendo um arquivo de anotação de rRNA para o pipeline, que calcula automaticamente o número de leituras de rRNA possíveis).

A terceira métrica de qualidade, igualmente importante, é a correlação entre réplicas da mesma amostra. Geralmente, a correlação entre réplicas de uma mesma amostra deve ser maior do que a correlação entre réplicas de amostras diferentes. Uma baixa correlação entre as réplicas pode indicar grande variabilidade biológica entre as réplicas, alta similaridade entre as amostras/condições/tecidos, outros efeitos de lote ou até mesmo troca de amostras ou rotulagem incorreta. O Rup calcula as correlações em pares entre todas as amostras e produz um mapa de calor agrupado das correlações da amostra. Além disso, uma análise de componentes principais (PCA) pode ser calculada nas amostras para identificar possíveis efeitos de lote que precisam ser reconhecidos em análises posteriores a jusante. Embora os efeitos do lote possam ser corrigidos em análises posteriores, pode ser melhor remover amostras problemáticas se a diferença medida for muito alta.

O material de entrada tem um impacto direto na qualidade do sequenciamento: amostragem de tecido, extração de RNA e preparação de biblioteca são etapas críticas para minimizar os problemas de qualidade do RNA-seq. Condições consistentes devem ser mantidas, se possível, por exemplo, usando câmaras de crescimento. As medidas de controle de pragas devem ser realizadas em tempo hábil, pois apenas indivíduos saudáveis devem ser selecionados para amostragem. Além disso, é aconselhável coletar amostras no mesmo dia e/ou ao mesmo tempo para reduzir a variação do transcriptoma circadiano. Numerosos kits de extração de RNA estão disponíveis, e a seleção de um kit apropriado para as espécies-alvo pode melhorar a qualidade do mRNA. Os protocolos de preparação da biblioteca devem incluir etapas para o enriquecimento de RNA poliA + para minimizar a fração de rRNA.

O Rup pode ser usado como uma etapa inicial de controle de qualidade na análise de expressão gênica diferencial. Esse controle de qualidade é necessário por vários motivos críticos, pois garante a qualidade dos dados de entrada (identifica réplicas/amostras de baixa qualidade, pode definir limites de qualidade para profundidade de leitura e taxas de mapeamento) e identifica problemas técnicos, como erros de sequenciamento, efeitos de lote ou rotulagem incorreta. Se a entrada de RNA for de qualidade suficiente, o Rup pode ajudar cortando leituras de baixa qualidade ou identificando amostras rotuladas incorretamente. No entanto, o Rup não pode compensar a baixa qualidade do RNA de entrada. Em casos de RNA de baixa qualidade ou dados de sequenciamento errôneos, pode ser necessário coletar novamente amostras, ajustar o protocolo de extração de RNA e/ou repetir o sequenciamento. O mesmo se aplica a amostras com uma grande fração de leituras não mapeadas que podem ter sido causadas por contaminação. No entanto, de acordo com nossa experiência, nem todas as extrações de RNA podem ser simplesmente repetidas devido à falta de disponibilidade de material. Nesse caso, eles são inadequados para algumas aplicações downstream: amostras com baixo número de leituras mapeadas únicas não devem ser analisadas em análises de expressão gênica diferencial, mas ainda contêm informações para a análise da presença de transcritos. Nesse caso, a ausência de transcritos em tecidos/tratamentos/condições não pode ser considerada para a análise e quantificação da abundância de transcritos.

O Rup inclui algumas limitações. Por exemplo, ele não testa diretamente a degradação do RNA, pois são necessárias medições diretas da integridade do RNA antes da preparação e sequenciamento da biblioteca. No entanto, um script para identificar a degradação de RNA usando os módulos RSeQC geneBodyCoverage.py e tin.py está incluído no repositório desse pipeline. Além disso, o Rup não testa o conteúdo do GC e o viés de comprimento da transcrição. O tamanho do genoma de referência está atualmente limitado a 4 Gb e, de acordo com nossa experiência, o módulo de mapeamento de leitura superestima as leituras multimapeadas em poliplóides, conforme exemplificado no Arquivo Suplementar em uma comparação do mapeamento de leitura para uma versão haploide vs. as versões do genoma diplóide de E. californica, um genoma haplóide é, portanto, a entrada preferida para este pipeline. A qualidade de saída do Rup depende em grande parte da qualidade do genoma de referência e da anotação, de modo que uma sequência do genoma incompleta ou altamente fragmentada pode levar a uma superestimação das leituras não mapeadas. Além disso, a anotação incompleta do modelo genético leva a uma superestimação de leituras não atribuídas. A completude e a duplicação da sequência do genoma e da anotação do gene podem ser inferidas com ferramentas como o BUSCO²⁰.

O Rup é uma ferramenta autônoma para todas as etapas iniciais de controle de qualidade essenciais em experimentos de RNA-seq. Ao contrário do RSeQC e do RNA-SeQC, o pré-processamento das leituras brutas de sequenciamento não é necessário. A análise da qualidade do sequenciamento não pode ser feita com RNA-SeQC, e os mapas de calor de correlação de amostra não são calculados por RSeQC e RNA-SeQC (Tabela 1). Além disso, a saída de normalização está em FPKM no RSeQC e a visualização de saída não é implementada como padrão para todos os módulos no RSeQC e RNA-SeQC. Assim, vários problemas de controle de qualidade não são testados nas duas ferramentas alternativas, que incluem baixo número de leituras, alta fração de leituras aparadas e a identificação de outliers replicados. Uma comparação de Rup, RSeQC e RNA-SeQC está disponível no Arquivo Suplementar 1. Além disso, o Rup pode ser usado de forma complementar ao RNA-SeQC ou RSeQC, por exemplo, os arquivos BAM ordenados produzidos por este pipeline podem ser usados como entrada para essas ferramentas.

Atualmente, o Rup é otimizado para um pequeno número de amostras, mas as etapas mais demoradas, como corte de qualidade e mapeamento de leitura, podem ser pré-computadas, por exemplo, em um cluster de computadores ou em uma infraestrutura de nuvem. Para genomas maiores que 4 Gb, isso é obrigatório, uma vez que o alinhador Rsubread é limitado a genomas menores que 4 Gb. A anotação de rRNA é feita com barrnap, que identifica os genes de rRNA altamente conservados. De acordo com nossa experiência, a anotação do gene rRNA não requer abrangência, porque mesmo que alguns genes anormais de rRNA estejam faltando, uma visão geral bruta da presença de rRNA no conjunto de dados é suficiente para a avaliação da qualidade do conjunto de dados. Versões futuras do Rup podem mudar para um método de mapeamento diferente, como a ferramenta de pseudo-alinhamento salmon²¹, para diminuir o tempo de execução. Além disso, mais métodos de normalização e correção, como viés de GC ou correção de viés de comprimento, podem ser adicionados ao Rup.

Em resumo, o Rup fornece informações essenciais para garantir a confiabilidade e reprodutibilidade da análise de dados de RNA-seq. Este pipeline relata de forma abrangente as principais métricas de qualidade de RNA-seq e produz arquivos de saída para uso direto em análises downstream. Ele foi projetado como uma ferramenta autônoma para pesquisadores com conhecimento mínimo de bioinformática para avaliar a qualidade de seus dados de sequenciamento primário com visualização intuitiva.

Materiais

Lista de materiais utilizados neste artigo

Nome	Empresa	Número de catálogo	Comentários
fastQCR	R	0.1.3
FUJITSU ESPRIMO D958 Desktop	FUJITSU		o pipeline foi desenvolvido e testado no Ubuntu Linux 24.02 (32 GB de RAM)
getopt	R	1.20.4
ggplot2	R	3.5.2
MacBook Pro & nbsp;	Maçã		o pipeline foi testado no maxOS 15.4.1 (16 GB de RAM)
pheatmap	R	1.0.13
R		4.4.3
Reshape2	R	1.4.4
RFASP	Biocondutor	1.16.0
rsamtools	Biocondutor	2.22.0
rsubread	Biocondutor	2.20.0