概要

この記事の目的では1回のRNA-Seqで、様々な条件で比較するようなサンプルが並んでいる場合における、なるべく条件を揃えた上でまとめて実行できないか、というのを試した例である。元データはsubread FeatureCountsで得られたカウントテーブルで、今回はedgeRにてそれぞれのボルケーノプロットと、それぞれの条件におけるテストのテーブルを出力する。¹

なお、なにかご意見ありましたらお気軽にコメントください。

条件

以下の様な条件を想定する。

GRCh38でgencodeのannotation (PAR_Yが入っている)
4グループあり、1-2, 1-3, 1-4 の様に独立に見ていきたい

本文

前準備1

HISAT2等でアライメントしたbamファイルをsortし、得られたそれぞれのファイルをカウントデータにする。

featureCounts -T 32 \
	-g gene_id \
	-a ~/db/GRCh38/annotation.gtf \
	-o counttab.txt \
	*.sort.bam

前準備2

gene_idでカウントしたためカウントテーブルではEnsemblのgene idが行名部分に入るが、これにGene Nameなどをつけたいため、EnsemblのBioMartからGene Stable ID, Gene Name, Gene Descriptionの順番でテーブルを作成しておく。(marttable.tsv) これは{biomaRt}でも可能に見えるが、こちらの環境では上手くいかなかったため泥臭い方法をとっている。

前準備3

メタデータテーブルを作る。ファイル名を仮定して、以下の様に作成する。(metadata.tsv)

sampleid	group
Ctrl_01	0
Ctrl_02	0
Ctrl_03	0
TR_A_01	A
TR_A_02	A
TR_A_03	A
TR_B_01	B
TR_B_02	B
TR_B_03	B
TR_C_01	C
TR_C_02	C
TR_C_03	C

なお、今回の例では、Ctrl対TR_A, TR_B, TR_Cをそれぞれ実行する例を示す。なお、この下に出てくるデータは全て例であり、実際のデータとは一切関係ない。

データの読込

ライブラリをそれぞれ読みこむ。

library(edgeR)
library(tidyverse) # dplyr, ggplot2, purrr, etc.
# library(datawizard)
# library(plotly) この2つに関しては後で使うがlibraryでは読み込まない。
# この記事においては上の4つのライブラリが入っているものを前提としている

その上でメタデータ、BioMartのテーブル、カウントテーブルをそれぞれ読み込む。

# メタデータ読みこみ
metadata <- read_tsv("./metadata.tsv") 

# BioMartのテーブル読込 取扱いのため列名変更
marttab <- read_tsv("./marttable.tsv") 
names(marttab) <- c("gene_id", "gene_name", "gene_description")

# カウントテーブルの読込
# バージョン情報は削除する
counttab <- read_tsv("./counttab.txt", skip =1 ) %>% 
	dplyr::rename(gene_id=Geneid) %>% 
	mutate(gene_id = str_remove(gene_id, ".\\d+$")) %>% print

以上で一旦データを読みこめた。しかしながらこのままではデータがとりあつかい辛い部分があるため、それぞれを良い形に変更していく。
一応のため、それぞれのlength等をとったテーブルを作成する。

count.ann <- counttab %>% 
  select(gene_id:Length) %>% 
  filter(!str_detect(gene_id,"_PAR_Y")) %>% 
  left_join(marttab) %>% 
  select(gene_id, gene_name, Chr:Length, gene_description) %>% 
  print

これで得られるテーブルは、Ensembl Gene IDとGene Name (Symbol)等が含まれたテーブルとなる。

次いで、メインのカウントテーブルを使い易いようにし、gencodeで入っているPAR_Yを今回は合計する。なお、この値に関しては実験に応じて取扱いが異なるそうなので注意する必要がある。以下のBiostarsの質問及びEnsemblのページを参考にしたい。

# まず読込を行い、gene_idを行名にする
suffix.bam <- ".sort.bam"
gex.counts <- counttab %>% 
  select(gene_id, ends_with(".bam")) %>% 
  column_to_rownames("gene_id") 
# 列名からsuffixを除去する
colnames(gex.counts) <- str_remove(colnames(gex.counts), suffix.bam) 

# PAR_Yを合計する処理を行う
gex.counts.enst <- gex.counts %>% 
  rownames_to_column("gene_id") %>% 
  mutate(gene_id = str_remove(gene_id, "\\..*_PAR_Y")) %>% 
  pivot_longer(
    -gene_id
  ) %>% 
  group_by(gene_id, name ) %>% 
  summarise(value = sum(value)) %>% 
  pivot_wider(id_cols = gene_id) %>% 
  print

これによって以下の様なテーブルが得られる。なお、値は適当な数字であることを留意いただきたい。

gene_id	Ctrl_1	Ctrl_2	…	TR_C_03
ENSG00000000003	2002	2502	…	5200
ENSG00000000938	0	0	…	0
︙	︙	︙	⋱	︙
ENSG00000067113	5	10	…	520

データを縦持ちにして、多数の条件に対応し易くする。

gex.longer.enst <- gex.counts.enst %>% 
  pivot_longer(
    -gene_id,
    names_to = "sampleid"
    ) %>% 
  left_join(metadata) %>% 
  print

gene_id	sampleid	value	group
ENSG00000000003	Ctrl_1	2002	0
ENSG00000000003	Ctrl_2	2502	0
︙	︙	︙	︙
ENSG00000000003	TR_C_03	5200	C
︙	︙	︙	︙
ENSG00000067113	Ctrl_1	5	0
︙	︙	︙	︙

これでようやく準備が完了したため、edgeRでのテストを一気に実行する。

テストパターンリストと関数の作成

それぞれの組み合わせでそれぞれ実行していくと、ミスも多くなり、変数を誤ってデータのとりまちがえの可能性が高くなる。よって組み合わせを全部先につくり、名前をつけて処理する。

list.tests <-
  list(
    test_A = c("0", "A"),
    test_B = c("0", "B"),
    test_C = c("0", "C")
  )

このリストと上で作成した縦持ちのテーブルを入力として、テスト結果とcpmをつけたテーブルとDGE Objectを両方もったリストを出力する関数を作成する。なお、edgeRの部分に関してはedgeR公式のUser's Guideを参考にしている。

edgeR

func.edger <- function(longertab, combination){
  # テーブルを必要なサンプルだけにフィルター
  # groupの部分をfactor型に変更
  tab <- longertab %>% 
    filter(group %in% combination) %>% 
    mutate(group = factor(group, levels = combination)) 
  
  # メタデータ部分だけを抽出 
  meta <- tab %>% ungroup() %>% distinct(sampleid, group) 
  
  # 横持ちのテーブルに変換してDGEListに読み込める形に変更
  mtx <- tab %>% 
    select(gene_id, sampleid, value) %>% 
    pivot_wider(
      id_cols = gene_id, 
      names_from = sampleid, 
      values_from = value,
      values_fill = 0
    ) %>% 
    column_to_rownames("gene_id")
  
  # DGEListを作成
  dge <- DGEList(
    counts = mtx,
    group  = meta$group
  )
  
  # マニュアルの通りに実行
  design <- model.matrix(~meta$group)
  dge <- dge[filterByExpr(dge), , keep.lib.sizes = FALSE]
  dge <- calcNormFactors(dge)
  dge <- estimateDisp(dge, design)
    
  # 今回はそれぞれのサンプルに関係がないとしてglmQLFit()を実行
  fit <- glmQLFit(dge, design)
  out <- glmQLFTest(fit)
  
  # cpmテーブルの作成 最後に".cpm"を付加することでわかりやすく
  cptab <- cpm(dge, log=FALSE)
  colnames(cptab) <- str_c(colnames(cptab), ".cpm")
  cptab <- cptab %>%  
    as.data.frame() %>% 
    rownames_to_column("id")
 
  # cpmテーブルとglmQLFTestを結合させる 
  # FDRをp.adjust()関数で計算
  # 手法は Benjamini-Hochberg 法 とする
  tabout<-out$table %>% 
    as.data.frame() %>% 
    rownames_to_column("id") %>% 
    mutate(padj = p.adjust(PValue , method = "BH")) %>% 
    left_join(cptab) %>% 
    as_tibble()
  
  # 返り値はlistとする
  return(list(table = tabout, dgeobj = dge))
}

上の関数をmap関数をつかって一気に実行する。

list.dges <- map(list.tests, ~{func.edger(gex.longer.enst, .x)})

この結果を、例えばlist.dges$test_A$tableなどで呼び出せばそれぞれを見ることができるが、ここでは一気にボルケーノプロットを描画し、テーブルを出力する。

# ボルケーノプロットの描画
map2(
  list.dges,
  names(list.dges), #ファイル名につけるために名前を呼び出す
  ~{
    .x$table %>%
      ggplot(aes(
         x=logFC, 
         y=-log10(PValue),
         color= padj < 0.05
       ))+
       geom_point()+
       labs(title = str_c("VP: ", .y))+
       theme_bw()
    
    ggsave(str_c("vp_gene_", .y, ".png"),
           width = 10, height= 7 
           )
       }
)

# テーブルの出力
map2(
  list.dges,
  names(list.dges),
  ~{
    tab <- .x$table 
    filename <- str_c("table_gene_test_", .y, ".tsv")
    
    tab %>% 
      left_join(marttab, by = c("id"="gene_id")) %>% 
      select(id, gene_name, logFC:padj, ends_with(".cpm")) %>% 
      write_tsv(filename, na = ".")
  }
)

これにて全てでの例におけるデータが得られた。ここからは解析に応じて非常に変わると考えられるため、この記事はここで終了となる。

蛇足: データの概観をPCA + Plotlyで見る

本来は最初にやるべきかもしれないが、データの全体像をつかむのにPCAを行ってみるのは悪くないかもしれない。これは追加のメモとして Plotly を用いて、インタラクティヴなプロットを作成する例である。
まず、最初の横持ちのデータをDGEListに格納、フィルターを行いcpmを計算する。なお、この際には全てグループは考えず計算しておく。その上でpcacomp()関数をつかってPCAを実行する。

dgeall <- DGEList(gex.counts.enst %>% column_to_rownames("gene_id"))
dgeall <- dgeall[filterByExpr(dgeall), , keep.lib.sizes=FALSE]
dgeall <- calcNormFactors(dgeall)
gex.cpm.all <- cpm(dgeall) %>% as.data.frame() 
cpmtab.pcacomp <- prcomp(datawizard::data_rotate(gex.cpm.all), scale = T) 
summary(cpmtab.pcacomp)

PCAの結果をもってPlotlyにて描画する。

fig <- cpmtab.pcacomp$x[, 1:2] %>% 
  as.data.frame() %>% 
  rownames_to_column("sampleid") %>% 
  left_join(metadata) %>% 
  mutate(group.vis = str_c("group_", group)) %>% 
  plotly::plot_ly(
    data = ., 
    x = ~PC1,
    y = ~PC2 ,
    text = ~sampleid,
    color = ~group.vis,
    colors = c("#888888",
                RColorBrewer::brewer.pal(4, "Set1"),
                RColorBrewer::brewer.pal(3, "Set2")
    ) 
  )
#htmlwidgets::saveWidget(fig, "PCAplot_interactive.html")
fig

おわりに

この様にして、多種の組み合わせを一気に実行する例を提示した。これはどちらかというと他の人がどの様に行っているか気になったため記事にした部分が大きい。ご意見頂けると幸いである。
なお、以上のコードについては実際に実行した例を編集している。その実行時はR 4.3.0, tidyverse 2.0.0, edgeR 3.42.4である。なお、編集した部分は主だって関数以外の部分が顕著であるため、参考にする場合には注意いただきたい。

実際のところ、ここまで問題が分けられていればChatGPTに聴いてしまえば同じような結果が得られるかもしれません。 ↩