geo差异表达分析_什么，你一定要基于FPKM标准化表达矩阵做单细胞差异分析

学徒和学员已经陆续出师，是时候把生信技能树的舞台交给后辈了！

下面是《生信入门第7期》学员的分享

最近看到有一个简单的数据挖掘文章，就做了一个单细胞表达矩阵的差异分析，然后强行解释一波。而且使用的是基于FPKM标准化表达矩阵载入seurat包进行分析，就随便使用我们学习班获得的技能复现一下，分享给广大读者。

• 1、概述

  • 计算基因长度

  • 1.1 单细胞差异分析pipeline

  • 1.2 count标准化

• 2、官方fpkm数据差异分析

  • 2.1 表达矩阵与分组信息

  • 2.2 ID转换

  • 2.3 创建seurat，质控，差异分析一键操作

  • 2.4 差异结果可视化

• 3、根据count矩阵转换fpkm并完成差异分析

  • 3.1 导入count矩阵

  • 3.2 计算fpkm矩阵

  • 3.3 差异分析

  • 3.4 可视化

前言：使用GSE81861提供的数据，比较CRC肿瘤上皮细胞与正常上皮细胞的差异。

GEO提供了count与fpkm两种数据。笔记内容先用官方的fpkm数据做差异分析，再利用counts数据手动计算fpkm矩阵，完成差异分析。最后比较两种方法的结果是否存在差异。

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE81861

注：因为有不少重复的步骤，故设置较多的函数。

1、概述

1.1 单细胞差异分析pipeline

简单来说分为三步：首先导入、制备规范的表达矩阵以及分组信息；然后利用Seurat包构建seurat对象，归一化；最后进行差异分析，以及结果的可视化。

1.2 count标准化

主要受测序文库(样本总read数)与基因长度的影响，测序的counts数据不能直接进行差异分析，需要进行标准化处理。常见的几种标准化方法简单介绍如下–

• rpkm：counts先对测序文库标准化，再对基因长度标准化；

• fpkm：FPKM同RPKM是一样的，只是RPKM用于单末端测序，而FPKM用于双末端测序；

• tpm：counts先对基因长度标准化，再对测序文库标准化；

• cpm：counts只对测序文库标准化。

测序文库相对容易计算，直接使用colSums()函数即可；而基因长度则比较难求，首先要了解基因长度有不同的定义标准，其次要知道哪些R包提供相关生物数据。我目前了解到了以下三种方法，以及根据与官方fpkm验证，最终选择第三种方法用于后续的分析。

计算基因长度

(1)TxDb.Hsapiens.UCSC.hg19.knownGene包

if (!require("TxDb.Hsapiens.UCSC.hg19.knownGene", quietly = TRUE))
  BiocManager::install("TxDb.Hsapiens.UCSC.hg19.knownGene")
txdb exon_txdb=exons(txdb)
genes_txdb=genes(txdb)

g_l.1–根据非冗余外显子之和定义

g_l_1 function(){
    
  o   t1=exon_txdb[queryHits(o)]
  t2=genes_txdb[subjectHits(o)]
  t1=as.data.frame(t1)
  t1$geneid=mcols(t2)[,1]
# 得到exon_id与geneid的对应关系
  g_l.1 function(x){
    
#按gene id拆分表格
    head(x)
    tmp=apply(x,1,function(y){
    
      y[2]:y[3]
    }) #根据每一个gene所有exon的区间，生成区间内的整数，返回的为list。
    length(unique(unlist(tmp)))
#计算共有多少种整数，即为最终的非冗余exon长度之和
  })
  head(g_l.1) #为一个list
  g_l.1=data.frame(gene_id=names(g_l.1),length=as.numeric(g_l.1))
  dim(g_l.1)
  head(g_l.1)
#为基因ID增添ENSEMBLE ID
library(org.Hs.eg.db)
  s2g=toTable(org.Hs.egENSEMBL)
  head(s2g)
  g_l.1=merge(g_l.1,s2g,by='gene_id')
#把g_l,s2g两个数据框以'gene_id'为连接进行拼接
  head(g_l.1)
return(g_l.1)
}



g_l.1 head(g_l.1)

##   gene_id length      ensembl_id
## 1       1   7255 ENSG00000121410
## 2      10   1317 ENSG00000156006
## 3     100   1532 ENSG00000196839
## 4    1000   4570 ENSG00000170558
## 5   10000   7458 ENSG00000117020
## 6   10000   7458 ENSG00000275199

g_l.2—-根据最长转录本定义

g_l_2 function(){
    
  t_l=transcriptLengths(txdb)
  head(t_l)
  t_l=na.omit(t_l)
#先按基因ID，再按转录本长度从大到小排序
  t_l=t_l[order(t_l$gene_id,t_l$tx_len,decreasing = T),]
  head(t_l);dim(t_l)
#根据gene_id去重，选择第一个，也就是最长的那个
  t_l=t_l[!