18. 差异表达基因探索

1. overView

  当我们在探索差异表达基因时,究竟在探索什么呢?我们可以探索差异表达基因的数量,用柱状图或韦恩图画出来.也可以比较对照组和实验组之间差异表达基因的表达情况,使用火山图绘制.但是注意的是,不是差异越大越有研究价值,如果当差异相近时,我们可能会更侧重fold change.因为p-value的值已经很小了.我们挑选出差异基因后,也可以绘制热图.行是样本,列是基因,中间值是表达量.

2. 导入数据

  首先需要导入所有差异表达分析的结果.从下图可以发现文件有点多,但有两个方法可以解决.

• 使用Linux语法合并为一个表格.
    可以在表头加一个gene_id.

#保留表头
head -1 DESeq2.12082.dir/genes.counts.matrix.Cs2_vs_Cs1.DESeq2.DE_results > merge.DE_results
#合并表
cat DESeq2.12082.dir/genes.counts.matrix.*.DE_results | grep -v '^sampleA' >> merge.DE_results
#grep -v 表示搜索出不带查找内容的行

  之后在RStudio导入.

deg_result = read.table(file='24.DE_analysis/merge.DE_results',header = T,sep='t')

• 使用R循环读入4个表格,然后合并成一个数据框

  发现R里很人性化的一点是,你合并数据框的时候,表头会自动合并.bind_rows是tidyverse提供的合并函数.

#获取文件路径
#返回的是字符串
deg_files = list.files(path='24.DE_analysis/DESeq2.12082.dir',pattern='DE_results$',full.names = T)
#
deg_test = NA
for (i in 1:length(deg_files)){
  if (i == 1) {
    deg_test = read.table(file=deg_files[i],header = T,sep='t')  
  }else{
    deg_test = bind_rows(deg_test,read.table(file=deg_files[i],sep='t'))
  }
}
#deg_test
deg_test = rownames_to_column(deg_test,var='gene_id')

3. 差异表达数据预处理

  这里主要是设置一个新列,标记是差异表达基因还是不是,如果是,需要进一步标注是上调基因还是下调基因.同样可以筛选一些列.

#select 选择列
select(deg_result,gene_id,sampleA,sampleB,log2FoldChange,pvalue,padj) %>% 
#mutate 新增列标注差异表达基因
  mutate(direction = if_else(abs(log2FoldChange)  0.05,'NS',
  #嵌套if esle
                             if_else(log2FoldChange >= 1,'Up','Down')
                             )) -> deg_result
deg_result

4. 统计差异基因数目

4.1 代码实现

  比如,我们可以按新加的列,将差异表达基因筛选出来.

filter(deg_result,direction != 'NS')
#Description:df [11,383 × 7]

  也可以按分组统计,每个样本相对于Cs1有多少差异表达基因.

filter(deg_result,direction != 'NS') %>%
  group_by(sampleA) %>% 
  summarise(count=n())

  也可以按两个条件进行分类.

filter(deg_result,direction != 'NS') %>%
  group_by(sampleA,direction) %>% 
  summarise(count=n())

4.2 韦恩图

  用韦恩图绘制前,需要先引入包.R里可用于画韦恩图的包很多,这里只举例一个

BiocManager::install('VennDiagram')
library(VennDiagram)

  韦恩图用于画各个样本之间的重复差异表达基因的数目.因此我们需要将各个样本之间差异表达基因的ID列出来.下面举一个group2group2group2和group1group1group1之间的例子.

#方法二:select(filter(deg_result,sampleA=='Cs2',sampleB=='Cs1',direction!='NS'),gene_id)
filter(deg_result,sampleA=='Cs2',sampleB=='Cs1',direction!='NS')$gene_id
#方法三:
filter(deg_result,sampleA=='Cs2',sampleB=='Cs1',direction!='NS') %>%
    pull(gene_id) -> Cs2#将该列提取出来

  接下来,处理好Cs3~5.
  VennDiagram::venn.diagram函数可以用于绘制韦恩图.它需要列表做输入.R的列表非常奇怪,可以容纳各种类型的数据.甚至可以给每个列表数据起一个名字.

l = list(a=2,b='23131',c=3.13213)
l
l['a']
l[1]#也可以提取2

  那么我们将差异基因的gene_id打包放在list.

deg_list = list(
  Cs2 = pull(filter(deg_result,sampleA=='Cs2',sampleB=='Cs1',direction!='NS'),gene_id),
  Cs3 = pull(filter(deg_result,sampleA=='Cs3',sampleB=='Cs1',direction!='NS'),gene_id),
  Cs4 = pull(filter(deg_result,sampleA=='Cs4',sampleB=='Cs1',direction!='NS'),gene_id),
  Cs5 = pull(filter(deg_result,sampleA=='Cs5',sampleB=='Cs1',direction!='NS'),gene_id)
)
deg_list
#绘制韦恩图
venn.diagram(x = deg_list,filename = 'deg_venn.tiff')
#此包实现画图不允许图直接出现在输出下.

  但是可以取文件路径找.

  也可以指定颜色.

library(RColorBrewer)
venn.diagram(x = deg_list,filename = 'deg_venn.tiff',fill=RColorBrewer::brewer.pal(4,'BrBG'))

5. 火山图

  可以画火山图的包也有很多,这里举例EnhancedVolcano.火山图只能展示一组比较,这里就拿Cs2和Cs1比.

library(EnhancedVolcano)

  提取Cs2相对于Cs1的差异结果.然后绘制火山图.会直接输出在Rmd里,比较难直接观察.可以考虑在命令行输入指令,然后从plots面板导出,注意pdf长宽不要默认,调大一点,比如15x15.

deg_vol = filter(deg_result,sampleA=='Cs2',sampleB=='Cs1')
EnhancedVolcano(
  #差异基因表格
  toptable = filter(deg_result,sampleA=='Cs2',sampleB=='Cs1'),
  #火山图点标签
  lab = gene_id,
  #表格里log2foldchange名
  x = 'log2FoldChange',
  y = 'padj',
)

  灰色的是非差异表达基因.红色是同时满足2个条件的差异表达基因.同时有一个很重要的参数selectlab可以选择要显示的标签.

key_gene = c('Cs_ont_9g023780','Cs_ont_9g002090','Cs_ont_7g024130','Cs_ont_1g026520')
EnhancedVolcano(
  #差异基因表格
  toptable = filter(deg_result,sampleA=='Cs2',sampleB=='Cs1'),
  #火山图点标签
  lab = filter(deg_result,sampleA=='Cs2',sampleB=='Cs1')$gene_id,
  #表格里log2foldchange名
  x = 'log2FoldChange',
  y = 'padj',
  selectLab = key_gene,
  pcutoff = 0.05,#padj的边界
  FCcutoff = 2,#log2FoldChange的边界
)

6. 热图

  差异表达基因的关系可视化还有一种方法是绘制热图,我们可以看这些基因的关系性.这里我们只从Cs2组挑出前20个genes.

filter(deg_result,sampleA=='Cs2',sampleB=='Cs1') %>% 
  arrange(desc(abs(log2FoldChange))) -> deg_sort
slice(deg_sort,1:20)  

  接下来,需要根据这些GenesGenesGenes在表达矩阵里选出对应行.

key_genes = slice(deg_sort,1:20) %>% pull(gene_id)
gen_exp = read.table(file='23.merge/genes.TMM.TPM.matrix')
#挑选出对应行,一定要加, 这是R语言基础
key_gen_exp = gen_exp[key_genes,]
key_gen_exp

  之后调用函数即可.之前我们画的是样本之间相关系数矩阵,这次是基因与样本的表达量矩阵.

library(pheatmap)
pheatmap(key_gen_exp,cluster_rows = F,cluster_cols = F)

  当我们测到的表达矩阵,如果存在基因在样本A表达量很高很高,而其他基因在样本B表达量很低很低,会导致热图的颜色映射区间很大,比如上图,因为最高值太高,导致几乎全是蓝色.

  第一个解决这个问题的方法就是取对数.

pheatmap(log10(key_gen_exp+1),cluster_rows = F,cluster_cols = F)

  还有一种方法是进行标准化.这里标准化的对象是一行,也就是genes.标准化后,就不能进行基因与基因之间的比较了.
热图是可以聚类的,只是上节不应该对相关性进行聚类,但是可以和表达量聚类

pheatmap(key_gen_exp+1,scale = 'row',cluster_rows = F,cluster_cols = F)

  也可以继续细化图.比如下面的参数:

• cutree_rows:将行分为x簇
• cutree_columns:将列分为x簇