きまぐれ研究室

座右の名は、"let it be"ですかね

Rでマルチプルアライメントを可視化しよう

マルチプルアライメントを可視化するRパッケージに、ggmsaを見つけた。

yulab-smu.top

可視化するために、気に入った見た目にできるツールがあまりなかったので良かった。

(雑な)使い方:とりあえず使ってみる

DNAの時

library(ggmsa)
library(ggplot2)

msa_file <- "file_path" #読み込むfastaファイルのファイル名

ggmsa(msa_file, font = NULL, color = "Chemistry_NT", show.legend = TRUE) + #font = NULLにすると、配列の文字が消えるので、動作が速くなる。
    facet_msa(field = 120) + #横に入りきらないとき、決まった長さごとに改行できる。
    labs(fill = NULL) + #レジェンド名を消す。
    theme(text = element_text(family = "helvetica", colour = "black", size = 10), axis.text = element_text(family = "helvetica", colour = "black", size = 5)) #文字サイズなどの調節。
ggsave("figure_name.png", width = 8, height = 8, dpi = 600, device = "png", bg="white") #サイズなどは適宜変更

アミノ酸の時（ほぼ一緒）

library(ggmsa)
library(ggplot2)

msa_file <- "file_path" #読み込むfastaファイルのファイル名

ggmsa(msa_file, font = NULL, color = "Zappo_AA", show.legend = TRUE) + #font = NULLにすると、配列の文字が消えるので、動作が速くなる。
    facet_msa(field = 120) + #横に入りきらないとき、決まった長さごとに改行できる。
    labs(fill = NULL) + #レジェンド名を消す。
    theme(text = element_text(family = "helvetica", colour = "black", size = 10), axis.text = element_text(family = "helvetica", colour = "black", size = 5)) #文字サイズなどの調節。
ggsave("figure_name.png", width = 8, height = 8, dpi = 600, device = "png", bg="white") #サイズなどは適宜変更

colorやfontは、好みによって変えること。

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

生物学 R シロイヌナズナ

注意

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

やりたいこと

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コードと発現量と機能が一つの表で見られるようになっているはず。

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