30、第三十计 反客为主
本是客人却用主人的口气说话。后指在一定的场合下采取主动措施，以声势压倒别人。
代码框较乱，点这里阅读原文
最近看到有一个简单的数据挖掘文章，就做了一个单细胞表达矩阵的差异分析，然后强行解释一波。而且使用的是基于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 <- TxDb.Hsapiens.UCSC.hg19.knownGeneexon_txdb=exons(txdb)genes_txdb=genes(txdb)
g_l.1–根据非冗余外显子之和定义
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 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[!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–根据最长转录本定义
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()
比较三种结果的差异
dim(g_l.1);dim(g_l.2);dim(g_l.3)
## [1] 23729 3
## [1] 25159 3
## [1] 67130 2
g_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 3811
summary(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.3
2、官方fpkm数据差异分析
2.1 表达矩阵与分组信息
#表达矩阵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 160
nm_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.000000
tumor_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 0
dim(tumor_FPKM)#56380个基因 #272个样本
## [1] 56380 272
exp_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
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.000000
id_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 432
exp_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.000000
2.3 创建seurat，质控，差异分析一键操作
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-22
dim(FPKM_diff)
## [1] 1421 5
FPKM_diff <- FPKM_diff[FPKM_diff$p_val<0.01 & abs(FPKM_diff$avg_logFC)>0.8,]dim(FPKM_diff)
## [1] 122 5
exp_FPKM_diff <- exp_FPKM[match(rownames(FPKM_diff),rownames(exp_FPKM)),]
2.4 差异结果可视化
热图
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
图片
箱图
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矩阵
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 0
dim(nm_COUNT)
## [1] 56380 160
tumor_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 0
dim(tumor_COUNT)
## [1] 56380 272
3.2 计算fpkm矩阵
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.0000000
dim(nm_my_FPKM)
## [1] 50813 160
tumor_my_FPKM <- my_FPKM(counts = tumor_COUNT, g_l = g_l)dim(tumor_FPKM)
## [1] 56380 272
colnames(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.0000000
3.3 差异分析
exp_my_FPKM <- cbind(nm_my_FPKM,tumor_my_FPKM)dim(exp_my_FPKM) #转换id前
## [1] 50813 432
exp_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.0000000
exp_my_FPKM <- id_change(exp_my_FPKM)
## [1] 50813 432## [1] 24931 432
dim(exp_my_FPKM) #转换id后
## [1] 24931 432
exp_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.0000000
my_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-22
dim(my_FPKM_diff)
## [1] 1231 5
my_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 5
exp_my_FPKM_diff <- exp_my_FPKM[match(rownames(my_FPKM_diff),rownames(exp_my_FPKM)),]
3.4 可视化
热图
my_heatmap(exp_my_FPKM_diff)
图片
箱图
my_boxplot(exp_my_FPKM_diff, 2)
图片
my_boxplot(exp_my_FPKM_diff, "PGK1")
图片
