visium解析引継ぎ資料

Visium解析引継ぎ資料

解析フローとしては大きく以下の流れで行った。

  1. Space Ranger(というよりかはSTAR)と互換性を保つようにGTFファイルを編集
  2. Space Rangerのmkrefコマンドでインデックスを作成
  3. Space Rangerのcountコマンドでリファレンスへのマッピングとカウント、および組織の検出
  4. Web Summaryによる概要の確認
  5. Loupe Browserによる解析
  6. Seuratによる解析
  7. その他下流解析

10X Genomicsによって提供されているヒトやマウス以外の非モデル生物のカスタムリファレンスを使用したい場合は、GTFファイルを編集し、STARとの互換性を保つようにする必要がある。

互換性を保つためには以下の2つの処理が必要(もしかしたら他にも必要かも)

  1. GTFファイルの9列目のAttribute列にgene_idtranscript_idが入っている必要がある。
  2. コメント行を削除

Before

# Predicted genes for sequence number 1 on both strands
# start gene g00001
scaffold10|size694159   AUGUSTUS        gene    19      26628   0.11    +       .       ID=pfu_aug170726_10_00001
scaffold10|size694159   AUGUSTUS        transcript      19      26628   0.11    +       .       ID=pfu_aug170726_10_00001.t1;Parent=pfu_aug170726_10_00001
scaffold10|size694159   AUGUSTUS        transcription_start_site        19      19      .       +       .       Parent=pfu_aug170726_10_00001.t1
scaffold10|size694159   AUGUSTUS        exon    19      284     .       +       .       Parent=pfu_aug170726_10_00001.t1
scaffold10|size694159   AUGUSTUS        start_codon     224     226     .       +       0       Parent=pfu_aug170726_10_00001.t1
scaffold10|size694159   AUGUSTUS        intron  285     6779    1       +       .       Parent=pfu_aug170726_10_00001.t1
scaffold10|size694159   AUGUSTUS        intron  6929    8036    1       +       .       Parent=pfu_aug170726_10_00001.t1
scaffold10|size694159   AUGUSTUS        intron  8154    8308    1       +       .       Parent=pfu_aug170726_10_00001.t1

After

scaffold10|size694159   AUGUSTUS        gene    19      26628   0.11    +       .       ID=pfu_aug170726_10_00001; transcript_id "pfu_aug170726_10_00001"; gene_id "pfu_aug170726_10_00001";
scaffold10|size694159   AUGUSTUS        transcript      19      26628   0.11    +       .       ID=pfu_aug170726_10_00001.t1;Parent=pfu_aug170726_10_00001; transcript_id "pfu_aug170726_10_00001.t1"; gene_id "pfu_aug170726_10_00001.t1";
scaffold10|size694159   AUGUSTUS        transcription_start_site        19      19      .       +       .       Parent=pfu_aug170726_10_00001.t1; transcript_id "pfu_aug170726_10_00001.t1"; gene_id "pfu_aug170726_10_00001.t1";
scaffold10|size694159   AUGUSTUS        exon    19      284     .       +       .       Parent=pfu_aug170726_10_00001.t1; transcript_id "pfu_aug170726_10_00001.t1"; gene_id "pfu_aug170726_10_00001.t1";
scaffold10|size694159   AUGUSTUS        start_codon     224     226     .       +       0       Parent=pfu_aug170726_10_00001.t1; transcript_id "pfu_aug170726_10_00001.t1"; gene_id "pfu_aug170726_10_00001.t1";
scaffold10|size694159   AUGUSTUS        intron  285     6779    1       +       .       Parent=pfu_aug170726_10_00001.t1; transcript_id "pfu_aug170726_10_00001.t1"; gene_id "pfu_aug170726_10_00001.t1";
scaffold10|size694159   AUGUSTUS        intron  6929    8036    1       +       .       Parent=pfu_aug170726_10_00001.t1; transcript_id "pfu_aug170726_10_00001.t1"; gene_id "pfu_aug170726_10_00001.t1";
scaffold10|size694159   AUGUSTUS        intron  8154    8308    1       +       .       Parent=pfu_aug170726_10_00001.t1; transcript_id "pfu_aug170726_10_00001.t1"; gene_id "pfu_aug170726_10_00001.t1";

Space Rangerのmkrefコマンドを使用して、リファレンスのインデックスを作成する。

spaceranger mkref --genome=(インデックス名) --fasta=(リファレンスの
のFASTAファイル) --genes=(リファレンスのGTFファイル)

実際に使用したコマンドは以下の通り

spaceranger mkref --genome=akoya_ref_20211015 --fasta=raw/pfu_170726_scaffold.gapclosed.fasta --genes=raw/pfu_aug170726.modID_modified_20201020_filtered.gtf

ちなみにSpace Rangerのダウンロードは以下のリンク先から可能である。

Space Ranger Software Download

また最新版でない可能性があるが、伊藤がインストールしたSpace Rangerは/suikou/files/m208/ito.takumi/work/visium/space_ranger/spaceranger-1.0.0/bin/spaceranger にあるのでそれを使用してもよい。

Space Rangerのcountコマンドを用いて、組織検出、リードのマッピングとカウントを行った。

spaceranger count --id=(作成するディレクトリの名前) --transcriptome=(リファレンスインデックスのPATH) --fastqs=(シーケンスデータが格納されたディレクトリのPATH) --sample=(サンプル名、FASTQファイルの接頭辞?おそらくなくてもOK) --image=(組織切片の画像) --slide=(Visium SlideのID、Slideに記載されている) --area=(Slide上の切片を置いた位置) --localcores=(コア数) --localmem=(メモリ)

実際に使用したコマンドは以下のとおり

spaceranger count --id=mantle_20211015 --transcriptome=/home/ito.takumi/work/pinctada/visium_spaceranger3/ref/akoya_ref_20211015/ --fastqs=/home/ito.takumi/work/pinctada/visium_spaceranger3/mantle/input/ --sample=akoya_mantle --image=/home/ito.takumi/work/pinctada/visium_spaceranger3/image/20210718155118_pinctada_mantle.tif --slide=V10T17-119 --area=A1 --localcores=8 --localmem=64

countコマンドの出力としてoutsディレクトリが作成され、以下のファイルおよびディレクトリがoutsディレクトリ内に作成される。

analysis → 下流解析の結果
cloupe.cloupe → Loupe Browserの入力ファイルとして使用
filtered_feature_bc_matrix →
filtered_feature_bc_matrix.h5
metrics_summary.csv
molecule_info.h5
possorted_genome_bam.bam
possorted_genome_bam.bam.bai
raw_feature_bc_matrix
raw_feature_bc_matrix.h5
spatial
web_summary.html → Visium実験のサマリーを表示するHTMLファイル

Space Rangerのcountコマンドの出力であるWeb Summaryをもとに実験概要の確認を行う。

①Alertの確認

②シーケンス情報の確認(スポット当たりのリード数が50000リード程度あれば十分なシーケンス量を確保できているといえる)

③マッピング情報の確認

③スポットの情報の確認

④サンプル情報の確認

またこれらをグラフで可視化したものや、簡易的なクラスタリング図などを「Analysis」タブから確認することもできる。

Loupe Browserを用いることでGUIベースで空間遺伝子発現解析の解析を行うことができる。各遺伝子の空間遺伝子発現分布の取得、組織のクラスタリング、UMAPの描画などが主な機能となる。ただしクラスタリングに関しては細かいパラメータの調整などが不可能であるため、CUIベースのSeuratで行うことを推奨する。一方各遺伝子の空間遺伝子発現分布の取得においてはどちらでもそこまで変わりはないためLoupe Browserで行った方が楽に思われる(伊藤の感想としては、色調のパラメータ調整がLegendの表示がLoupe Browserの方が楽と感じた)。

入力データとしてはSpace Rangerの出力ファイルである、cloupe.cloupe を使用する。

①「Categories」を選択することでクラスタリングを行う。(Spatialで組織クラスタリング結果図を、UMAPでUMAP図を表示する)

②「Gene/Feature Expression」から各遺伝子の空間遺伝子発現分布を取得できる。

Seuratを用いて、各遺伝子の空間遺伝子発現分布や組織のクラスタリング、UMAPの描画、クラスターのマーカー遺伝子の取得、他の組織の空間遺伝子発現解析データとのマージなどを行うことができる。

Seurat解析手順

1.環境構築

①R 4.1.2をインストール (https://cran.r-project.org/index.html)

②R studioをインストール (https://www.rstudio.com/products/rstudio/download/#download)

③Seuratをインストール

> install.packages('Seurat')

④hdf5rをインストール

> install.packages("hdf5r")

2.解析

①必要なライブラリをロードする

> library(Seurat)
> library(ggplot2)
> library(patchwork)
> library(dplyr)

②Space Rangerの解析結果を読み込む

> mantle <- Load10X_Spatial("C:/kougaku/pinctada/visium2/Visium2/Result/mantle/outs_20211015", filename = "filtered_feature_bc_matrix.h5", assay = "mantle", slice = "C:/kougaku/pinctada/visium2/Visium2/Result/Image/Visium_20210718/20210718155118_pinctada_mantle.tif", filter.matrix = TRUE,to.upper = FALSE)
> hinge <- Load10X_Spatial("C:/kougaku/pinctada/visium2/Visium2/Result/hinge/outs_20211022", filename = "filtered_feature_bc_matrix.h5", assay = "hinge", slice = "C:/kougaku/pinctada/visium2/Visium2/Result/Image/20210902105928_pinctada_hinge.tif", filter.matrix = TRUE,to.upper = FALSE)

③バイオリンプロットの作成

> plot1 <- VlnPlot(mantle, features = "nCount_mantle", pt.size = 0.1) + NoLegend()
> plot2 <- VlnPlot(hinge, features = "nCount_hinge", pt.size = 0.1) + NoLegend()
> plot1 + plot2

④組織内のリードカウント数の分布を表示

> plot3 <- SpatialFeaturePlot(mantle, features = "nCount_mantle") + theme(legend.position = "right")
> plot4 <- SpatialFeaturePlot(hinge, features = "nCount_hinge") + theme(legend.position = "right")
> plot3 + plot4

⑤カウント数が0のスポットのデータを除去する。

> mantle_filtered <-  subset(mantle, subset =  nCount_mantle > 1)
> hinge_filtered <-  subset(hinge, subset =  nCount_hinge > 1)

⑥データの正規化を行う。

> mantle_scaled <- SCTransform(mantle_filtered, assay = "mantle", verbose = FALSE)
> hinge_scaled <- SCTransform(hinge_filtered, assay = "hinge", verbose = FALSE)

⑦mantle, hingeのデータをマージ

> akoya.merge <- merge(mantle_scaled, hinge_scaled)

⑧RNA発現データ上での次元削減とクラスタリング (1組織で行う場合は⑦をスキップし、⑥で得られたオブジェクトを引数として以下の解析を行えばよい)

> DefaultAssay(akoya.merge) <- "SCT"
> VariableFeatures(akoya.merge) <- c(VariableFeatures(mantle_scaled), VariableFeatures(hinge_scaled))
> akoya.merge <- RunPCA(akoya.merge, verbose = FALSE)
> akoya.merge <- FindNeighbors(akoya.merge, dims = 1:15)
> akoya.merge <- FindClusters(akoya.merge, verbose = FALSE)
> akoya.merge <- RunUMAP(akoya.merge, dims = 1:30)

⑨mergeデータのUMAPを作成

> DimPlot(akoya.merge, reduction = "umap")

⑩merge後の組織のクラスタリング図を作成(2組織とも表示)

> SpatialDimPlot(akoya.merge)

⑪merge後の組織のクラスタリング図を作成(1組織ずつラベルと一緒に表示)

> SpatialDimPlot(akoya.merge,label = TRUE, combine = FALSE)

⑫各クラスターのポジティブマーカーを抽出(閾値は適宜設定)

> pbmc.markers <- FindAllMarkers(akoya.merge, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
> View(pbmc.markers)
> write.csv(pbmc.markers, "cluster_markers.csv")

⑬各遺伝子の空間遺伝子発現を表示(MSI60=pfu-aug170726-71-03468.t1の場合)(_が-に変換されていることに注意)、本日は行いませんでしたが、Loupe Browserで行っているような解析をSeuratでも行う事ができます。

> SpatialFeaturePlot(akoya.merge, features = "pfu-aug170726-71-03468.t1")

⑭各遺伝子が発現しているスポットをUMAP上で表示(以下は複数遺伝子を同時に表示したい場合)

> FeaturePlot(akoya.merge, features = c("pfu-aug170726-71-03468.t1", "pfu-aug170726-869-31402.t1", "pfu-aug170726-1189-31713.t1", "pfu-aug170726-1547-25368.t1", "pfu-aug170726-3143-12855.t1", "pfu-aug170726-102-06870.t1"), pt.size = 0.2, ncol = 3)

※2021/12/10追記:特定のクラスタのみを抽出し、そのデータで再度PCA解析を行う。 ⑮特定のクラスター(この場合はクラスター7,11)のデータを抽出する。

> hinge_selected <- subset(hinge_scaled, idents = c(7,11))

再度PCA解析などを行う

> hinge_selected <- RunPCA(hinge_selected, verbose = FALSE)
> hinge_selected <- FindNeighbors(hinge_selected, dims = 1:15)
> hinge_selected <- FindClusters(hinge_selected, verbose = FALSE)
> hinge_selected <- RunUMAP(hinge_selected, dims = 1:30)

グラフを描画(crop=TRUEにしないと組織の一部分しか表示されないので注意)

> SpatialDimPlot(hinge_selected,label = TRUE,crop = FALSE)

その他の解析として、各クラスターのマーカー遺伝子の機能解析を行った。

手順としては以下の通り

①daimond(blastxの高速版)でuniprotデータベースに対して相同性検索。

/suikou/tool-all/bin/diamond blastx -d /home/ito.takumi/work/pinctada/diamond_db/uniprot_sprot_2021_04 -q /home/ito.takumi/work/pinctada/pinctada_id/tmp/${geneNam}.fa

②各クラスターの相同性を得られた遺伝子群を用いて、Enrichrでエンリッチメント解析。(https://maayanlab.cloud/Enrichr/

ここでUniprotに相同性検索を行った際に得られるsp|C1D8N3|SERC_LARHHの出力からSERCのみを抽出してEnrichrの入力としてあげる必要があることに注意。

③diamondでnrデータベースに対して相同性検索。

④Alphafoldで構造予測を行い、Dali(http://ekhidna2.biocenter.helsinki.fi/dali/)で構造に基づいた相同性検索を行う。Alphafoldを使用する際はtranseqを用いて塩基配列をアミノ酸配列に変換したのちに、ストップコドン(*)を除去してからAlphafoldにかけてやる必要がある。具体的なAlphafoldの使用方法はこちらを参照。

  • visium解析引継ぎ資料.1645782250.txt.gz
  • 最終更新: 2022/02/25 09:44
  • by 133.11.144.12