はじめに

RNA-seqなどのデータを眺めていて、得られた発現変動遺伝子のGC含量やUTRの長さなどの情報を調べたくなることがあったのでその方法をまとめる。

動機についてもう少し詳しく

遺伝子の発現制御には様々なメカニズムが関わっていることが知られているが、単純に遺伝子の長さやGC含量の違いで制御が分かれているという報告がいくつかある。例えば、転写物の長さが老化に伴う遺伝子発現の変化と関連していることを示した論文や(参考1)、AT-richな転写物とGC-richなmRNAでは細胞内での主要な分解系や局在が変わることを示した論文もある(参考2)。また、3'UTRの長さの違いはmRNAの安定性や翻訳の効率に影響することも知られている(参考3)。

以上を踏まえて、自分の解析したデータセットでも発現変動遺伝子群と転写物長、GC含量、UTR長に何かしらの傾向が見られるのではないかと期待した。

使用するライブラリ

biomaRt
BioMartはEuropean Bioinformatics Institute(EBI)などが管理している遺伝子IDやGOなどのアノテーションを網羅しているデータベースで、異なるデータベース間で使用されているIDの変換やタンパク質のドメイン情報など遺伝子の機能に関連するアノテーション情報の取得の他、遺伝子の配列やタンパク質の配列情報の取得も行える。biomaRtはRからBioMartにアクセスするためのインターフェース。
dplyr
テーブルデータの処理用。

テスト用の遺伝子IDの準備

biomaRtで一度に大量の遺伝子のアノテーションを取ろうとすると、ここに書かれているように通信関連のエラーが出る。そのため、今回はテスト用にヒトのprotein_coding geneのensembl_gene_idを200個、biomaRtを使って抽出する。

library(biomaRt)

# useEnsemblを用いてEnsemblのデータベースにアクセスする。
# この時に種やEnsemblのreleaseのバージョンも指定できる。
# 今回はヒトのrelease 106を使用する。

ensembl <- useEnsembl(biomart = "ensembl",
                      dataset = 'hsapiens_gene_ensembl',
                      version = 106)

# getBMを用いてensembl_gene_idを取得する。
# attributesで取得したいattributeを指定する。
# filtersで指定したattributeが、valuesで指定した値と一致してるかどうかで取ってくる情報を選抜できる。

human_protein_coding_gene_list <- getBM(attributes = 'ensembl_gene_id',
                                        filters = 'biotype',
                                        values = 'protein_coding',
                                        mart = ensembl)

# 得られたヒトのprotein_coding geneのリストからランダムに200個抽出
set.seed(5)
gene_list <- sample(human_protein_coding_gene_list$ensembl_gene_id, 200)

テスト用のtranscript IDの取得とフィルタリング

欲しいのは遺伝子の転写物が持つUTRの配列長なので、前項で引っ張ってきた200遺伝子に対応するtranscriptのidが必要。しかしここで問題になるのが、一つの遺伝子がAlternative splicingにより複数のtranscriptを持つ場合である。Alternative splicingによって起こる変化に興味があるなら別だが、遺伝子ごとに代表的なtranscriptの情報だけ欲しい場合は、各遺伝子が持つ複数のisoformから代表的なtranscriptを選択しなくてはならない。

遺伝子ごとの代表的なtranscriptを選ぶ方法は色々あり、その選択に使えるアノテーションもbiomaRで取得することが可能である。以下が使えそうなアノテーションの一例(自分が見つけてないだけで他にも有用なアノテーションがあるかもしれません)。

ccds: コンセンサスCDS（CCDS）プロジェクトによって、ヒトとマウスのprotein coding transcriptのコアセットに割り振られているID。
transcript_tsl: Transcript Support Level (TSL)は、UCSCとEnsemblが提供するmRNAとESTのアラインメントに基づいて、転写物の構造がどれほどサポートできているかを示す指標。tsl1がベスト。
transcript_is_canonical: Ensembl Canonical transcriptは各遺伝子に割り振られた1つの代表的なtranscriptのこと。解析の際に1つの転写物のみが必要な場合に一貫性を持たせるために、遺伝子座ごとに単一の代表的な転写物をEnsemblが定めたもの。選抜アルゴリズムには配列長、保存度、発現量など様々な項目が用いられているらしい。

代表的なtranscriptのみでUTR長を計算したい場合はこれらのアノテーション情報をもとにtranscriptをフィルタリングすると良い。今回は'transcript_is_canonical == 1'の、Ensembl Canonical transcriptでアノテーションを取得する。

library(dplyr)

tmp0 <- getBM(attributes = c('ensembl_gene_id',
                             'ensembl_transcript_id',
                             'transcript_is_canonical'),
                filters = "ensembl_gene_id",
                values = gene_list,
                mart = ensembl) %>% dplyr::as_tibble()
  
transcript_list <- tmp0 %>%
    dplyr::filter(transcript_is_canonical == 1) %>%
    dplyr::pull(ensembl_transcript_id)

アノテーションの取得

フィルタリングして得られた200のtranscript_idを用いて、まずは転写物の長さの情報を取得する。

tmp1 <- getBM(attributes = c('ensembl_transcript_id',
                               'transcript_biotype',
                               'transcript_length'),
                filters = "ensembl_transcript_id",
                values = transcript_list,
                mart = ensembl) %>% dplyr::as_tibble()

続いて5'UTRのstartとendの座標を取得し、差を取ることでUTR長を計算する。

tmp2 <- getBM(attributes = c('ensembl_transcript_id',
                               '5_utr_start',
                               '5_utr_end'),
                filters = "ensembl_transcript_id",
                values = transcript_list,
                mart = ensembl) %>% 
    dplyr::as_tibble() %>%
    dplyr::mutate(`5_utr_length` = `5_utr_end` - `5_utr_start` + 1)

上記のコードで得られるテーブルを確認すると、一つのtranscriptに複数の5_utr_startと5_utr_endが紐ついている場合がある。これは、UTRが複数のexonにまたがっている場合に発生する(参照1, 参照2)。そのため、5'UTRの全長を計算するには5_utr_startと5_utr_endの差をtranscriptごとにグループ化して合計する必要がある。また、5_utr_startと5_utr_endがNAとなっている空白行も含まれているため、それらは除去する。

tmp2 <- getBM(attributes = c('ensembl_transcript_id',
                               '5_utr_start',
                               '5_utr_end'),
                filters = "ensembl_transcript_id",
                values = transcript_list,
                mart = ensembl) %>% 
    dplyr::as_tibble() %>%
    dplyr::mutate(`5_utr_length` = `5_utr_end` - `5_utr_start` + 1) %>%
    dplyr::filter(!is.na(`5_utr_length`)) %>%
    dplyr::group_by(ensembl_transcript_id) %>%
    dplyr::mutate(`5_utr_total_length` = sum(`5_utr_length`)) %>%
    dplyr::distinct(ensembl_transcript_id,`5_utr_total_length`)

同様に3'UTRの情報も取得する。

tmp3 <- getBM(attributes = c('ensembl_transcript_id',
                               'chromosome_name',
                               '3_utr_start',
                               '3_utr_end'),
                filters = "ensembl_transcript_id",
                values = transcript_list,
                mart = ensembl) %>%
    dplyr::as_tibble() %>%
    dplyr::mutate(`3_utr_length` = `3_utr_end` - `3_utr_start` + 1) %>%
    dplyr::filter(!is.na(`3_utr_length`)) %>%
    dplyr::group_by(ensembl_transcript_id) %>%
    dplyr::mutate(`3_utr_total_length` = sum(`3_utr_length`)) %>%
    dplyr::distinct(ensembl_transcript_id,`3_utr_total_length`)

最後に、cds長と遺伝子のGC含量の情報をそれぞれ取得する。一度に全てのアノテーションを取らないのは、アノテーションの組み合わせが複雑になり行数が増えるのを防ぐため。

tmp4 <- getBM(attributes = c('ensembl_transcript_id',
                               'cds_length'),
                filters = "ensembl_transcript_id",
                values = transcript_list,
                mart = ensembl) %>% dplyr::as_tibble()
  
tmp5 <- getBM(attributes = c('ensembl_gene_id',
                               'percentage_gene_gc_content'),
                filters = "ensembl_gene_id",
                values = gene_list,
                mart = ensembl) %>% dplyr::as_tibble()

最後に、得られたテーブルをgene_idとtranscript_idで統合する。

full_data <-  tmp0 %>%
    dplyr::filter(transcript_is_canonical == 1) %>%
    dplyr::left_join(tmp1, by = "ensembl_transcript_id")%>%
    dplyr::left_join(tmp2, by = "ensembl_transcript_id") %>%
    dplyr::left_join(tmp3, by = "ensembl_transcript_id") %>%
    dplyr::left_join(tmp4, by = "ensembl_transcript_id") %>%
    dplyr::left_join(tmp5, by = "ensembl_gene_id")

結果のバリデーション

最後に、集計したutrの長さとcds_lengthの合計がtranscript_lengthと一致するか、念の為確認する。

full_data %>%
    dplyr::mutate(total_length = `5_utr_total_length` + `3_utr_total_length`+ cds_length,
                  test = transcript_length == total_length) %>%
    dplyr::count(test)

# A tibble: 2 × 2
#  test      n
#  <lgl> <int>
# 1 TRUE    190
# 2 NA       10

190個のtranscriptで計算は合っているが、10個のtranscriptでNAとなっている。

full_data %>%
    dplyr::mutate(total_length = `5_utr_total_length` + `3_utr_total_length`+ cds_length,
                  test = transcript_length == total_length) %>%
    dplyr::filter(is.na(test))

#  A tibble: 10 × 13
#   ensembl_gene_id ensembl_transc…¹ ccds  trans…² trans…³ trans…⁴ trans…⁵ 5_utr…⁶ 3_utr…⁷ cds_l…⁸ perce…⁹ total…˟ test 
#   <chr>           <chr>            <chr> <chr>     <int> <chr>     <int>   <dbl>   <dbl>   <int>   <dbl>   <dbl> <lgl>
# 1 ENSG00000203661 ENST00000641363  "CCD… ""            1 protei…    2525      NA    1577     948    42.3      NA NA   
# 2 ENSG00000224320 ENST00000552498  ""    "tslNA"       1 protei…    5422      NA    4237    1185    46.4      NA NA   
# 3 ENSG00000233046 ENST00000446761  ""    "tslNA"       1 protei…    1266     300      NA     966    36.7      NA NA   
# 4 ENSG00000257062 ENST00000593147  ""    "tsl5"        1 nonsen…     519      NA     228     291    35.2      NA NA   
# 5 ENSG00000274699 ENST00000632763  ""    "tslNA"       1 protei…     960      NA      NA     960    40.3      NA NA   
# 6 ENSG00000275754 ENST00000616594  ""    "tslNA"       1 protei…     975      NA      NA     975    52.9      NA NA   
# 7 ENSG00000282307 ENST00000632258  ""    "tsl3"        1 protei…     254     125      NA     129    54.8      NA NA   
# 8 ENSG00000282318 ENST00000633447  ""    "tsl1"        1 protei…    8796      NA    7070    1726    47.2      NA NA   
# 9 ENSG00000282956 ENST00000634409  ""    ""            1 protei…    4389      NA      NA    4389    40.3      NA NA   
# 10 ENSG00000283961 ENST00000640795  ""    ""            1 protei…    1656      NA     405    1251    48.0      NA NA   
# … with abbreviated variable names ¹ensembl_transcript_id, ²transcript_tsl, ³transcript_is_canonical,
#   ⁴transcript_biotype, ⁵transcript_length, ⁶`5_utr_total_length`, ⁷`3_utr_total_length`, ⁸cds_length,
#   ⁹percentage_gene_gc_content, ˟total_length

上記のコードでは、集計する行のどれか一つにNAが入っている場合は合計もNAとなってしまう。そのため、5'UTR長と3'UTR長がNAになっているtranscriptが10個いることが、NAが10個いた原因だった。

5'UTR長がNAになっているENST00000641363の配列情報をensemblで確認すると、この遺伝子は5'UTRの情報が欠けていることがわかる。5'UTR長と3'UTR長の両方がNAになっているENST00000632763では5'UTRと3'UTR両方の情報が欠けている。UTRの開始と終点のアノテーションがまだついてないものでUTR長が計算できないのは仕方ないので、この挙動は合っている。

改めて、集計したutrの長さとcds_lengthの合計がtranscript_lengthと一致するか、NAがレコードに含まれる場合も考慮して確認する。

full_data %>%
    dplyr::rowwise() %>% 
    dplyr::mutate(total_length = sum(`5_utr_total_length`, `3_utr_total_length`, cds_length, na.rm=TRUE),
                  test = transcript_length == total_length) %>%
    dplyr::count(test)

# A tibble: 1 × 2
# Rowwise: 
#  test      n
#  <lgl> <int>
# 1 TRUE    200

これで今回アノテーションを取得した200のtranscript全てでtranscript_lengthと合計値があっていることが確認できた。

最後に

今回取得したGC含量(percentage_gene_gc_content)はtranscriptごとではなく遺伝子ごとに計算されているらしい(参照3)。TranscriptごとのGC含量を全長、5utr、cds、3utrそれぞれで計算した結果もマージできると良いかも。

また、今回取得したアノテーション(遺伝子長、GC含量、UTR長)に加えて、G/C ending codonの割合が発現変動と相関するパターンもあるらしい (参考4)。ゆくゆくは遺伝子のIDをもとにGC/ΑΤ ending codonの割合なども一度に取得できるようにしたい。

2023年8月1日更新
TranscriptごとのGC含量を全長、5utr、cds、3utrそれぞれで計算する方法と、G/C ending codonの割合を計算する方法をまとめた記事を7月28日にアップしました。