Rup(RNA-seq 사용성 평가 파이프라인) - 진핵생물의 대량 RNA-seq 실험을 위한 품질 관리

요약

이 프로토콜은 생물정보학 경험이 제한된 습식 실험실 생물학자를 위한 RNA-seq 실험에 대한 초기 품질 관리를 허용합니다.

초록

분자 식물 과학의 현대적 접근 방식에는 예를 들어 처리 시 전사체의 전반적인 변화를 추적하거나 조절 경로의 주요 구성 요소를 식별하기 위해 대량 RNA-seq 실험이 필요한 경우가 많습니다. 결과적으로 식물 과학의 다양한 분야는 과학적 발전을 위해 고품질의 재현 가능한 벌크 RNA-seq 데이터에 의존하고 있습니다. 그러나 우리의 경험에 비추어 볼 때 RNA-seq 데이터 세트의 품질 관리 조치에 대한 지식과 적용은 종종 부족합니다. 여기에서는 Rup에 대한 기본 지식을 가진 습식 실험실 생물학자에게 독립형으로 쉽게 적용할 수 있는 후속 유전자 발현 분석을 위한 Rup(RNA-seq 사용성 평가 파이프라인)을 소개합니다. Rup은 다운스트림 유전자 발현 실험에 적합한 고품질 시퀀싱 데이터와 일반적인 추가 분석에 적합하지 않은 데이터를 구별하는 데 도움이 됩니다. Rup에는 불충분한 판독 수 또는 매핑, 오염 식별, 전체 RNA-seq 데이터에서 rRNA 분획의 정량화, 복제 유사성 테스트, 데모를 위한 실제 데이터 사용 및 직관적인 시각화 제공과 같이 일반적으로 발생하는 몇 가지 문제에 대한 테스트가 포함됩니다. Rup은 표준화된 전사체 분석 전에 실험적 단점을 식별할 수 있는 도구 모음을 제공하여 개별 연구자와 현장의 데이터 품질을 향상시킵니다. 이는 대량 RNA-seq 데이터 분석에 대한 신뢰도를 높이고 최소 품질 관리 기준을 정의하는 향후 지침의 기반을 제공하여 게시된 RNA-seq 데이터의 신뢰성과 투명성을 향상시킵니다. Rup 및 테스트 데이터는 https://github.com/oliverrupp/rup 에서 확인할 수 있습니다.

서론

전사체학(RNA-seq) 실험은 표현형을 형성하는 전사 시그니처를 종합적으로 조사합니다. 이 접근법은 예를 들어 발달, 식물 병원체 상호 작용 또는 비생물적 스트레스 저항성과 관련된 단일 또는 공동 발현 유전자 및 생물학적 과정을 식별하기 위해 식물 분자 유전학에서 대체할 수 없게 되었습니다^{1, 2 ,3}. RNA-seq 기술의 최근 발전으로 특이성이 향상되고 단일 염기 분해능에서 다양한 이소형 및 변이체를 검출할 수 있게 되어 더 큰 인델에서 단일 뉴클레오티드 다형성(SNP)에 이르는 서열 변이를 식별할 수 있게 되었습니다. RNA-seq로 얻은 데이터는 넓은 동적 범위를 특징으로 하여 매우 풍부하고 저발현된 전사체를 모두 쉽게 검출할 수 있으며 일관성을 위해 적절한 실험 조건이 필요합니다. 또한 RNA-seq 데이터 세트는 규모가 커서 분석 및 저장에 계산 집약적이므로 효율적인 데이터 관리와 광범위한 컴퓨팅 리소스가 필요합니다. RNA-seq는 모든 단계의 약점이 결과를 전파하고 손상시킬 수 있으므로 워크플로의 모든 단계에서 고품질 입력이 필요합니다. 라이브러리 준비 중 RNA 무결성이 좋지 않거나 기술적 문제가 발생하면 편향이 발생하고 정확도가 감소합니다. 고품질 원시 판독을 얻기 위한 매개변수는 차세대 염기서열 분석(NGS) 시설마다 다를 수 있습니다.

우리의 경험에 따르면, RNA-seq 라이브러리 준비를 위한 입력으로 RNA 무결성 수(RIN)가 7 이상인 RNA만 사용하는 것이 좋습니다. 성공적인 시퀀싱을 위해서는 표준 라이브러리 준비 프로토콜에 50–200 ng/μL 농도의 총 RNA 약 2 μg이 필요합니다. RNA 순도는 분광 광도계를 사용하여 1.8에서 2.1 사이의 OD_260/280 비율과 1.5보다 큰 OD_260/230 비율로 확인해야 합니다. 불충분한 시퀀싱 깊이, 낮은 매핑률 또는 참조 게놈에 대한 정렬 불량은 유전자 발현 및 스플라이싱 정량화를 더욱 왜곡할 수 있습니다. 더욱이, 서로 다른 조직이나 처리의 샘플 간의 차별 부족, 반복 간의 높은 변동성과 같은 부적절한 실험 설계는 노이즈를 유발하고 재현성을 감소시킬 수 있습니다. 많은 실험실에서 전사체를 일상적으로 분석하지만 첫 번째 필수 분석 단계에 대한 엄격한 품질 관리는 종종 보고되지 않거나 간행물에 전혀 없을 수 있습니다. 이로 인해 전사체 분석에서 파생된 결과가 과도하게 해석되어 결과적으로 재현할 수 없는 결과가 발생할 수 있습니다.

여기에서는 전사 변경을 측정하기 위해 mRNA 프로파일의 고품질 전사체 분석에 필요한 초기 단계의 품질 관리를 위한 워크플로를 제공합니다. 우리의 목표는 생물정보학에 대한 지식이 제한된 습식 실험실 생물학자가 1차 전사체학 데이터를 평가할 수 있도록 하는 것입니다. Rup(그림 1)은 R에 대한 기본 지식에 익숙한 연구자가 액세스할 수 있습니다. 여기에 제시된 워크플로를 실행하면 연구자는 후속 분석에 대한 잠재적인 한계를 포함하여 기본 데이터에 대한 자세한 이해를 얻을 수 있습니다. 우리가 아는 한, 고품질 데이터와 저품질 데이터를 구별하기 위한 지침과 결합된 실용적인 1차 전사체 데이터 평가 파이프라인은 지금까지 부족합니다.

그림 1: RNA-seq in silico 품질 관리 파이프라인의 워크플로. 품질 관리를 위한 입력 파일은 식물 재료 시퀀싱을 통해 생성된 RNA-seq 데이터와 공개적으로 사용 가능한 데이터 세트에서 파생됩니다. 제공된 R 파이프라인은 시퀀싱 품질, 매핑 품질 및 복제 품질의 세 가지 주요 접근 방식을 통해 시퀀스 품질을 평가합니다. 다양한 품질 메트릭이 계산되고 ggplot2 및 pheatmap과 같은 일반적인 R 패키지(예: 막대 차트 및 히트맵)를 사용하여 통계 결과를 시각화합니다. 이 그림의 더 큰 버전을 보려면 여기를 클릭하십시오.

이전에 보고된 평가 파이프라인은 사전 처리된 데이터(예: .bam 파일로 정렬 읽기)가 필요하거나, 모든 메트릭을 다루지 않거나, 더 이상 유지되지 않습니다 ^4,5,6. 여기에 제시된 워크플로의 장점은 가장 일반적인 품질 관리 문제를 단일 파이프라인으로 통합함으로써 포괄적이라는 것입니다. 또한 모든 다운스트림 분석 애플리케이션에 적합한 고품질 데이터에 대한 예와 저품질 데이터에 대한 예도 제공하고 추가 분석을 위한 특정 제한 사항에 대해 논의합니다. 이전에 발표된 여러 보고서는 RNA-seq 분석 도구의 목적을 설명하고 성능에 대한 비교 평가를 제공합니다 ^7,8. 그러나 RNA-seq 품질 관리 및 방법론 보고 표준은 표준화되지 않았으며 전사체학 실험의 재현성과 생물학적으로 의미 있는 해석을 악화시킵니다.

Rup은 표준 고품질 도구, 품질 관리 분석 및 결과 시각화를 통합합니다. Rup은 원시 시퀀싱 읽기와 주석이 달린 게놈을 입력으로 필요로 하며 Mac OS, Linux 및 Windows 시스템의 "Linux용 Windows 하위 시스템"(WSL)에서 실행됩니다. 시퀀싱 품질을 위한 도구를 통합하고 게놈에 대한 읽기 매핑을 정량화하고, 샘플에 rRNA 함량 측정을 포함하고, 복제 상관 관계를 평가하기 위해 샘플 상관 관계를 제공합니다. 테스트 데이터는 라이브러리 준비를 위한 초저 입력 프로토콜인 식물 종 Eschscholzia californica의 두 가지 다른 단계의 중앙 분열 조직 조직의 레이저 미세 해부를 사용하여 얻었고 Novaseq 6000에서 시퀀싱되었습니다.

프로토콜

Rup은 최소한의 입력 파일이 필요한 단일 스크립트로 실행할 수 있습니다. 여기에 표시된 대로 파이프라인에는 Illumina 쌍 말단 RNA 시퀀싱 데이터가 필요합니다. 동일한 분석 단계를 사용하는 단일 엔드 시퀀싱을 위해 조정된 파이프라인도 GitHub 리포지토리에 보관됩니다. fasta 파일로서의 게놈 서열, gtf 파일로서의 유전자 모델 및 rRNA 주석, fastq.gz 파일로서의 원시 시퀀싱 판독만 필요합니다. 게놈 및 주석 파일이 reference_folder 변수에 지정된 폴더에 genome.fa, annotation.gtf 및 rRNA.gtf로 제공되는지 확인합니다. 모든 시퀀싱 파일을 .fq.gz 파일로 read_file_folder 변수에 저장하면 다음과 같이 Rup을 실행할 수 있습니다.

1. 준비:

1. Bioconductor를 사용하여 필요한 모든 R 패키지를 설치합니다.

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

2. 필요한 모든 파일을 설정합니다.
   1. 필요한 모든 파일을 포함할 소스 폴더를 만듭니다. "reference/genome.fa"라는 fasta 파일의 참조 게놈 서열을 소스 폴더에 추가합니다. 유전자 모델 주석을 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))

   
   
   참고: 최소 읽기 길이는 25bp보다 작아서는 안 되며, 더 작은 읽기를 하면 매핑이 중단될 수 있습니다. 값이 클수록 더 고유하게 매핑된 읽기가 생성되지만 시퀀싱 품질에 따라 트리밍 중에 더 많은 읽기가 삭제됩니다. 이러한 매개변수에 대한 자세한 내용은 보충 자료에 설명되어 있습니다.

2. 시퀀싱 품질 평가

참고: 이 단계는 각 샘플을 트리밍하기 전과 후의 판독 수를 보여주는 막대 플롯을 생성합니다.

1. 고처리량 염기서열 데이터의 일반적인 품질 관리를 위한 도구인 fastqc⁹를 실행하여 총 염기서열 판독 수, 판독 길이, 평균 GC 함량, 염기서열 분석 어댑터 오염 및 과도하게 표현된 염기서열(예: rRNA 판독)에 대한 정보를 제공합니다. 염기 및 염기서열 품질 보고서별로 염기서열 분석 품질을 분석합니다. R 패키지 fastqcr¹⁰ 을 사용하여 원시 시퀀싱 파일에서 fastqc를 실행하고 모든 입력 파일을 동일한 디렉토리(read_file_folder)에 배치합니다. 함수 qc_aggregate를 사용하여 fastqc_folder의 출력 파일을 R 개체로 요약합니다.
   참고: 결과는 디렉토리 fastqc_folder에 저장됩니다. 출력 폴더에서 fastqc가 만든 .html 파일에 액세스하여 추가 통계를 얻을 수 있습니다.
    

   # 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를 실행합니다. 어댑터 서열, 프라이머, 저품질 염기, 폴리-A/T 서열을 제거하기 위한 트림은 R 패키지 Rfastp¹²에서 사용할 수 있는 Fastp¹¹에 의해 실행됩니다. 트리밍된 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. 매핑 품질

참고: 이 단계는 막대 플롯을 계산하여 단일 매핑된 판독(예: 단백질 코딩 유전자의 전사체)과 다중 매핑된 판독(예: rRNA 판독) 및 매핑되지 않은 판독(예: 오염)을 구별합니다.

1. Rsubread¹⁴ 패키지를 사용하여 판독값을 참조 게놈에 매핑합니다. 먼저 참조 게놈 fasta 파일(genome_fasta_file)의 인덱스를 구축합니다.
   참고: 이 단계는 각 참조 게놈에 대해 한 번 수행됩니다.
    

   # load the Rsubread library for read mapping and read counting
library(Rsubread)

# create a subread index from the reference genome sequence
buildindex(file.path(reference_folder,"subread.index"), genome_fasta_file)

2. 모든 샘플을 반복하고 Rsubread 패키지의 align() 함수를 사용하여 참조 게놈에 읽기를 정렬합니다. 이 함수는 출력 폴더에 .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. Rsubread 패키지의 featureCounts() 함수를 사용하여 각 유전자에 대한 판독 횟수를 계산합니다. 주석 파일(annotation_file)이 GTF 형식인지 확인합니다. 게놈과 단일 일치가 있는 판독만 계산합니다.
    

   # collect all bam files
bam_files <- list.files(bam_folder, pattern = "*.bam$")

# count the number of reads for each gene
     gene_feature_counts = featureCounts(file.path(bam_folder, bam_files),
     annot.ext = annotation_file, # the gene models
     isGTFAnnotationFile = TRUE,
     countMultiMappingReads = FALSE, # only count unique reads
     strandSpecific = stranded, # for strand specific data
     isPairedEnd = TRUE, # paired-end sequencing
# the following paramter allow to control for false-positive read assignments
     requireBothEndsMapped = TRUE,
     checkFragLength = TRUE,
     minFragLength = minFragLength,
     maxFragLength = maxFragLength,
     nthreads = n_threads)

4. 샘플의 rRNA 함량을 평가하기 위해 rRNA 유전자에 매핑되는 판독 수를 계산합니다. 여기에서 다중 매핑된 읽기를 계산할 수 있습니다. GTF 파일(rrna_file)에서 rRNA 유전자 유전자좌를 지정합니다.
    

   # 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. rRNA 할당의 통계를 수집합니다.
    

   # collect the rrna assignment statistics
rrna_count_stats <- rrna_feature_counts$stat
# set the sample names as column names
colnames(rrna_count_stats) <- gsub(".bam", "", colnames(rrna_count_stats))

# only the number of assigned reads are important in this case
rrna_counts = as.data.frame(t(rrna_count_stats[rrna_count_stats$Status == "Assigned",2:ncol(rrna_count_stats)]))
colnames(rrna_counts) = "Alignments"

# setup names and categories for plotting
rrna_counts$Sample = rownames(rrna_counts)
rrna_counts$Group = "rRNA"
rrna_counts$Reference = "rRNA"

7. 읽기 매핑 통계량을 막대 그림으로 표시합니다.
    

   # combine the proteind coding and rRNA gene statistics
mapping_stats = rbind(transformed_stats, rrna_counts)

# plot the assignment statistics
mapping_stats_plot = ggplot(data = mapping_stats, aes(x = Sample, y = Alignments)) +
 geom_col(aes(fill = Group), width = 0.7) +
 theme_bw() + facet_wrap(~Reference) +
 ylab("Number of Alignments") + xlab("Samples") +
 ggtitle("Read Mapping Numbers") +
 theme(text = element_text(size = 18)) +
 theme(axis.text.x = element_text(angle = 90, hjust = 0))
print(mapping_stats_plot)

8. 유전자에 할당된 판독 수에 따라 유전자를 그룹으로 분류하고 결과를 막대 그림으로 표시합니다.
    

   # get the reads per gene counts
gene_count_matrix = gene_feature_counts$counts
# set the sample names as column names
colnames(gene_count_matrix) = gsub(".bam", "", colnames(gene_count_matrix))

# group and count genes in classes of specific read counts
read_count_classes = data.frame("no_reads"=colSums(gene_count_matrix == 0),
     "at_least_1_read"=colSums(gene_count_matrix >= 1 & gene_count_matrix < 10),
     "at_least_10_reads"=colSums(gene_count_matrix >= 10 & gene_count_matrix < 100),
     "at_least_100_reads"=colSums(gene_count_matrix >= 100 & gene_count_matrix < 1000),
     "at_least_1000_reads"=colSums(gene_count_matrix >= 1000))

# setup data.frame for plotting
read_count_classes$Sample = rownames(read_count_classes)
melt_rcc = melt(read_count_classes)
melt_rcc$variable = as.character(melt_rcc$variable)

# sort the gene groups
melt_rcc$variable = factor(melt_rcc$variable, levels=c("no_reads",
         "at_least_1_read",
         "at_least_10_reads",
         "at_least_100_reads",
         "at_least_1000_reads"))

# plot the data as bar plot
gene_coverage_plot <- ggplot(melt_rcc, aes(x=Sample, y=value, fill=variable)) +
geom_bar(stat="identity") +
ylab("Number of Genes") + xlab("Samples") +
ggtitle("Number of Reads per Gene") +
guides(fill=guide_legend(title="Number of assigned reads")) +
theme(text = element_text(size = 18)) +
theme(axis.text.x = element_text(angle = 90, hjust = 0))
print(gene_coverage_plot)

그림 3: 매핑 개요 및 통계를 읽습니다. 대부분의 샘플은 할당된 판독 횟수가 많고(왼쪽, 갈색) rRNA 판독 횟수가 적습니다(오른쪽, 파란색). 샘플 s2_r3에는 높은 rRNA 판독 수(오른쪽)에 해당하는 비정상적으로 많은 양의 다중 매핑 판독(왼쪽, 녹색)이 있습니다. 샘플 s2_r4는 예상 rRNA 판독 수와 함께 매핑되지 않은 판독 수가 많은지 보여주며, 이는 다른 유기체의 판독으로 오염되었음을 암시합니다. 이 그림의 더 큰 버전을 보려면 여기를 클릭하십시오.

4. 품질을 복제합니다.

1. 동일한 샘플 또는 조직의 복제본이 다른 샘플 또는 조직의 복제본보다 서로 더 높은 상관관계를 나타내는지 확인하십시오. 복제 상관 관계의 클러스터형 히트맵을 플로팅하여 원시 읽기 횟수 또는 TPM(Transcripts Per Million)¹⁵ 정규화된 횟수를 사용하여 데이터를 시각화합니다.
    

   # 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: 유전자 판독 수 분류. 샘플의 모든 유전자는 할당된 판독 횟수에 따라 5가지 범주 중 하나로 분류됩니다. 빨간색으로 표시된 것은 판독이 할당되지 않은 유전자의 수입니다. 총 할당된 리드 수가 적은(s2_r1 - s1_r4) 10-100개의 할당된 리드를 가진 더 많은 수의 유전자를 가지고, 1000개 이상의 리드를 가진 더 적은 수의 유전자를 갖습니다. 발현이 낮은 유전자는 이러한 샘플에서 누락될 수 있습니다. 이 그림의 더 큰 버전을 보려면 여기를 클릭하십시오.

결과

파이프라인은 R 스크립트로 완전히 구현되며 Linux 및 Mac OS 운영 체제에서 테스트되었습니다. Windows 사용자는 WSL(Linux용 Windows 하위 시스템)을 사용할 수 있습니다. 코드 및 테스트 데이터는 GitHub 리포지토리(https://github.com/oliverrupp/rup)로 사용할 수 있습니다. 시퀀싱 데이터는 ENA EBI 프로젝트 PRJEB96400에서 사용할 수 있습니다.

분석을 위해 Rup을 사용하여 벌크 RNA-seq의 품질 관리 중에 발생할 수 있는 다양한 문제를 예시하기 위해 두 개의 실제 샘플에서 10개의 샘플을 인위적으로 만들었습니다. 샘플 s2_r1은 낮은 총 판독 수를 나타내도록 설계되었으며, 샘플 s2_r2에는 트리밍 프로세스에 의해 폐기될 저품질 판독의 상당 부분이 포함되어 있습니다. 샘플 s2_r3에는 rRNA 판독의 상당 부분이 포함되며 샘플 s2_r4에는 참조 게놈에 매핑되지 못한 오염 물질 판독이 포함됩니다. 낮은 반복 상관 관계를 설명하기 위해 s1_r5와 s2_r5의 샘플 이름을 바꿨습니다.

프로토콜 섹션 2를 사용하면 트리밍 전후에 판독 수가 낮은 샘플을 식별할 수 있습니다(그림 2). 막대 그림은 트리밍 전후의 샘플 s1_r1에서 더 낮은 판독 수를 보여줍니다. 여기서 초기 읽기 수는 낮았습니다. 트리밍 후 샘플 s1_r2의 판독 수가 적다는 것은 트리밍 공정 중에 제거된 많은 양의 어댑터 서열, 시퀀싱 오류, 프라이머 서열, 시퀀싱된 폴리-A/T 스트레치를 시사합니다. 입력 RNA의 분해는 고품질 판독 횟수를 감소시킬 수도 있습니다.

프로토콜 섹션 3은 읽기 할당의 문제를 식별합니다. 그림 3 은 샘플 s2_r3(녹색)에서 증가된 다중 매핑 판독 수와 높은 수의 rRNA 판독을 보여줍니다. 샘플 s2_r4의 오염은 참조 게놈(분홍색)에 매핑되지 않은 판독의 상당 부분에 의해 명백합니다. 이러한 판독은 rRNA와 상관관계가 없지만 비표적 유기체의 서열과 상관관계가 있습니다. 그림 4 는 전사체에서 할당된 판독 수가 적은 것과 관련된 일반적인 문제를 보여줍니다. 참조 게놈(s2_r1에서 r4까지)에 고유하게 매핑된 판독률이 낮은 샘플에서는 유전자의 약 3분의 1만이 100개 이상의 판독이 할당된 반면, 다른 샘플에서는 유전자의 약 절반이 이러한 클래스에 속합니다. 발현이 매우 낮은 유전자의 판독은 r4 s2_r1 샘플에서 발견되지 않을 수 있으므로 이러한 샘플과의 비교 발현 분석은 매우 신뢰할 수 없으며 피해야 합니다.

프로토콜 섹션 4는 복제 이상값을 식별하는 데 사용할 수 있습니다. 동일한 샘플/조건/조직의 복제는 다른 샘플/조건/조직의 복제보다 서로 더 높은 상관관계를 나타낼 것으로 예상됩니다. 그림 5 는 각각 5개의 반복실험이 있는 두 샘플(S1 및 S2)의 상관 히트맵을 보여줍니다. 플롯의 상단과 측면에 있는 덴드로그램은 각 클러스터에 5개의 반복이 있는 두 개의 클러스터를 보여줍니다. 왼쪽 클러스터에는 샘플 1의 반복실험 4개와 샘플 2의 반복실험 1개(s2_r5)가 포함되고, 오른쪽 클러스터에는 샘플 2의 반복실험 4개와 샘플 1의 반복실험 1개(s1_r5)가 포함되어 있습니다. 이 경우 샘플 이름 s1_r5와 s2_r5가 다시 교체되면 각 클러스터에는 한 샘플의 모든 반복이 포함되며, 이는 반복 레이블 지정 오류를 나타낼 수 있습니다. 복제물이 함께 클러스터링되지 않는 다른 이유는 샘플/조건/조직 간의 차별화 부족 또는 시퀀싱 품질상의 이유로만 복제의 클러스터링일 수 있습니다. 후자는 트리밍 및 매핑 후 읽기 수가 매우 낮은 모든 복제 또는 읽기 수가 매우 높은 복제가 클러스터를 형성할 때 발생할 수 있습니다.

그림 5: 샘플 상관 관계 히트맵. 샘플 상관 관계 히트맵은 로그₂ 변환된 TPM 값을 기반으로 하며 각각 5개의 샘플로 구성된 두 개의 개별 클러스터를 보여줍니다. 생물학적/기술적 복제는 다른 조직/치료의 복제보다 서로 더 높은 상관관계를 보일 것으로 예상됩니다. 왼쪽 클러스터에는 샘플 1의 반복실험 4개와 샘플 2의 반복실험 1개(s2_r5)가 포함되고, 오른쪽 클러스터에는 샘플 2의 반복실험 4개와 샘플 1의 반복실험 1개(s1_r5)가 포함되어 있습니다. 단일 반복실험은 히트맵에서 클러스터링을 설명하기 위해 가능한 샘플 스와핑 또는 배치 효과에 대해 확인해야 합니다. 이 그림의 더 큰 버전을 보려면 여기를 클릭하십시오.

표 1: RNA-seq 품질 평가 파이프라인의 비교. 자세한 분석은 보충 파일 1을 참조하십시오. 이 표를 다운로드하려면 여기를 클릭하십시오.

보충 파일 1: 매개변수 및 참조 선택. 분석은 매개변수 및 참조 선택이 전체 QC 결과에 미치는 영향을 보여줍니다. 이 파일을 다운로드하려면 여기를 클릭하십시오.

토론

차등 유전자 발현 분석의 품질은 두 가지 요인, 즉 각 샘플에서 시퀀싱된 판독 횟수 및 복제^{횟수 16} 및 샘플당 반복 횟수^17,18에 크게 좌우됩니다. 여기에서는 각 반복에서 유전자 발현 정량화에 적합한 판독 수를 결정하기 위해 사용자 친화적인 Rup 파이프라인을 제시합니다. 다양한 메트릭을 통해 연구자는 반복실험에서 할당된 판독 수가 낮은 이유를 이해하고 샘플 상관 관계의 문제를 식별할 수 있습니다. Rup은 식물 RNA-seq 샘플의 품질 평가를 위해 개발되었지만 다른 진핵 생물에도 동일하게 적합하므로 비식물 샘플에 대한 추가 조정이 필요하지 않습니다(보충 파일 1).

Rup의 첫 번째 단계는 읽기 트리밍 전후의 총 읽기 수를 결정합니다. 시퀀싱된 판독 횟수에 따라 검출 가능한 유전자의 수가 결정되며, 판독 횟수가 너무 낮으면 차등적으로 발현되는 많은 유전자가 검출되지 않은 상태로 남아 있습니다. 품질 트리밍 및 필터링 과정에서 폐기된 판독의 상당 부분은 RNA 분해를 나타낼 수 있으며, 이는 고도로 분해된 유전자에 대한 발현을 과소평가할 수 있습니다. 그러나 샘플을 다운스트림 응용 분야에 사용할 수 있는지 여부는 표적 유기체와 연구 목적에 따라 다릅니다. 예를 들어, 대부분의 식물 유전자의 발현을 포착하려면 우리의 경험에 따르면 3천만 내지 5천만 개의 읽기가 필요한 반면, 곰팡이의 경우 1천만 개만 있으면 충분할 수 있습니다. 따라서 Rup은 문제가 있는 샘플을 제외하기 위한 품질 임계값을 정의하지 않고 발생할 수 있는 다양한 문제를 식별하기 위한 메트릭을 제공합니다.

두 번째 단계(매핑 품질)는 참조 게놈에 대한 읽기 정렬의 정확성과 신뢰성을 평가하여 다운스트림 분석에 대한 적합성을 지정합니다. 모든 시퀀싱된 판독이 유전자 풍부도 계산에 사용될 수 있는 것은 아닙니다. 게놈 또는 전사체에 정렬되지 않는 판독 또는 게놈의 여러 위치에 매핑되는 판독은 대부분의 경우 무시되며¹⁹ 전체 판독 횟수에 기여하지 않습니다. 파이프라인의 두 번째 단계의 결과는 풍부도 계산에 읽기가 사용되지 않는 이유를 이해하는 데 사용할 수 있습니다. 매핑되지 않은 판독 수가 많으면 RNA 추출 중 오염(예: 식물 병원체 또는 초식 동물) 또는 불완전한 참조 게놈을 나타낼 수 있습니다. 참조 게놈의 불완전한 유전자 모델은 많은 수의 "특징 없음" 판독을 초래할 수 있습니다. 다중 매핑된 판독의 상당 부분은 시퀀싱 라이브러리 준비 과정에서 rRNA 제거가 불충분함을 나타내는 라이브러리의 많은 수의 rRNA 유전자로 인해 발생할 수 있습니다(다중 매핑된 판독이 실제로 rRNA에서 파생된 경우, 파이프라인에 rRNA 주석 파일을 제공하여 추가 분석할 수 있으며, 파이프라인은 가능한 rRNA 판독 수를 자동으로 계산합니다).

세 번째, 똑같이 중요한 품질 지표는 동일한 샘플의 반복 간의 상관 관계입니다. 일반적으로 동일한 샘플의 반복실험 간의 상관관계는 서로 다른 샘플의 반복실험 간의 상관관계보다 높아야 합니다. 반복 간의 낮은 상관관계는 반복 간의 큰 생물학적 변동성, 샘플/조건/조직 간의 높은 유사성, 기타 배치 효과 또는 샘플 교체 또는 잘못된 라벨링을 나타낼 수 있습니다. Rup은 모든 표본 간의 쌍별 상관관계를 계산하고 표본 상관관계의 군집형 히트맵을 생성합니다. 또한 샘플에 대해 주성분 분석(PCA)을 계산하여 추가 다운스트림 분석에서 인정해야 하는 가능한 배치 효과를 식별할 수 있습니다. 배치 효과는 다운스트림 분석에서 보정할 수 있지만 측정된 차이가 너무 높으면 문제가 있는 샘플을 제거하는 것이 더 나을 수 있습니다.

입력 물질은 시퀀싱 품질에 직접적인 영향을 미칩니다: 조직 샘플링, RNA 추출 및 라이브러리 준비는 RNA-seq 품질 문제를 최소화하는 데 중요한 단계입니다. 예를 들어 성장 챔버를 사용하여 가능하면 일관된 조건을 유지해야 합니다. 해충 방제 조치는 건강한 개인만 샘플링을 위해 선택해야 하므로 적시에 수행되어야 합니다. 일주기 전사체 변이를 줄이기 위해 같은 날 및/또는 동시에 샘플을 수집하는 것이 좋습니다. 수많은 RNA 추출 키트를 사용할 수 있으며 표적 종에 적합한 키트를 선택하면 mRNA 품질을 향상시킬 수 있습니다. 라이브러리 준비 프로토콜에는 rRNA 분획을 최소화하기 위한 polyA+ RNA 농축 단계가 포함되어야 합니다.

Rup은 차등 유전자 발현 분석의 초기 품질 관리 단계로 사용할 수 있습니다. 이러한 품질 관리는 입력 데이터의 품질을 보장하고(저품질 복제/샘플 식별, 판독 깊이 및 매핑 속도에 대한 품질 임계값 설정 가능) 시퀀싱 오류, 배치 효과 또는 잘못된 라벨링과 같은 기술적 문제를 식별하기 때문에 몇 가지 중요한 이유로 필요합니다. RNA 입력의 품질이 충분한 경우 Rup은 품질이 낮은 판독을 트리밍하거나 잘못 라벨링된 샘플을 식별하여 도움을 줄 수 있습니다. 그러나 Rup은 열악한 입력 RNA 품질을 보상할 수 없습니다. RNA 품질이 낮거나 잘못된 시퀀싱 데이터의 경우 샘플을 다시 수집하고, RNA 추출 프로토콜을 조정하고, 시퀀싱을 반복해야 할 수 있습니다. 오염으로 인해 발생했을 수 있는 매핑되지 않은 판독 중 많은 부분이 있는 샘플에도 동일하게 적용됩니다. 그러나 우리의 경험에 따르면 재료 가용성 부족으로 인해 모든 RNA 추출을 단순히 반복할 수 있는 것은 아닙니다. 이 경우 일부 다운스트림 응용 분야에는 적합하지 않습니다: 단일 매핑 판독 수가 적은 샘플은 차등 유전자 발현 분석에서 분석해서는 안 되지만 전사체 존재 분석을 위한 정보는 여전히 보유합니다. 이 경우, 조직/치료/조건에 전사체가 없다는 것은 전사체 풍부도의 분석 및 정량화를 위해 고려할 수 없습니다.

Rup에는 몇 가지 제한 사항이 있습니다. 예를 들어, 라이브러리 준비 및 시퀀싱 전에 RNA 무결성을 직접 측정해야 하므로 RNA 분해를 직접 테스트하지 않습니다. 그러나 RSeQC 모듈 geneBodyCoverage.py 및 tin.py 를 사용하여 RNA 분해를 식별하는 스크립트는 이 파이프라인의 저장소에 포함되어 있습니다. 또한 Rup은 GC 함량 및 전사체 길이 편향을 테스트하지 않습니다. 참조 게놈 크기는 현재 4Gb로 제한되어 있으며, 우리의 경험에 따르면, 읽기 매핑 모듈은 일배체에 대한 읽기 매핑과 E . californica의 이배체 게놈 버전을 비교하는 보충 파일에 예시된 바와 같이 배수체의 다중 매핑된 읽기를 과대평가하며, 따라서 반수체 게놈은 이 파이프라인에 선호되는 입력입니다. Rup 출력 품질은 참조 게놈 및 주석의 품질에 크게 좌우되므로 불완전하거나 고도로 단편화된 게놈 서열은 매핑되지 않은 읽기를 과대평가할 수 있습니다. 더욱이, 불완전한 유전자 모델 주석은 할당되지 않은 읽기의 과대 평가로 이어집니다. 게놈 서열 및 유전자 주석의 완전성과 중복은 BUSCO²⁰과 같은 도구를 사용하여 추론할 수 있습니다.

Rup은 RNA-seq 실험에 필수적인 모든 초기 품질 관리 단계를 위한 독립형 도구입니다. 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은 소수의 샘플에 최적화되어 있지만 품질 트리밍 및 읽기 매핑과 같은 가장 시간이 많이 걸리는 단계는 예를 들어 컴퓨터 클러스터 또는 클라우드 인프라에서 미리 계산할 수 있습니다. 4Gb보다 큰 게놈의 경우 Rsubread 정렬기가 4Gb보다 작은 게놈으로 제한되기 때문에 이는 필수입니다. rRNA 주석은 고도로 보존된 rRNA 유전자를 식별하는 barrnap으로 수행됩니다. 우리의 경험에 따르면 rRNA 유전자 주석은 일부 비정상적인 rRNA 유전자가 누락되더라도 데이터 세트의 rRNA 존재에 대한 총체적인 개요만으로도 데이터 세트 품질 평가에 충분하기 때문에 포괄성이 필요하지 않습니다. Rup의 향후 버전에서는 런타임을 줄이기 위해 의사 정렬 도구 salmon²¹과 같은 다른 매핑 방법으로 전환될 수 있습니다. 또한 GC 바이어스 또는 길이 바이어스 보정과 같은 더 많은 정규화 및 보정 방법을 Rup에 추가할 수 있습니다.

요약하면, Rup은 RNA-seq 데이터 분석의 신뢰성과 재현성을 보장하는 데 필수적인 정보를 제공합니다. 이 파이프라인은 주요 RNA-seq 품질 지표를 포괄적으로 보고하고 다운스트림 분석에 직접 사용할 출력 파일을 생성합니다. 최소한의 생물정보학 지식을 가진 연구자가 직관적인 시각화를 통해 기본 시퀀싱 데이터 품질을 평가할 수 있는 독립형 도구로 설계되었습니다.

재료

이 논문에 사용된 재료 목록

이름	회사	카탈로그 번호	댓글
fastQCR	R	0.1.3
후지쓰 ESPRIMO D958 데스크탑	후지쯔		파이프라인은 Ubuntu Linux 24.02(32GB RAM)에서 개발 및 테스트되었습니다
겟옵트	R	1.20.4
GGPLOT2	R	3.5.2
맥북 Pro	애플		파이프라인은 maxOS 15.4.1(16GB RAM)에서 테스트되었습니다
히트맵	R	1.0.13
R		4.4.3
리셰이프2	R	1.4.4
rfastp	생체도체	1.16.0
rsamtools	생체도체	2.22.0
rsubread	생체도체	2.20.0