KaryoploteRで染色体を描く

KaryoploteR

karyoploteRはRを用いて核型を描画できるパッケージである。
https://github.com/bernatgel/karyoploteR

GRCh38やGRCh37などの頻繁に利用されるreference配列の描画は下記のようなパッケージなどを用いて対応されており、ワンコマンドで済む。
https://www.bioconductor.org/packages/release/data/annotation/html/BSgenome.Hsapiens.UCSC.hg38.html

一方で、T2Tや非モデル生物などは未対応なので、自分でカスタムゲノムを作成する必要がある。今回はT2Tのゲノムを描画してみることにする。

GRCh38に関する図を描画する

kp <- plotKaryotype(genome="hg38")

kp <- plotKaryotype(genome="hg38", chromosomes = "chr1")

biomaRtを噛ませることで、遺伝子名も自由自在に表記できる。
参照: https://bernatgel.github.io/karyoploter_tutorial/

library(biomaRt)
library(regioneR)

ensembl <- useEnsembl(biomart="ensembl", dataset="hsapiens_gene_ensembl")
genes <- toGRanges(getBM(attributes=c('chromosome_name', 'start_position', 'end_position', 'hgnc_symbol'),
               filters = 'hgnc_symbol', values =gene.symbols, mart = ensembl))
seqlevelsStyle(genes) <- "UCSC"

head(genes)

GRanges object with 1 range and 1 metadata column:
    seqnames            ranges strand | hgnc_symbol
       <Rle>         <IRanges>  <Rle> | <character>
  1     chr1 50960745-50974634      * |      CDKN2C
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths

kp <- plotKaryotype(genome="hg38", chromosomes = "chr1")
kpPlotMarkers(kp, data=genes, labels=genes$hgnc_symbol)

今回は記載しないが、場所とそれに付随する情報があれば基本的にデータとして利用できる。例えば、bamファイルや統計値などをまとめたtsvファイルがあれば、染色体毎のリード深度やマンハッタンプロットなども作成可能である。

T2Tのkpを作成する

T2Tの別名hs1で動かしてみると下記のようなエラーが出る。

> kp <- plotKaryotype(genome="hs1")
'getOption("repos")' replaces Bioconductor standard repositories, see 'help("repositories", package =
"BiocManager")' for details.
Replacement repositories:
    CRAN: https://cran.rstudio.com/
'getOption("repos")' replaces Bioconductor standard repositories, see 'help("repositories", package =
"BiocManager")' for details.
Replacement repositories:
    CRAN: https://cran.rstudio.com/
'getOption("repos")' replaces Bioconductor standard repositories, see 'help("repositories", package =
"BiocManager")' for details.
Replacement repositories:
    CRAN: https://cran.rstudio.com/
Error in value[[3L]](cond) : 
  It was not possible to identify or load the requested genome. Error in .stopOnAvailablePkg(genome): BSgenome.Hsapiens.UCSC.hs1 package is not currently installed.
  You first need to install it, which you can do with:
      library(BiocManager)
      install("BSgenome.Hsapiens.UCSC.hs1")

エラーに従い、以下のように進めてみる。

library(BiocManager)
install("BSgenome.Hsapiens.UCSC.hs1")

kp <- plotKaryotype(genome="hs1", hromosomes="all")

しかし、線しか出て来ず、GRCh38のようなかっこいいkaryotypeにはならない。

Custom genomeとしてT2Tを描画する

そこでcustom genomeを描画してみる。
参照: https://bernatgel.github.io/karyoploter_tutorial//Tutorial/CustomGenomes/CustomGenomes.html

行うべきこととしては、mygenome.txtおよびmycytobands.txtの作成である。

mygenome.txtはchr, start, endからなるtxtファイルが必要である。
https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_009914755.1/
このページのDownloadからtsvファイルを開けるので、それをそのまま利用する。

custom.genome <- toGRanges("~/Downloads/mygenome_T2T.tsv")

mycytobands.txtは以下のサイトのCytobandsに落ちていたので、それをダウンロードして加工する。具体的にはchr start, end, name, gieStainの5行からなるtxtファイルが必要なので、最初に行を足せば良い。
https://github.com/marbl/CHM13?tab=readme-ov-file

custom.cytobands <- toGRanges("~/Downloads/chm13v2.0_cytobands_allchrs.txt")

最後に上で作成したファイルを利用して、plotKaryotypeする。

kp <- plotKaryotype(genome = custom.genome, chromosomes = "all", cytobands = custom.cytobands)

終わりに

P.S. 以下にファイルを置いておきました。
https://github.com/geedrn/T2T_KaryoploteR/tree/main