下面是《生信入门第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包
代码语言：javascript复制if (!require("TxDb.Hsapiens.UCSC.hg19.knownGene", quietly = TRUE))
BiocManager::install("TxDb.Hsapiens.UCSC.hg19.knownGene")
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
exon_txdb=exons(txdb)
genes_txdb=genes(txdb)g_l.1–根据非冗余外显子之和定义
代码语言：javascript复制g_l_1 <- function(){
o <- findOverlaps(exon_txdb,genes_txdb)
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 <- lapply(split(t1,t1$geneid),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 <- 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 ENSG00000275199g_l.2—-根据最长转录本定义
代码语言：javascript复制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[!duplicated(t_l$gene_id),]
head(t_l);dim(t_l)
g_l.2=t_l[,c(3,5)]
library(org.Hs.eg.db)
s2g=toTable(org.Hs.egENSEMBL)
g_l.2=merge(g_l.2,s2g,by='gene_id')
head(g_l.2)
return(g_l.2)
}
g_l.2 <- g_l_2()
head(g_l.2)## gene_id tx_len ensembl_id
## 1 1 3419 ENSG00000121410
## 2 10 1317 ENSG00000156006
## 3 100 1532 ENSG00000196839
## 4 1000 4367 ENSG00000170558
## 5 10000 7082 ENSG00000275199
## 6 10000 7082 ENSG00000117020(2)biomaRt包
g_l.3–根据最长转录本定义
代码语言：javascript复制g_l_3 <- function(){
if (!require("biomaRt", quietly = TRUE))
BiocManager::install("biomaRt")
ensembl <- useMart("ensembl") #connect to a specified BioMart database
ensembl = useDataset("hsapiens_gene_ensembl",mart=ensembl)
#use the hsapiens(人类) dataset.或者直接如下设置
#ensembl = useMart("ensembl",dataset="hsapiens_gene_ensembl")
test <- getBM(attributes=c('ensembl_gene_id', 'start_position',
'end_position','ensembl_transcript_id',
'transcript_length'),
mart = ensembl)
test <- test[order(test$ensembl_gene_id,test$transcript_length,
decreasing = T),]
g_l.3 <- test[!duplicated(test$ensembl_gene_id),]
g_l.3 <- g_l.3[,c(1,5)]
head(g_l.3)
return(g_l.3)
}
g_l.3 <- g_l_3()比较三种结果的差异
代码语言：javascript复制dim(g_l.1);dim(g_l.2);dim(g_l.3)## [1] 23729 3## [1] 25159 3## [1] 67130 2g_l.1 <- g_l.1[,-1]
colnames(g_l.1) <- c("g_l.1","ensembl_id")
g_l.2 <- g_l.2[,-1]
colnames(g_l.2) <- c("g_l.2","ensembl_id")
colnames(g_l.3) <- c("ensembl_id","g_l.3")
g_l_all <- merge(g_l.1, g_l.2, by="ensembl_id")
g_l_all <- merge(g_l_all, g_l.3, by="ensembl_id")
head(g_l_all,10)## ensembl_id g_l.1 g_l.2 g_l.3
## 1 ENSG00000000003 2486 2069 3796
## 2 ENSG00000000005 3031 3031 1205
## 3 ENSG00000000419 1069 1069 1161
## 4 ENSG00000000457 3426 3090 6308
## 5 ENSG00000000460 5654 3852 4355
## 6 ENSG00000000938 3470 2485 2637
## 7 ENSG00000000971 4386 4127 6985
## 8 ENSG00000001036 4415 2749 2385
## 9 ENSG00000001084 6221 3812 3785
## 10 ENSG00000001167 6589 6385 3811summary(g_l_all)## ensembl_id g_l.1 g_l.2 g_l.3
## Length:23906 Min. : 20 Min. : 20 Min. : 27
## Class :character 1st Qu.: 1636 1st Qu.: 1472 1st Qu.: 1660
## Mode :character Median : 2964 Median : 2515 Median : 2902
## Mean : 3760 Mean : 3083 Mean : 3564
## 3rd Qu.: 4934 3rd Qu.: 4104 3rd Qu.: 4743
## Max. :117846 Max. :109223 Max. :109224#最终选择第三种结果
g_l <- g_l.32、官方fpkm数据差异分析
2.1 表达矩阵与分组信息
代码语言：javascript复制#表达矩阵
nm_FPKM <- read.csv("GSE81861_CRC_NM_epithelial_cells_FPKM.csv/GSE81861_CRC_NM_epithelial_cells_FPKM.csv")
data_input <- function(data){
#data <- nm_FPKM
row.names(data) <- data[,1]
data <- data[,-1]
data <- as.matrix(data)
rownames(data) <- sapply(strsplit(rownames(data),"_"),"[",3)
rownames(data) <- sapply(strsplit(rownames(data),"[.]"),"[",1)
colnames(data) <- sapply(strsplit(colnames(data),"_"),"[",1)
data <- data[grep("ENSG",rownames(data)),]
return(data)
}
nm_FPKM <- data_input(nm_FPKM)
dim(nm_FPKM) #56380个基因 #160个样本## [1] 56380 160nm_FPKM[1:4,1:4]## RHC4104 RHC6087 RHL2880 RHC5949
## ENSG00000000003 0.000 185.315 323.203 22.311700
## ENSG00000000005 0.000 0.000 0.000 0.000000
## ENSG00000000419 0.000 231.756 0.000 0.851514
## ENSG00000000457 100.135 0.000 0.000 0.000000tumor_FPKM <- read.csv("GSE81861_CRC_tumor_epithelial_cells_FPKM.csv/GSE81861_CRC_tumor_epithelial_cells_FPKM.csv")
tumor_FPKM <- data_input(tumor_FPKM)
tumor_FPKM[1:4,1:4]## RHC4075 RHC5563 RHC5552 RHC4874
## ENSG00000000003 0 0.000 0 0
## ENSG00000000005 0 0.000 0 0
## ENSG00000000419 0 0.000 0 0
## ENSG00000000457 0 193.167 0 0dim(tumor_FPKM)#56380个基因 #272个样本## [1] 56380 272exp_FPKM <- cbind(nm_FPKM,tumor_FPKM)
dim(exp_FPKM) # 56380个基因 432个样本## [1] 56380 432#分组信息
group_dat <- data.frame(group=c(rep('normal',160),
rep('tumor',272)),
row.names = colnames(exp_FPKM))2.2 ID转换
因为seurat质控需要过滤线粒体基因，所以需要把ensembl ID转换为 symbol ID
代码语言：javascript复制exp_FPKM[1:4,1:4]## RHC4104 RHC6087 RHL2880 RHC5949
## ENSG00000000003 0.000 185.315 323.203 22.311700
## ENSG00000000005 0.000 0.000 0.000 0.000000
## ENSG00000000419 0.000 231.756 0.000 0.851514
## ENSG00000000457 100.135 0.000 0.000 0.000000id_change <- function(data){
print(dim(data))
library(org.Hs.eg.db)
ids1 <- data.frame(ID=c(1:nrow(data)),
ensembl_id=rownames(data))
ids2 <- merge(toTable(org.Hs.egENSEMBL),
toTable(org.Hs.egSYMBOL),by="gene_id")
ids <- merge(ids1, ids2, by="ensembl_id")
data <- data[ids$ID,]
rownames(data) <- ids$symbol
print(dim(data))
return(data)
}
exp_FPKM <- id_change(exp_FPKM) #有2w 个ensembl ID没配对到 symbol## [1] 56380 432
## [1] 24941 432exp_FPKM[1:4,1:4]## RHC4104 RHC6087 RHL2880 RHC5949
## TSPAN6 0.000 185.315 323.203 22.311700
## TNMD 0.000 0.000 0.000 0.000000
## DPM1 0.000 231.756 0.000 0.851514
## SCYL3 100.135 0.000 0.000 0.0000002.3 创建seurat，质控，差异分析一键操作
代码语言：javascript复制scRNA_deg <- function(exp,group){
library(Seurat)
print("创建seurat对象...")
scRNA <- CreateSeuratObject(counts=exp,
meta.data=group)
print("质控...")
scRNA[["percent.mt"]] <- PercentageFeatureSet(scRNA, pattern = "^MT-")
minGene=500;maxGene=4000;pctMT=15
scRNA <- subset(scRNA, subset = nFeature_RNA > minGene & nFeature_RNA < maxGene & percent.mt < pctMT)
print("归一化...")
scRNA <- NormalizeData(scRNA, normalization.method = "LogNormalize", scale.factor = 10000)
print("差异分析...")
diff_dat <- FindMarkers(scRNA,ident.1="normal",ident.2="tumor",
group.by='group')
}
FPKM_diff <- scRNA_deg(exp=exp_FPKM, group=group_dat)## [1] "创建seurat对象..."
## [1] "质控..."
## [1] "归一化..."
## [1] "差异分析..."head(FPKM_diff)## p_val avg_logFC pct.1 pct.2 p_val_adj
## GUCA2A 1.353485e-30 2.576937 0.742 0.101 3.375728e-26
## CLCA4 3.106166e-30 1.112995 0.688 0.067 7.747089e-26
## MT1E 7.460402e-28 1.733571 0.789 0.162 1.860699e-23
## CKB 1.723524e-27 2.241367 0.883 0.285 4.298641e-23
## ZG16 1.936670e-27 2.753875 0.664 0.073 4.830248e-23
## PHGR1 5.482718e-27 1.882931 0.938 0.531 1.367445e-22dim(FPKM_diff)## [1] 1421 5FPKM_diff <- FPKM_diff[FPKM_diff$p_val<0.01 & abs(FPKM_diff$avg_logFC)>0.8,]
dim(FPKM_diff)## [1] 122 5exp_FPKM_diff <- exp_FPKM[match(rownames(FPKM_diff),rownames(exp_FPKM)),]2.4 差异结果可视化
热图
代码语言：javascript复制my_heatmap <- function(exp_dat){
n <- t(scale(t(exp_dat)))
n[n>2]=2;n[n< -2]= -2
library(pheatmap)
pheatmap(n, show_rownames = F,
show_colnames = F,
annotation_col = group_dat)
}
p.exp_FPKM_diff <- my_heatmap(exp_FPKM_diff)
p.exp_FPKM_diff
箱图
代码语言：javascript复制my_boxplot <- function(data,i){
library(ggpubr)
if(!is.numeric(i)){
i=match(i,rownames(data))
}
df <- data.frame(gene=data[i,], group=group_dat$group)
ggboxplot(df,x="group",y="gene",
color = "group",add = "jitter",
ylab = rownames(data)[i])
theme_bw()
}
my_boxplot(exp_FPKM_diff, 2)
my_boxplot(exp_FPKM_diff, "PGK1")
3、根据count矩阵转换fpkm并完成差异分析
3.1 导入count矩阵
代码语言：javascript复制nm_COUNT <- read.csv("GSE81861_CRC_NM_epithelial_cells_COUNT.csv/GSE81861_CRC_NM_epithelial_cells_COUNT.csv")
nm_COUNT <- data_input(nm_COUNT)
nm_COUNT[1:4,1:4]## RHC4104 RHC6087 RHL2880 RHC5949
## ENSG00000000003 0 428 66 141
## ENSG00000000005 0 0 0 0
## ENSG00000000419 0 179 0 1
## ENSG00000000457 465 0 0 0dim(nm_COUNT)## [1] 56380 160tumor_COUNT <- read.csv("GSE81861_CRC_tumor_epithelial_cells_COUNT.csv/GSE81861_CRC_tumor_epithelial_cells_COUNT.csv")
tumor_COUNT <- data_input(tumor_COUNT)
tumor_COUNT[1:4,1:4]## RHC4075 RHC5563 RHC5552 RHC4874
## ENSG00000000003 0 0 0 0
## ENSG00000000005 0 0 0 0
## ENSG00000000419 0 0 0 0
## ENSG00000000457 0 133 0 0dim(tumor_COUNT)## [1] 56380 2723.2 计算fpkm矩阵
代码语言：javascript复制my_FPKM <- function(counts,g_l){
####根据有基因长度的基因，筛选矩阵子集
#counts=nm_COUNT
ng=intersect(rownames(counts),g_l$ensembl_id)
length(ng)
lengths=g_l[match(ng,g_l$ensembl_id),2]
names(lengths) <- g_l[match(ng,g_l$ensembl_id),1]
head(lengths)
#### 计算样本文库大小，以及最后的fpkm计算
counts <- counts[names(lengths),]
counts[1:4,1:4]
total_count <- colSums(counts)
head(total_count)
#根据counts、length、total_count计算fpkm
FPKM <- t(do.call( rbind,
lapply(1:length(total_count),
function(i){
10^9*counts[,i]/lengths/total_count[i]
#lengths向量自动遍历
}) ))
FPKM[1:4,1:4]
return(FPKM)
}
nm_my_FPKM <- my_FPKM(counts = nm_COUNT,
g_l = g_l)
colnames(nm_my_FPKM) <- colnames(nm_COUNT)
nm_my_FPKM[1:4,1:4]## RHC4104 RHC6087 RHL2880 RHC5949
## ENSG00000000003 0.0000 160.5715 168.6839 22.7021034
## ENSG00000000005 0.0000 0.0000 0.0000 0.0000000
## ENSG00000000419 0.0000 219.5694 0.0000 0.5264304
## ENSG00000000457 124.2016 0.0000 0.0000 0.0000000dim(nm_my_FPKM)## [1] 50813 160tumor_my_FPKM <- my_FPKM(counts = tumor_COUNT,
g_l = g_l)
dim(tumor_FPKM)## [1] 56380 272colnames(tumor_my_FPKM) <- colnames(tumor_COUNT)
nm_my_FPKM[1:4,1:4]## RHC4104 RHC6087 RHL2880 RHC5949
## ENSG00000000003 0.0000 160.5715 168.6839 22.7021034
## ENSG00000000005 0.0000 0.0000 0.0000 0.0000000
## ENSG00000000419 0.0000 219.5694 0.0000 0.5264304
## ENSG00000000457 124.2016 0.0000 0.0000 0.00000003.3 差异分析
代码语言：javascript复制exp_my_FPKM <- cbind(nm_my_FPKM,tumor_my_FPKM)
dim(exp_my_FPKM) #转换id前## [1] 50813 432exp_my_FPKM[1:4,1:4] #转换id前## RHC4104 RHC6087 RHL2880 RHC5949
## ENSG00000000003 0.0000 160.5715 168.6839 22.7021034
## ENSG00000000005 0.0000 0.0000 0.0000 0.0000000
## ENSG00000000419 0.0000 219.5694 0.0000 0.5264304
## ENSG00000000457 124.2016 0.0000 0.0000 0.0000000exp_my_FPKM <- id_change(exp_my_FPKM)## [1] 50813 432
## [1] 24931 432dim(exp_my_FPKM) #转换id后## [1] 24931 432exp_my_FPKM[1:4,1:4] #转换id后## RHC4104 RHC6087 RHL2880 RHC5949
## TSPAN6 0.0000 160.5715 168.6839 22.7021034
## TNMD 0.0000 0.0000 0.0000 0.0000000
## DPM1 0.0000 219.5694 0.0000 0.5264304
## SCYL3 124.2016 0.0000 0.0000 0.0000000my_FPKM_diff <- scRNA_deg(exp=exp_my_FPKM, group=group_dat)## [1] "创建seurat对象..."
## [1] "质控..."
## [1] "归一化..."
## [1] "差异分析..."head(my_FPKM_diff)## p_val avg_logFC pct.1 pct.2 p_val_adj
## CLCA4 2.640706e-30 1.6080767 0.688 0.067 6.583544e-26
## GUCA2A 3.289287e-30 2.4478445 0.742 0.101 8.200520e-26
## ZG16 1.629645e-27 2.9466466 0.664 0.073 4.062867e-23
## MT1E 2.707239e-27 1.4704522 0.789 0.162 6.749417e-23
## CKB 5.633316e-27 1.9409546 0.883 0.285 1.404442e-22
## PADI2 2.153941e-26 0.5522112 0.734 0.128 5.369989e-22dim(my_FPKM_diff)## [1] 1231 5my_FPKM_diff <- my_FPKM_diff[my_FPKM_diff$p_val<0.01 & abs(my_FPKM_diff$avg_logFC)>0.8,]
dim(my_FPKM_diff)## [1] 112 5exp_my_FPKM_diff <- exp_my_FPKM[match(rownames(my_FPKM_diff),rownames(exp_my_FPKM)),]3.4 可视化
热图
代码语言：javascript复制my_heatmap(exp_my_FPKM_diff)
##### 箱图
代码语言：javascript复制my_boxplot(exp_my_FPKM_diff, 2)
my_boxplot(exp_my_FPKM_diff, "PGK1")
