バルクRNAセクシークシーデータに基づくトランスクリプトミクス解析

Summary

本プロトコルは、生データから機能的濃縮解析までのバルクRNA-seqプロセスの解析のための完全なパイプラインを確立します。

Abstract

非アルコール性脂肪肝(NAFL)は通常良性の病気と考えられています。しかし、非アルコール性脂肪性肝炎(NASH)に進行すると、患者は末期肝疾患を発症するリスクが大幅に高まります。多くの研究がNAFLからNASHへの移行の分子メカニズムを解明しようとしています。高スループットシーケンシング技術(バルクRNA-seqなど)は、トランスクリプトームを調べることで、分子の発現、シグナル伝達経路の活性化、疾患進行に関連するその他の要因を明らかにすることで、研究者により深い理解をもたらしました。研究者が病気治療の潜在的な標的を特定するために分析できる豊富なオープンソースデータがあります。しかし、トランスクリプトームの上流解析のための効率的かつ信頼性の高いプロセスが不足しているため、関連研究は制約を受けています。ここでは、再現性が高くユーザーフレンドリーな上流解析と関連する差異遺伝子解析パイプラインを提供し、プライベートまたは公開データの標準化された処理と詳細な解析を実現します。パイプラインは4つのステップに分かれています:(1) データの品質管理;(2) 遺伝子マッピング;(3) 差異遺伝子解析;および(4)関数解析。このプロセスは、疾患変換の分子メカニズムを明らかにし、Bulk RNA-seqデータの解析を通じて潜在的な薬物標的や治療アプローチのスクリーニングを支援することを目的としています。

Introduction

非アルコール性脂肪肝疾患(NAFLD)は、世界的に最も多く見られる慢性肝疾患であり、人口の4分の1以上に影響を及ぼしています。近年、その発生率は劇的に増加しています。^1,2,3。増加する病気の負担、特にそのより進行した形態である非アルコール性脂肪肝炎(NASH)は、世界的な健康上の大きな課題であり、重い経済的負担となっています^。NAFLDの第一段階は非アルコール性脂肪肝(NAFL)で、炎症や線維化を伴い、NASHへと進行することがあります。後者は肝硬変や肝細胞癌(^HCC)を含む末期肝疾患への進行リスクを大幅に高めます^5,6,7。HCCの発症率と死亡率は^{NASHの増加}と関連しており^、2030年までに肝移植の主要な適応となると予想^{されています。}しかし、NAFLDの臨床進行は非常に異質^{であり、関連}薬剤の開発を著しく妨げる^{ため、関与}する分子メカニズムを正確に探求することが特に重要です。

細胞組成情報のバルクRNA-seqベースの獲得は、さまざまな疾患の病因を著しく明らかにすることができます。近年、モデル生物とヒトを対象に多数の大量のRNA-seq研究が行われ、NASH進行^13、14、15の遺伝子発現の違いを明らかにし、新たな治療標的を特定しようとしています。Xiongらは、RNAセクシークレーション解析に基づき、肝臓の非実質細胞(NPC)が細胞外マトリックスの形成や細胞接着などのプロセスに関与し、NASH¹⁶の進行に寄与していることを発見しました。Liらは肝細胞内のウィルムズ腫瘍1結合タンパク質(WTAP)が異所性脂質の蓄積と炎症を調節し、NASH形成を促進することを示しました¹⁷。バルクRNA-seq解析はNASHのメカニズムを解明する強力なツールですが、その結果は上流データの質に非常に敏感です。上流の実験操作や解析プロセスの異質性はデータの信頼性を著しく損なう可能性があり、その結果、真の生物学的情報を隠し、その後の解析の精度を妨げます。したがって、標準化された上流分析手順のセットを確立することが重要です。

単一細胞RNAシーケンシング(scRNA-seq)と比較して、バルクRNA-seqは実験設計と実用的応用の両面でいくつかの明確な利点を提供します。scRNA-seqは単一細胞レベルでの細胞異質性の同定や細胞型特異的な転写特徴の精密解析を可能にしますが、高コスト、複雑なデータ処理要件、低存量転写産物検出の感度の制限と関連しています¹⁸。対照的に、バルクRNA-seqはより深いシーケンス深度、低コスト、高いサンプルスループットを提供し、集団レベルの差異遺伝子発現解析や分子メカニズムの探求に特に適^{しています。}したがって、標準化された分析ワークフローに導かれれば、バルクRNA-seqは複雑な疾患の分子基盤を調査する効率的でコスト効率的かつ堅牢なアプローチであり続けます。

このプロトコルは、高いRNA完全性(RIN ≥ 7.0)かつ十分な入力RNA(サンプルあたり≥500 ng)を持つヒト組織から得られたバルクRNA-seqデータセット向けに特別に設計されています。アライメントおよび数値化ステップの信頼性を確保するためには、少なくとも10コアCPU、32GBのRAM、最低200GBの空きディスク容量を備えたローカルワークステーションが推奨されます。これらの要件を基に、プロトコルは詳細な運用指示や標準化されたパラメータ設定を含む効率的かつ使いやすい分析ワークフローを提供し、大規模なトランスクリプトミックデータを分析する研究者のニーズに応えます。

Protocol

デモンストレーション目的で、Lan Baiらが生成した公開データセットPRJNA1023502を用いて、上流および下流解析の各ステップを示すために用いられ^ました。このデータセットはオープンアクセスのNCBI SRAデータベースから提供されるため、追加の権限や倫理的承認は必要ありません。必要なソフトウェアおよびRパッケージのバージョンをすべて確認するには 、材料表 を参照してください。公開されているデータセットPRJNA1023502、6つの非NASH、6つのNAFL、6つのNASH肝RNA-seqサンプルで構成されています。このプロトコルでは、SRAデータベースからのデータ取得、品質管理(fastp)、アラインメント(HISAT2)、定量化(featureCounts)、下流の差分発現および機能豊か解析を含むバルクRNA-seqワークフローのすべてのステップを実証しました。

1. SRAツールキットのインストール

1. SRAツールキットの公式ウェブサイトを訪れ、バージョン3.2.1をダウンロードしてください。

2. 公開データのダウンロード

1. SRA番号を取得してください。
   1. 記事補足²⁰、データ可用性セクション、またはNCBI SRAデータベース内のキーワード検索でアクセスできます。
   2. ターミナルでprefetchと入力してダウンロードしてください。

3. 遺伝子カウントマトリックスの生成

1. 以下に説明する経路と設定を定義します。
   REFERENCE=~/reference/human/GRCh38/GRCh38.primary_assembly.genome.fa
GTF=~/reference/human/GRCh38/gencode.v44.annotation.gtf
INDEX=~/reference/human/GRCh38/GRCh38_index
FASTQ_DIR=~/SRA_tutorial/fastq
OUT_FASTP=~/RNAseq/fastp
OUT_HISAT2=~/RNAseq/hisat2
OUT_COUNTS=~/RNAseq/counts
​mkdir -p $FASTQ_DIR $OUT_FASTP $OUT_HISAT2 $OUT_COUNTS
   1. 公式GENCODEデータベース(https://www.gencodegenes.org/human/)からヒトリファレンスゲノム(GRCh38プライマリーアセンブリ)および対応する遺伝子注釈ファイル(gencode.v44、包括的な遺伝子注釈セット)をダウンロードしてください。
      for f in SRR*; do [[ ! $f =~ \.sra$ ]] && mv "$f" "$f.sra"; done
​for f in SRR*; do [[ ! $f =~ \.sra$ ]] && mv "$f" "$f.sra"; done
2. SRAファイルの名前を変更してください。
   1. すべてのSRAファイルが.sraファイル拡張子を使用していることを確認し、下流ツールによる正確な認識と処理を容易にしてください。
   2. 以下のコマンドを実行します。
      for f in SRR*; do [[ ! $f =~ \.sra$ ]] && mv "$f" "$f.sra"; done
3. SRAファイルをFASTQ形式に変換するために以下のコマンドを実行します:
   for f in *.sra; do fasterq-dump "$f" --split-files -O $FASTQ_DIR - e 20; done
4. 参照ゲノムからHISAT2インデックスを作成します。
   hisat2-build $REFERENCE $INDEX
5. 視覚的な進捗監視でFASTQファイルを処理します。
   for fq in $FASTQ_DIR/*.fastq; do
​sample=$(basename "$fq" .fastq)
   1. ペアエンドシーケンスデータの場合、次のコマンドを実行します
      for fq1 in $FASTQ_DIR/*_1.fastq; do
sample=$(basename "$fq1" _1.fastq)
​fq2=$FASTQ_DIR/${sample}_2.fastq
   2. fastpを使って品質管理とフィルタリングを行い、その結果を$OUT_FASTPディレクトリに出力します。単一端シーケンスデータの場合、以下のコマンドを実行します:
      fastp \
-i "${fq}" \
-o $OUT_FASTP/${sample}.clean.fastq \
-h $OUT_FASTP/${sample}.html \
-j $OUT_FASTP/${sample}.json \
-w 20
      ペアエンドシーケンスデータの場合、以下のコマンドを実行します。
      fastp \
-i "${fq}" \ -I "$fq2" \
-o $OUT_FASTP/${sample}_1.clean.fastq \
-O $OUT_FASTP/${sample}_2.clean.fastq \
-h $OUT_FASTP/${sample}.html \
-j $OUT_FASTP/${sample}.json \
​-w 20
   3. 各サンプルについては、HTMLレポートで視覚的な品質管理を行い、JSONレポートで読み取り品質、GC分布、リード冗長性、アダプター汚染などの構造化指標を確認してください。すべてのサンプルが以下の品質基準を満たすことを確認してください:Q30≥85%、アダプター汚染率<5%、異常なGCパターンなし。すべての出力ファイルを$OUT_FASTPディレクトリに保存します。
   4. HISAT2を使ってシーケンスアライメントを行い、SAMファイルを生成して$OUT_HISAT2ディレクトリに出力します。単一端シーケンスデータの場合、以下のコマンドを実行します:
      hisat2 -p 20 \ -x $INDEX \-U $OUT_FASTP/${sample}.clean.fastq \
-S $OUT_HISAT2/${sample}.sam
      ペアエンドシーケンスデータの場合、以下のコマンドを実行します。
      hisat2 -p 20 \-x $INDEX \-1 $OUT_FASTP/${sample}_1.clean.fastq \
-2 $OUT_FASTP/${sample}_2.clean.fastq \
​-S $OUT_HISAT2/${sample}.sam
   5. GRCh38プライマリーアセンブリのFASTAファイルから作成されたインデックスを用いて、品質フィルタ済みリードをヒト参照ゲノムにHISAT2アライメントします。シングルエンドおよびペアエンドのシーケンスデータの両方に対して例のコマンドを提供します。
   6. SAMをBAMに変換し、ソートし、インデックス化します。
      samtools view -@ 20 -bS $OUT_HISAT2/${sample}.sam \
| samtools sort -@ 20 -o $OUT_HISAT2/${sample}.sorted.bam
samtools index $OUT_HISAT2/${sample}.sorted.bam
​done
   7. 各サンプルごとにSAMファイルを生成し、SAMtoolsを使ってソートされインデックス化されたBAMファイルに変換します。高品質なヒトトランスクリプトミックサンプル(RNA完全性数、RIN ≥ 8.0)では、全体のアラインメント率が85%を超えていることを確認し、一般的なバルクRNA-seqサンプル(RIN ≥ 7.0)では、アラインメント率≥70%許容範囲を考慮します。
   8. featureCountを使って遺伝子定量を行います。
      featureCounts -T 20 -p -s 0 \
-a $GTF \
-o $OUT_COUNTS /${sample}.counts.txt \
$OUT_HISAT2/${sample}.sorted.bam
Done
   9. 各サンプルごとにfeatureCountsによって生成されるタブ区切りの出力ファイル(*.counts.txt)とその要約レポート(*.counts.txt.summary)を確認してください。ヒトバルクRNAセクシークシングのリード割り当て率が典型的な閾値≥70%を満たすこと;著しく低い率は、糸の向きの不一致、注釈の問題、またはアラインメント品質の低さを示すことがあります。この非鎖特異的RNA-seqデータセットにはパラメータ-s 0を使用します。ストランド固有のライブラリでは、コマンドで-s 0を-s 1または-s 2に置き換えます。
6. 以下に説明する遺伝子カウントマトリックスを作成します。
   1. 最初のサンプルの遺伝子IDとカウントでカウントマトリックスを初期化します。
      cut -f1 $(ls $OUT_COUNTS/*.counts.txt | head -1) > all_counts.txt
   2. 各サンプルからカウントを繰り返し加算し、単一の行列にします。
      for f in $OUT_COUNTS/*.counts.txt; do
cut -f7 "$f" | paste all_counts.txt - > tmp && mv tmp
all_counts.txt
​done
   3. カウント行列にサンプルID付きのヘッダー行を追加します。
      samples=$(ls *.counts.txt | sed 's/.counts.txt//' | paste -sd "\t")
echo -e "Geneid\t$samples" | cat - all_counts.txt > counts_matrix.txt
   4. GTFファイルから遺伝子長(遺伝子ごとのエクソン長の合計)を抽出します。
      awk '$3=="exon"{match($0,/gene_id "([^"]+)"/,a); if(a[1]!=""){len=$5-$4+1; gene_len[a[1]]+=len}} END{print "GENE_ID\tLENGTH"; for(g in gene_len) print g"\t"gene_len[g]}' \$GTF > gene_length.txt

4. 生計数マトリックス処理と遺伝子注釈

1. R環境の初期化および該当Rパッケージの読み込み。
   1. install.packages() 関数を使って tidyverse、ggplot2、ggrepel、RColorBrewer、ggridges、FactoMineR をインストールしてください。tidyverseパッケージはデータ操作やプロット作成に用いられています。GGPLOT2パッケージは可視化に使われます。Ggrepelパッケージは、重複しないテキストラベルに用いられます。カラーパレットにはRColorBrewerパッケージが使用されています。Ggridgesパッケージはリッジ区画に使われます。FactoMineRパッケージはPCAおよび多変量解析に使用されています。
   2. BiocManager::install() 関数を使ってbiomaRt、DESeq2、clusterProfilerをインストールしてください。biomaRtパッケージはEnsemblの遺伝子注釈に使用されます。DESeq2パッケージは微分式解析に用いられています。clusterProfilerパッケージは機能的エンリッチメント解析に使用されます。
   3. ステップ2で生成された元のカウントマトリックスを読み取り、biomaRtを使ってEnsembl IDをHGNC遺伝子名にマッピングし、重複したGeneidカラムを削除し、クリーンになったカウントマトリックスを作成します。元のカウント行列(counts_matrix.csv)を入力として使い、Ensemblの遺伝子IDを行、サンプルを列にします。クリーンなカウントマトリックスを出力(clean_counts_SRA.csv)として保存し、HGNC遺伝子記号を行、サンプルを列にします。
      mart <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
id_map <- getBM(attributes = c("ensembl_gene_id", "hgnc_symbol"),
filters = "ensembl_gene_id",
values = exprSet$GeneID,
mart = mart)
exprSet <- exprSet %>%
left_join(id_map, by = c("GeneID" = "ensembl_gene_id")) %>%
filter(!is.na(hgnc_symbol), hgnc_symbol != "") %>%
distinct(hgnc_symbol, .keep_all = TRUE) %>%
column_to_rownames("hgnc_symbol")

5. 遺伝子発現定量化

注:詳細なスクリプトについては 補足ファイル1 を参照してください。

1. 以下のコマンドを実行し、1キロベースあたりの転写本断片数、100万マッピングリードあたりの断片数(FPKM)および100万あたりの転写本数(TPM)を計算します。
   counts <- read.csv("output/clean_counts_SRA.csv", header=TRUE, row.names=1)
gene_len <- read.delim("data/gene_length.txt", header=FALSE, col.names=c("gene_symbol","length"))
gene_len <- gene_len %>% distinct(gene_symbol, .keep_all=TRUE)
rownames(gene_len) <- gene_len$gene_symbol
gene_len <- gene_len[match(rownames(counts), gene_len$gene_symbol),]
length_bp <- gene_len$length
fpkm <- (counts / length_bp) * 1e9 / colSums(counts)
write.csv(fpkm, "output/clean_fpkm_SRA.csv")
tpm <- (counts / length_bp) / colSums(counts / length_bp) * 1e6
write.csv(tpm, "output/clean_tpm_SRA.csv")

6. サンプルクラスタリングと差分可視化

1. PCAによるサンプルクラスタリングを検証してデータの質を評価します。詳細なスクリプトについては 補足ファイル1 を参照してください。
2. PCAの可視化には、以下のコードをご利用ください。
   gene.pca <- PCA(exprSet, ncp = 2, scale.unit = TRUE, graph = FALSE)
ggplot(pca_sample, aes(x = Dim.1, y = Dim.2)) +
geom_point(aes(color = group)) +
labs(x = paste('PC1:', pca_eig1, '%'),
y = paste('PC2:', pca_eig2, '%'))

7. 微分表現解析および結果の可視化

注:詳細なスクリプトについては 補足ファイル1 を参照してください。

1. 以下のコマンドを実行してDESeq2データセットを構築し、閾値を設定します。すべてのサンプルで総数が1>の遺伝子のみがモデル化前に保持されました。調整済みp<値の有意閾値を用いて差分発現遺伝子(DEGs)を特定します。0.05および|log2FC|> 0.5。
   dds <- DESeq(DESeqDataSetFromMatrix(countData = exprSet, colData = colData, design = ~group)); sizeFactors(dds); res <- results(dds); dds <- dds[rowSums(counts(dds)) > 1,]
dd1 <- results(dds, contrast = contrast, alpha = 0.05)
dd2 <- lfcShrink(dds, contrast = contrast, res = dd1, type = "ashr")
2. 火山プロットで差異発現遺伝子を可視化します。
   ggplot(data = data, aes(x = log2FoldChange, y = -log10(padj))) +
geom_point(aes(color = group), alpha = 1, size = 1.2) +
geom_hline(yintercept = -log10(0.05), lty = 4) +
geom_vline(xintercept = c(-0.5, 0.5), lty = 4) +
geom_text_repel(data = subset(data, abs(log2FoldChange) >= 1.5 & padj < 0.05),
aes(label = gene_id))

8. 機能的濃縮分析および可視化を行う

注:詳細なスクリプトについては 補足ファイル1 を参照してください。

1. KEGG濃縮結果を可視化します。
   EGG <- enrichKEGG(gene = gene$ENTREZID, organism = 'hsa',
pvalueCutoff = 0.05, qvalueCutoff = 0.05)
ggplot(symboldata, aes(richFactor, Description)) +
geom_point(aes(color = p.adjust, size = Count))
2. GOの濃縮効果を視覚化してください。実験仕様に応じてアーク解像度(n=500)、色マッピング、ファセット配置を修正して可視化をカスタマイズできます。
   ego <- enrichGO(gene = gene$ENTREZID, OrgDb = "org.Hs.eg.db", ont = "ALL",
pvalueCutoff = 0.05, qvalueCutoff = 0.05, pAdjustMethod = "BH")
ggplot(df) +
ggforce::geom_link(aes(x = 0, y = Description, xend = -log10(p.adjust),
yend = Description, color = ONTOLOGY), n = 500, show.legend = FALSE) +
facet_wrap(~ONTOLOGY, scales = "free", ncol = 1)
3. 遺伝子セット豊か解析(GSEA)の結果を可視化します。異なる効果サイズ分布に対応するために、経路数、密度リッジのスケーリング、層の透明度を調整してください。
   genelist <- sort(res$log2FoldChange, decreasing = TRUE)
names(genelist) <- rownames(res)
hallmarks <- read.gmt('resource/h.all.v2023.2.Hs.symbols.gmt')
y <- GSEA(genelist, TERM2GENE = hallmarks, pvalueCutoff = 0.05)
gsearesult <- yd %>% arrange(desc(NES)) %>% slice_head(n = 10)
ggplot(gsearesult, aes(x = logFC, y = Description, fill = -log10(pvalue))) +
geom_density_ridges(alpha = 0.8, scale = 0.8) +
geom_point(aes(size = abs(NES), x = -0.4, color = NES)) +
scale_fill_distiller(palette = 'Spectral') +
scale_color_distiller(palette = 'Reds') +
scale_size_continuous(range = c(2, 6))

Results

バルクRNA-seqの上流解析ワークフローは 図1Aに示されています。このワークフローはLinuxプラットフォーム上で以下の重要なステップを順次実行します。まず、fastpを使って生のシーケンスデータの厳格な品質管理を行い、低品質のリードやアダプターシーケンスを除去します。その後、HISAT2は高品質なリードをリファレンスゲノムにアラインメントし、Samtoolsがアラインメントファイルを変換・ソートします。最後に、FeatureCountsは遺伝子レベルの定量化を行い、遺伝子発現マトリックスを生成するため、後続解析のための高品質な入力を提供します。その後の処理と統計解析はR環境内で行われ、関連するワークフローと必要なソフトウェアパッケージは 図1Bに示されています。分析のために発表された研究から選ばれ、非NASH6件、NAFL6件、NASH^{サンプル6} 件(図1C)から構成されています(図1C)。非NASH対照サ...

Discussion

バルクRNA-seqデータ解析は、ゲノミクス、バイオインフォマティクス、統計学、コンピュータサイエンスを統合した学際的な課題として特徴づけられます。完全な分析ワークフローは、生データの前処理、品質管理、配列アラインメント、遺伝子レベルの定量化、データ正規化、差異発現解析、生物学的解釈など、複数の上流および下流のステップを含みます。これらのステップの中で、生のシーケンシングリードを高品質な遺伝子発現マトリックスに正確に変換することは特に重要であり、上流処理中に導入された誤りが下流の生物学的結論に影響しやすいためです。したがって、トランスクリプトミック研究の再現性向上には、透明かつ標準化された上流解析ワークフローの確立が不可欠です。

このプロトコルは、fastp(読み取りトリミングや品質管理用)、HISAT2(スプライス認識アラインメント用)、featureCounts(遺伝子レベルの定量化用)など広く使われているツールを統合した、効率的で完全なスクリプトベースのワークフローを提供します。これらのツールは、Tuxedo...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
バイオマルト	バイオコンダクター	2.64.0	Ensemblによる遺伝子注釈
clusterProfiler	バイオコンダクター	4.16.0	機能豊化解析
DESeq2	バイオコンダクター	1.48.1	微分表現解析
ファクトマインR	アグロパリテック	2.11.0	PCAおよび多変量解析
ファストプ	オープンジーン	1.0.1	FASTQデータの品質管理とフィルタリング
フィーチャカウント	ウォルター&エリザ・ホール医療研究所バイオインフォマティクス部門	2.0.0	遺伝子発現定量のために各遺伝子にマッピングされたリード数を数えます
ggplot2	ポジット	3.5.2	データ可視化
グレッペル	カミル・スワヴィコフスキ	0.9.6	重複しないテキストラベル
グリッジズ	クラウス・O・ウィルケ	0.5.6	尾根線の区画を作成
HISAT2	ジョンズ・ホプキンス大学	2.2.1	フィルタリング済みの高品質リードを参照ゲノムにアラインメントします
R	Rコアチーム	4.5.0	データ計算、分析、可視化のための環境
RColorBrewer	エーリッヒ・ノイヴィルト	1.1.3	プロット用のカラーパレット
サムツールズ	大規模ゲノミクスの作業分野	1.22.0	効率的な取得とアクセスのためにSAMファイルの変換と処理
SRAツールキット	国立バイオテクノロジー情報センター	3.2.1	NCBI SRAデータベースから生のシーケンスデータを取得・前処理