シロイヌナズナでedgeRの結果に、TAIR10のfunctional descriptionを加えるには？

注意

一応うまく行っているか確認しながらやってください。

やりたいこと

edgeRで出力した発現変動遺伝子の結果に、TAIR10からダウンロードした"functional description"を加えて一つの表にする。

手順

1. edgeRで発現変動遺伝子を取り出す。

サンプルコード

library(edgeR)
library(tidyverse)

# カウントデータを読み込む（列名をAGIとして読み込む。ただし、トランスクリプトは無視する。）
data_count <- read_csv("count_data.csv", col_names = c("A1", "A2", "A3", "B1", "B2", "B3")) |> column_to_rownames(var = "gene_id")

group = factor(c("A", "A", "A", "B", "B", "B"))
time = factor(c("1", "2", "3", "1", "2", "3"))

data_count_DGE <- DGEList(counts = data_count, group = group)
keep <- filterByExpr(data_count_DGE)
data_count_DGE <- data_count_DGE[keep, , keep.lib.sizes = FALSE] |> normLibSizes()


design <- model.matrix(~species)
rownames(design) <- colnames(data_count_DGE)
estimate <- estimateDisp(data_count_DGE, design, robust = TRUE)
fit <- glmQLFit(estimate, design, robust = TRUE)

qlf <- glmQLFTest(fit)
data_qlf <- qlf |> topTags(n = nrow(data_count)) |> as.data.frame() |> rownames_to_column(var = "gene_id")

data_cpm <- cpm(data_count_DGE) |> as.data.frame() |> rownames_to_column(var = "gene_id") |> full_join(data_qlf)

2. TAIRからTAIR10_functional_descriptionsをダウンロード

TAIR10_functional_descriptions

3. TAIR10_functional_descriptionsをRで読み込む

# ただし、トランスクリプトを無視する場合。遺伝子名を1文字目から9文字目までを取り出している。
description <- read_tsv("./TAIR10_functional_descriptions.txt") |>
    rename("gene_id" = Model_name) |>
    mutate(gene_id = str_sub(gene_id, end = 9)) |>
    distinct(gene_id, .keep_all = TRUE)

4. 2つのデータを結合する。

data_cpm <- data_cpm |> full_join(description, by = "gene_id")
data_cpm |> write_csv("data_cpm.csv")

data_cpm.csvを見ると、AGIコードと発現量と機能が一つの表で見られるようになっているはず。

間違いあったらコメントください。