首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >R语言利用转录组基因表达矩阵做基因共表达分析的学习资料推荐

R语言利用转录组基因表达矩阵做基因共表达分析的学习资料推荐

作者头像
用户7010445
发布2023-01-06 20:29:26
3240
发布2023-01-06 20:29:26
举报

参考资料链接

https://github.com/cxli233/SimpleTidy_GeneCoEx/tree/v1.0.1

提供完整的示例数据和代码,非常好的学习材料

做基因共表达比较常用的是WGCNA那个R包,这个链接里提供的代码不是用WGCNA这个R包实现的,而是利用表达量数据计算不同基因之间的相关性,这种方法也挺常用的在论文里见过

表达量数据是来源于论文

High-resolution spatiotemporal transcriptome mapping of tomato fruit development and ripening

https://www.nature.com/articles/s41467-017-02782-9

数据是不同发育阶段的转录数据,表达量数据的下载链接是 https://zenodo.org/record/7117357#.Y0WB13ZBzic

关于样本的一些分组信息在链接里提供了,大家如果感兴趣可以自己下载数据然后跟着这个链接完全重复一下

接下来的内容我重复一下资料中利用表达量数据做PCA的内容

代码

setwd("data/20221012/")
list.files()

#library(data.table)
library(readr)
Exp_table <- read_csv("Shinozaki_tpm_all.csv", col_types = cols())
head(Exp_table)
dim(Exp_table)

library(readxl)
Metadata <- read_excel("SimpleTidy_GeneCoEx-1.0.1/Data/Shinozaki_datasets_SRA_info.xlsx")
head(Metadata)
dim(Metadata)

library(tidyverse)
Exp_table_long <- Exp_table %>% 
  rename(gene_ID = `...1`) %>% 
  pivot_longer(cols = !gene_ID, names_to = "library", values_to = "tpm") %>% 
  mutate(logTPM = log10(tpm + 1)) 

head(Exp_table_long)


Exp_table_log_wide <- Exp_table_long %>% 
  select(gene_ID, library, logTPM) %>% 
  pivot_wider(names_from = library, values_from = logTPM)

head(Exp_table_log_wide)

my_pca <- prcomp(t(Exp_table_log_wide[, -1]))
pc_importance <- as.data.frame(t(summary(my_pca)$importance))
head(pc_importance, 20)


PCA_coord <- my_pca$x[, 1:10] %>% 
  as.data.frame() %>% 
  mutate(Run = row.names(.)) %>% 
  full_join(Metadata %>% 
              select(Run, tissue, dev_stage, `Library Name`, `Sample Name`), by = "Run")

head(PCA_coord)


PCA_coord <- PCA_coord %>% 
  mutate(stage = case_when(
    str_detect(dev_stage, "MG|Br|Pk") ~ str_sub(dev_stage, start = 1, end = 2),
    T ~ dev_stage
  )) %>% 
  mutate(stage = factor(stage, levels = c(
    "Anthesis",
    "5 DPA",
    "10 DPA",
    "20 DPA",
    "30 DPA",
    "MG",
    "Br",
    "Pk",
    "LR",
    "RR"
  ))) %>% 
  mutate(dissection_method = case_when(
    str_detect(tissue, "epidermis") ~ "LM",
    str_detect(tissue, "Collenchyma") ~ "LM",
    str_detect(tissue, "Parenchyma") ~ "LM",
    str_detect(tissue, "Vascular") ~ "LM",
    str_detect(dev_stage, "Anthesis") ~ "LM",
    str_detect(dev_stage, "5 DPA") &
      str_detect(tissue, "Locular tissue|Placenta|Seeds") ~ "LM",
    T ~ "Hand"
  ))

head(PCA_coord)

library(viridis)
library(RColorBrewer)

PCA_by_method <- PCA_coord %>% 
  ggplot(aes(x = PC1, y = PC2)) +
  geom_point(aes(fill = dissection_method), color = "grey20", shape = 21, size = 3, alpha = 0.8) +
  scale_fill_manual(values = brewer.pal(n = 3, "Accent")) +
  labs(x = paste("PC1 (", pc_importance[1, 2] %>% signif(3)*100, "% of Variance)", sep = ""), 
       y = paste("PC2 (", pc_importance[2, 2] %>% signif(3)*100, "% of Variance)", "  ", sep = ""),
       fill = NULL) +  
  theme_bw() +
  theme(
    text = element_text(size= 14),
    axis.text = element_text(color = "black")
  )


PCA_by_method

PCA_by_tissue <- PCA_coord %>% 
  ggplot(aes(x = PC1, y = PC2)) +
  geom_point(aes(fill = tissue), color = "grey20", shape = 21, size = 3, alpha = 0.8) +
  scale_fill_manual(values = brewer.pal(11, "Set3")) +
  labs(x = paste("PC2 (", pc_importance[2, 2] %>% signif(3)*100, "% of Variance)", sep = ""), 
       y = paste("PC3 (", pc_importance[3, 2] %>% signif(3)*100, "% of Variance)", "  ", sep = ""),
       fill = "tissue") +  
  theme_bw() +
  theme(
    text = element_text(size= 14),
    axis.text = element_text(color = "black")
  )

PCA_by_tissue

PCA_by_stage <- PCA_coord %>% 
  ggplot(aes(x = PC2, y = PC3)) +
  geom_point(aes(fill = stage), color = "grey20", shape = 21, size = 3, alpha = 0.8) +
  scale_fill_manual(values = viridis(10, option = "D")) +
  labs(x = paste("PC2 (", pc_importance[2, 2] %>% signif(3)*100, "% of Variance)", sep = ""), 
       y = paste("PC3 (", pc_importance[3, 2] %>% signif(3)*100, "% of Variance)", "  ", sep = ""),
       fill = "stage") +  
  theme_bw() +
  theme(
    text = element_text(size= 14),
    axis.text = element_text(color = "black")
  )

PCA_by_stage 


library(patchwork)

PCA_by_method+PCA_by_tissue+PCA_by_tissue

image.png

以上用到的代码和示例数据都可以在推文开头提到链接里找到

上面的代码有一步是对TPM值 加1然后取log10,他的实现方式是先将宽格式数据转换为长格式,然后把取log10后的长格式再转换为宽格式,这里我没能还可以借助mutate_at()函数

Exp_table %>% select(1,2,3) %>% 
  rename("gene_id"="...1") %>% 
  mutate_at(vars(starts_with("SRR")),
            function(x){log10(x+1)})
本文参与 腾讯云自媒体分享计划,分享自微信公众号。
原始发表:2022-10-12,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 小明的数据分析笔记本 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体分享计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 参考资料链接
  • 代码
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档