Rup(RNA-seq Usability Assessment Pipeline) - 真核生物におけるバルクRNA-seq実験の品質管理

サマリー

このプロトコルにより、バイオインフォマティクスの経験が限られているウェットラボ生物学者のためのRNA-seq実験の初期品質管理が可能になります。

要約

分子植物科学における最新のアプローチでは、たとえば、処理時のトランスクリプトームの全体的な変化を追跡したり、調節経路の主要な構成要素を特定したりするために、バルクRNA-seq実験が必要になることがよくあります。その結果、植物科学のさまざまな分野が、科学の進歩のために高品質で再現性のあるバルク RNA-seq データに依存しています。しかし、私たちの経験から、RNA-seqデータセットにおける品質管理手段に関する知識と適用はしばしば不足しています。ここでは、Rupの基本的な知識を持つウェットラボ生物学者がスタンドアロンですぐに適用できる、バルクRNA-seqデータの品質管理、その後の遺伝子発現解析のための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の約2 μgのトータルRNAが必要です。RNA純度は、分光光度計を使用して、1.8〜2.1のOD_260/280 比と1.5を超えるOD_260/230 比によって確認する必要があります。配列決定の深さが不十分であったり、マッピング率が低いり、リファレンスゲノムへの位置がずれていたりすると、遺伝子発現やスプライシングの定量化がさらに歪む可能性があります。さらに、異なる組織や処理のサンプル間の識別の欠如や反復間のばらつきが大きいなど、実験計画が不十分であると、ノイズが発生し、再現性が低下する可能性があります。多くのラボがトランスクリプトームを日常的に分析していますが、最初の重要な分析ステップの厳格な品質管理は報告されていないか、出版物にまったくないことがよくあります。これにより、トランスクリプトーム解析から得られた結果が過剰に解釈され、その結果、再現性のない結果が生じる可能性があります。

ここでは、転写変化を測定するためのmRNAプロファイルの高品質のトランスクリプトーム解析に必要な初期ステップの品質管理のためのワークフローを提供します。私たちの目的は、バイオインフォマティクスの知識が限られているウェットラボ生物学者が、主要なトランスクリプトミクスデータを評価できるようにすることです。Rup(図1)は、Rの基本的な知識に精通している研究者がアクセスできます。ここで紹介するワークフローを実行すると、研究者は、その後の分析の潜在的な制限を含め、一次データを詳細に理解することができます。私たちの知る限り、高品質のデータと低品質のデータを区別するためのガイドラインと組み合わせた実用的な一次トランスクリプトームデータ評価パイプラインは、これまでのところ不足しています。

図1:RNA-seq in silico 品質管理パイプラインのワークフロー。 品質管理のための入力ファイルは、植物材料のシーケンシングによって生成されたRNA-seqデータ、および公開されているデータセットから派生しています。提供されているRパイプラインは、シーケンス品質、マッピング品質、および複製品質の3つの主要なアプローチを通じてシーケンス品質を評価します。さまざまな品質指標が計算され、ggplot2 や pheatmap などの一般的な R パッケージ (棒グラフやヒートマップなど) を使用して統計結果が視覚化されます。 この図の拡大版を表示するには、ここをクリックしてください。

以前に報告された評価パイプラインでは、前処理されたデータ (たとえば、.bam ファイルとしてのアライメントの読み取り) が必要であり、すべてのメトリックをカバーしていないか、もはや維持されていません ^4,5,6。ここで紹介するワークフローの利点は、最も一般的な品質管理の問題を単一のパイプラインに統合することによる包括性にあります。さらに、すべての下流分析アプリケーションに適した高品質のデータの例だけでなく、低品質データの例も示し、さらなる分析のための特定の制限について説明します。以前に発表されたいくつかのレポートでは、RNA-seq解析ツールの目的が説明されており、その性能の比較評価が提供されています^7,8。しかし、RNA-seqの品質管理と方法論の報告基準は標準化されておらず、トランスクリプトミクス実験の再現性と生物学的に意味のある解釈を悪化させます。

Rupは、標準的な高品質ツール、品質管理分析、結果の視覚化を統合しています。Rupは、生のシーケンシングリードと注釈付きゲノムを入力として必要とし、Mac OS、Linux、およびWindowsシステム上の「Windows Subsystem for Linux」(WSL)で実行されます。シーケンシングの品質のためのツールを統合し、ゲノムへのリードマッピングを定量化し、サンプル中のrRNA含有量測定を含め、反復相関を評価するためのサンプル相関を提供します。試験データは、植物種の2つの異なる段階の中央分裂組織のレーザー顕微解剖を使用して得られました Eschscholzia californica ライブラリ調製のための超低入力プロトコル、およびNovaseq 6000で配列決定されました。

プロトコル

Rup は、最小限の入力ファイルで 1 つのスクリプトとして実行できます。ここに示すように、パイプラインにはイルミナのペアエンドRNAシーケンシングデータが必要です。同じ分析ステップでシングルエンドシーケンス用に調整されたパイプラインも、GitHubリポジトリに保管されます。必要なのは、ゲノム配列をfastaファイルとして、遺伝子モデルとrRNAアノテーションをgtfファイルとして、生のシーケンシングリードを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」というファイルでソースフォルダに追加します。遺伝子モデルアノテーションを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 bp 未満にしないでください。読み取りが小さいとマッピングが中断される可能性があるためです。値を大きくすると、より一意にマップされた読み取りになりますが、シーケンスの品質に応じて、トリミング中に破棄される読み取りも増えます。これらのパラメータの詳細については、サプリメントに記載されています。

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 12}で利用可能な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¹⁴ パッケージを使用して、リードをリファレンスゲノムにマッピングします。まず、リファレンスゲノムファスタファイル(genome_fasta_file)のインデックスを構築します。
   注:このステップは、参照ゲノムごとに1回実行されます。
    

   # 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形式であることを確認してください。カウントは、ゲノムと1回一致する読み取りのみを行います。
    

   # 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. 同じサンプルまたは組織からの複製が、異なるサンプルまたは組織からの複製よりも相互に高い相関を示すことを確認します。反復相関のクラスター化ヒートマップをプロットし、生の読み取りカウントまたは Transcripts Per Million (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:遺伝子リード数の分類。 サンプル内のすべての遺伝子は、割り当てられたリードの数に基づいて5つのカテゴリのいずれかに分類されます。赤で示されているのは、リードが割り当てられていない遺伝子の数です。割り当てられたリードの合計数が少ないサンプル(s2_r1からs1_r4)は、10〜100の割り当てられたリードを持つ遺伝子の数が多く、1000を超える遺伝子の数が少ないです。発現率の低い遺伝子は、これらのサンプルで欠落する可能性があります。 この図の拡大版を表示するには、ここをクリックしてください。

結果

パイプラインは完全に R スクリプトとして実装されており、Linux および Mac OS オペレーティング システムでテストされています。Windows ユーザーは、Windows Subsystem for Linux (WSL) を使用できます。コードとテストデータは、GitHub リポジトリ (https://github.com/oliverrupp/rup) として利用できます。シーケンスデータは、ENA EBIプロジェクトPRJEB96400で入手できます。

Rupを分析に使用してバルクRNA-seqの品質管理中に遭遇する可能性のあるさまざまな問題を例示するために、2つの実際のサンプルから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 は、2つのサンプル(S1とS2)の相関ヒートマップをそれぞれ5回の反復で示しています。プロットの上部と側面にある樹状図は、各クラスターに5つの反復がある2つのクラスターを示しています。左側のクラスターには、サンプル 1 の 4 つの反復とサンプル 2 の 1 つの反復 (s2_r5) が含まれ、右側のクラスターには、サンプル 2 の 4 つの反復とサンプル 1 の 1 つの反復 (s1_r5) が含まれています。この場合、サンプル名s1_r5とs2_r5を再度入れ替えると、各クラスターには 1 つのサンプルのすべての反復が含まれ、反復のラベル付けミスを示している可能性があります。反復がクラスター化されないその他の理由としては、サンプル/条件/組織間の区別の欠如、またはシーケンス品質上の理由のみによる反復のクラスタリングが考えられます。後者は、トリミングとマッピング後に読み取りカウントが非常に少ないすべてのレプリケート、または読み取りカウントが非常に多いすべてのレプリケートがクラスターを形成する場合に発生する可能性があります。

図5:サンプル相関ヒートマップ。 サンプル相関ヒートマップは、ログ₂ 変換された TPM 値に基づいており、それぞれ 5 つのサンプルからなる 2 つの異なるクラスターを示します。生物学的/技術的複製は、他の組織/治療からの複製よりも相互に高い相関を示すことが期待されます。左側のクラスターには、サンプル 1 の 4 つの反復とサンプル 2 の 1 つの反復 (s2_r5) が含まれ、右側のクラスターには、サンプル 2 の 4 つの反復とサンプル 1 の 1 つの反復 (s1_r5) が含まれています。ヒートマップでのクラスタリングを説明するために、単一の反復でサンプルのスワッピングまたはバッチ効果の可能性をチェックする必要があります。 この図の拡大版を表示するには、ここをクリックしてください。

表1:RNA-seq品質評価パイプラインの比較。 詳細な分析については、 補足ファイル 1 を参照してください。 この表をダウンロードするには、ここをクリックしてください。

補足ファイル1:パラメータと参照の選択。 分析は、パラメーターと参照の選択が全体的なQC結果に及ぼす影響を示しています。 このファイルをダウンロードするには、ここをクリックしてください。

ディスカッション

差次的遺伝子発現解析の質は、各サンプルで配列決定され、複製¹⁶のリード数とサンプルあたりの複製数^17,18の2つの要因に大きく依存します。ここでは、各反復における遺伝子発現定量に適したリード数を決定するためのユーザーフレンドリーなRupパイプラインを紹介します。さまざまな指標により、研究者は反復で割り当てられたリード数が低い理由を理解し、サンプル相関の問題を特定することができます。Rupは植物RNA-seqサンプルの品質評価のために開発されましたが、他の真核生物にも同様に適しており、植物以外のサンプルについてはそれ以上の調整を必要としません(補足ファイル1)。

Rup の最初のステップでは、読み取りトリミングの前後の合計読み取りカウントを決定します。シーケンスされたリードの数によって検出可能な遺伝子の数が決まり、リード数が少なすぎると、多くの差次的に発現する遺伝子が検出されないままになります。品質のトリミングおよびフィルタリングプロセス中に廃棄されたリードの大部分は、RNAの分解を示している可能性があり、高度に分解された遺伝子の発現が過小評価される可能性があります。ただし、サンプルを下流の用途に使用できるかどうかは、対象生物と研究の目的によって異なります。たとえば、ほとんどの植物遺伝子の発現をキャプチャするには、私たちの経験では、3,000万から5,000万の読み取りが必要ですが、真菌の場合は1,000万で十分かもしれません。したがって、Rupは問題のあるサンプルを除外するための品質しきい値を定義せず、発生する可能性のあるさまざまな問題を特定するための指標を提供します。

2番目のステップ(マッピング品質)では、参照ゲノムに対するリードアライメントの精度と信頼性を評価し、下流の分析への適合性を指定します。すべてのシーケンスされたリードを遺伝子存在量の計算に使用できるわけではありません。ゲノムまたはトランスクリプトームにアラインメントしないリード、またはゲノム上の複数の場所にマッピングされたリードは、ほとんどの場合無視され¹⁹ 、全体的なリードカウントには寄与しません。パイプラインの 2 番目のステップの結果を使用して、存在量計算で読み取りが使用されない理由を理解できます。マップされていないリードの数が多い場合は、RNA抽出中の汚染(植物病原体や草食動物など)または不完全な参照ゲノムを示している可能性があります。参照ゲノムの遺伝子モデルが不完全な場合、多数の「特徴なし」リードが発生する可能性があります。マルチマップリードの大部分は、ライブラリ内の多数のrRNA遺伝子が原因であり、シーケンシングライブラリ調製プロセス中のrRNA除去が不十分であることを示しています(マルチマップリードが本当にrRNAに由来する場合、パイプラインにrRNAアノテーションファイルを提供することでさらに分析でき、パイプラインは可能なrRNAリードの数を自動的に計算します)。

同様に重要な3番目の品質指標は、同じサンプルの反復間の相関関係です。一般に、同じサンプルの繰り返し間の相関は、異なるサンプルの繰り返し間の相関よりも高くなる必要があります。反復間の相関が低い場合は、反復間の生物学的変動が大きい、サンプル/条件/組織間の類似性が高い、その他のバッチ効果、さらにはサンプルの交換や誤標識を示している可能性があります。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のすべてのモジュールの標準として出力の可視化が実装されているわけではありません。したがって、2つの代替ツールでは、読み取り数の少なさ、トリミングされたリードの割合が高い、反復外れ値の特定など、いくつかの品質管理の問題はテストされていません。Rup、RSeQC、RNA-SeQCの比較は、 補足ファイル1で入手できます。さらに、RupはRNA-SeQCまたはRSeQCを補完的に使用でき、たとえば、このパイプラインによって生成されたソートされたBAMファイルをこれらのツールの入力として使用できます。

現在、Rup は少数のサンプル用に最適化されていますが、品質トリミングや読み取りマッピングなど、最も時間のかかる手順は、コンピューター クラスターやクラウド インフラストラクチャなどで事前に計算できます。4 Gb を超えるゲノムの場合、Rsubread アライナーは 4 Gb 未満のゲノムに限定されるため、これは必須です。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(32 GB RAM)上で開発・テストされました
ゲットプット	R	1.20.4
ggplot2	R	3.5.2
MacBook Pro	アップル		パイプラインはmaxOS 15.4.1(16 Gb 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