はじめに

前回書いた記事で遺伝子ごとにUTR長やGC含量を取得する方法を書いた。その時に利用したデータベースの情報では、GC含量(percentage_gene_gc_content)はtranscriptごとではなく遺伝子ごとに計算されていたため、transcriptごとに再計算する方法を知りたいと考えていた。

また、前回取得したアノテーション(遺伝子長、GC含量、UTR長)に加えて、G/C ending codonの割合が発現変動と相関するパターンもあること(参考4)から、G/C ending codonの割合も計算したいと考えていた。

そこで今回は、Transcriptごとに各領域のGC含量とG/C ending codonの割合を算出する方法を記事にまとめる。

使用するライブラリ

前回と同様

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

テスト用の遺伝子IDの準備とtranscript IDの取得とフィルタリング

前回と同様、テスト用にヒトのprotein_coding geneのensembl_gene_idを200個、biomaRtを使って抽出する。そのあと、各遺伝子が持つEnsembl Canonical transcriptのIDを取得する。

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)


library(dplyr)
tmp0 <- getBM(attributes = c('hgnc_symbol',
                             '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)

アノテーションの取得

TranscriptごとにGC含量を計算するには、各transcriptの配列情報が必要になる。今回はbiomaRtのgetSequence関数を用いて配列情報を取得し、GC含量を計算することにする。

tmp1_3utr <- getSequence(type = 'ensembl_transcript_id',
                    seqType = '3utr',
                    id = transcript_list,
                    mart = ensembl) %>%
  dplyr::as_tibble() %>%
  dplyr::mutate(
    length = case_when(
      `3utr` == "Sequence unavailable" ~ 0,
      `3utr` != "Sequence unavailable" ~ as.numeric(stringr::str_length(`3utr`))),
    gc_content = (stringr::str_count(`3utr`, "G") +
                    stringr::str_count(`3utr`, "C")) / length * 100
    )


tmp1_5utr <- getSequence(type = 'ensembl_transcript_id',
                    seqType = '5utr',
                    id = transcript_list,
                    mart = ensembl) %>%
  dplyr::as_tibble() %>%
  dplyr::mutate(
    length = case_when(
    `5utr` == "Sequence unavailable" ~ 0,
    `5utr` != "Sequence unavailable" ~ as.numeric(stringr::str_length(`5utr`))),
    gc_content = (stringr::str_count(`5utr`, "G") +
                    stringr::str_count(`5utr`, "C")) / length * 100
  )


tmp1_cds <- getSequence(type = 'ensembl_transcript_id',
                    seqType = 'coding',
                    id = transcript_list,
                    mart = ensembl) %>%
  dplyr::as_tibble() %>%
  dplyr::mutate(
    length = case_when(
      coding == "Sequence unavailable" ~ 0,
      coding != "Sequence unavailable" ~ as.numeric(stringr::str_length(coding))),
    gc_content = (stringr::str_count(coding, "G") +
                    stringr::str_count(coding, "C")) / length * 100,
  )


tmp1_cdna <- getSequence(type = 'ensembl_transcript_id',
                    seqType = 'cdna',
                    id = transcript_list,
                    mart = ensembl) %>%
  dplyr::as_tibble() %>%
  dplyr::mutate(
    length = case_when(
      cdna == "Sequence unavailable" ~ 0,
      cdna != "Sequence unavailable" ~ as.numeric(stringr::str_length(cdna))),
    gc_content = (stringr::str_count(cdna, "G") +
                    stringr::str_count(cdna, "C")) / length * 100
  )

これでtranscriptの各領域でのGC含量が計算できた。次はCDSからG/C ending codonの割合を算出する。G/C ending codonの割合の計算するため、まずは各codonの3番目の塩基を取得し配列を作る。その後、3番目の塩基をまとめた配列を用いて、GC塩基の割合を計算する。

check_3rd_letter_of_codon <- function(coding) {
  
  vcoding <- strsplit(coding, "")[[1]]
  bool_vec <- rep(c(FALSE,FALSE,TRUE), stringr::str_length(coding)/3)
  third_letters <- paste(vcoding[bool_vec], collapse = "")
  gc_ending_codon_rate = (
    stringr::str_count(third_letters, "G") + 
    stringr::str_count(third_letters, "C")
    ) / stringr::str_length(third_letters) * 100
  
  return(gc_ending_codon_rate)
}

vcheck_3rd_letter_of_codon <- Vectorize(check_3rd_letter_of_codon)


tmp1_cds <- tmp1_cds %>%
    dplyr::mutate(
    gc_ending_codon_rate = vcheck_3rd_letter_of_codon(coding)
)

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

full_data <-  tmp0 %>%
    dplyr::filter(transcript_is_canonical == 1) %>%
    dplyr::left_join(tmp1_5utr, by = "ensembl_transcript_id")%>%
    dplyr::left_join(tmp1_3utr, by = "ensembl_transcript_id",
                     suffix = c(".5utr", ".3utr")) %>%
    dplyr::left_join(tmp1_cds, by = "ensembl_transcript_id") %>%
    dplyr::left_join(tmp1_cdna, by = "ensembl_transcript_id",
                     suffix = c(".cds", ".cDNA"))

結果のバリデーション

full_dataをensembl_gene_idでソートした場合、先頭に来るのはPRSS21のレコード。　　
PRSS21のレコードを確認してみる。　　

full_data <- full_data %>%
    arrange(ensembl_gene_id)

full_data[1,] %>%
    dplyr::select(hgnc_symbol:transcript_is_canonical,
                  contains('length'),
                  contains('content'),
                  contains('rate')) %>% t()

#                         [,1]             
# hgnc_symbol             "PRSS21"         
# ensembl_gene_id         "ENSG00000007038"
# ensembl_transcript_id   "ENST00000005995"
# transcript_is_canonical "1"              
# length.5utr             "32"             
# length.3utr             "114"            
# length.cds              "945"            
# length.cDNA             "1091"           
# gc_content.5utr         "78.125"         
# gc_content.3utr         "51.75439"       
# gc_content.cds          "59.25926"       
# gc_content.cDNA         "59.02841"       
# gc_ending_codon_rate    "71.42857"

Ensemblのレコードを用いて結果のバリデーションを行う。　　
上のリンク先のページの"cDNA"タブにいき、そこの"Download sequence"から5UTR, 3UTR, CDS, cDNAの配列がFASTA形式で取得できる。

>ENST00000005995.8 PRSS21-201 cdna:protein_coding
AGAGGGGGCGTCAGGCCGCGGGAGAGGAGGCCATGGGCGCGCGCGGGGCGCTGCTGCTGG
CGCTGCTGCTGGCTCGGGCTGGACTCAGGAAGCCGGAGTCGCAGGAGGCGGCGCCGTTAT
CAGGACCATGCGGCCGACGGGTCATCACGTCGCGCATCGTGGGTGGAGAGGACGCCGAAC
TCGGGCGTTGGCCGTGGCAGGGGAGCCTGCGCCTGTGGGATTCCCACGTATGCGGAGTGA
GCCTGCTCAGCCACCGCTGGGCACTCACGGCGGCGCACTGCTTTGAAACCTATAGTGACC
TTAGTGATCCCTCCGGGTGGATGGTCCAGTTTGGCCAGCTGACTTCCATGCCATCCTTCT
GGAGCCTGCAGGCCTACTACACCCGTTACTTCGTATCGAATATCTATCTGAGCCCTCGCT
ACCTGGGGAATTCACCCTATGACATTGCCTTGGTGAAGCTGTCTGCACCTGTCACCTACA
CTAAACACATCCAGCCCATCTGTCTCCAGGCCTCCACATTTGAGTTTGAGAACCGGACAG
ACTGCTGGGTGACTGGCTGGGGGTACATCAAAGAGGATGAGGCACTGCCATCTCCCCACA
CCCTCCAGGAAGTTCAGGTCGCCATCATAAACAACTCTATGTGCAACCACCTCTTCCTCA
AGTACAGTTTCCGCAAGGACATCTTTGGAGACATGGTTTGTGCTGGCAATGCCCAAGGCG
GGAAGGATGCCTGCTTCGGTGACTCAGGTGGACCCTTGGCCTGTAACAAGAATGGACTGT
GGTATCAGATTGGAGTCGTGAGCTGGGGAGTGGGCTGTGGTCGGCCCAATCGGCCCGGTG
TCTACACCAATATCAGCCACCACTTTGAGTGGATCCAGAAGCTGATGGCCCAGAGTGGCA
TGTCCCAGCCAGACCCCTCCTGGCCGCTACTCTTTTTCCCTCTTCTCTGGGCTCTCCCAC
TCCTGGGGCCGGTCTGAGCCTACCTGAGCCCATGCAGCCTGGGGCCACTGCCAAGTCAGG
CCCTGGTTCTCTTCTGTCTTGTTTGGTAATAAACACATTCCAGTTGATGCCTTGCAGGGC
ATTCTTCAAAA
>PRSS21-201 cds:protein_coding
ATGGGCGCGCGCGGGGCGCTGCTGCTGGCGCTGCTGCTGGCTCGGGCTGGACTCAGGAAG
CCGGAGTCGCAGGAGGCGGCGCCGTTATCAGGACCATGCGGCCGACGGGTCATCACGTCG
CGCATCGTGGGTGGAGAGGACGCCGAACTCGGGCGTTGGCCGTGGCAGGGGAGCCTGCGC
CTGTGGGATTCCCACGTATGCGGAGTGAGCCTGCTCAGCCACCGCTGGGCACTCACGGCG
GCGCACTGCTTTGAAACCTATAGTGACCTTAGTGATCCCTCCGGGTGGATGGTCCAGTTT
GGCCAGCTGACTTCCATGCCATCCTTCTGGAGCCTGCAGGCCTACTACACCCGTTACTTC
GTATCGAATATCTATCTGAGCCCTCGCTACCTGGGGAATTCACCCTATGACATTGCCTTG
GTGAAGCTGTCTGCACCTGTCACCTACACTAAACACATCCAGCCCATCTGTCTCCAGGCC
TCCACATTTGAGTTTGAGAACCGGACAGACTGCTGGGTGACTGGCTGGGGGTACATCAAA
GAGGATGAGGCACTGCCATCTCCCCACACCCTCCAGGAAGTTCAGGTCGCCATCATAAAC
AACTCTATGTGCAACCACCTCTTCCTCAAGTACAGTTTCCGCAAGGACATCTTTGGAGAC
ATGGTTTGTGCTGGCAATGCCCAAGGCGGGAAGGATGCCTGCTTCGGTGACTCAGGTGGA
CCCTTGGCCTGTAACAAGAATGGACTGTGGTATCAGATTGGAGTCGTGAGCTGGGGAGTG
GGCTGTGGTCGGCCCAATCGGCCCGGTGTCTACACCAATATCAGCCACCACTTTGAGTGG
ATCCAGAAGCTGATGGCCCAGAGTGGCATGTCCCAGCCAGACCCCTCCTGGCCGCTACTC
TTTTTCCCTCTTCTCTGGGCTCTCCCACTCCTGGGGCCGGTCTGA
>PRSS21-201 utr3:protein_coding
GCCTACCTGAGCCCATGCAGCCTGGGGCCACTGCCAAGTCAGGCCCTGGTTCTCTTCTGT
CTTGTTTGGTAATAAACACATTCCAGTTGATGCCTTGCAGGGCATTCTTCAAAA
>PRSS21-201 utr5:protein_coding
AGAGGGGGCGTCAGGCCGCGGGAGAGGAGGCC

この配列を用いてGC含量・配列長・G/C ending codonの割合を計算した結果は自分が計算した結果と合っていた。

最後に

これでmRNAの配列の特徴はかなり集計できるようになった。あとはintronやexonの特徴を集められるようになりたい。自分が調べた限りでは、getSequence()関数ではtranscriptごとのexon/intron配列やその位置情報を取得することはできない。そのためexon/intronで似たような情報を取るには、exon/intronに関する座標情報をtransciptごとに調べて、その座標情報を使って配列を取ってくる必要がありそう。