参考资料链接

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)})

image.png

欢迎大家关注我的公众号

小明的数据分析笔记本

小明的数据分析笔记本 公众号 主要分享：1、R语言和python做数据分析和数据可视化的简单小例子；2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记；3、生物信息学入门学习资料及自己的学习笔记！