学徒和学员已经陆续出师，是时候把生信技能树的舞台交给后辈了！下面是《生信入门第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

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

#按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

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

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

if (!require('biomaRt', quietly = TRUE))

BiocManager::install('biomaRt')

ensembl

ensembl = useDataset('hsapiens_gene_ensembl',mart=ensembl)

#use the hsapiens(人类) dataset.或者直接如下设置

#ensembl = useMart('ensembl',dataset='hsapiens_gene_ensembl')

test

'end_position','ensembl_transcript_id',

'transcript_length'),

mart = ensembl)

test

decreasing = T),]

g_l.3

g_l.3

head(g_l.3)

return(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

colnames(g_l.1)

g_l.2

colnames(g_l.2)

colnames(g_l.3)

g_l_all

g_l_all

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

2、官方fpkm数据差异分析

2.1 表达矩阵与分组信息#表达矩阵

nm_FPKM

data_input

#data

row.names(data)

data

data

rownames(data)

rownames(data)

colnames(data)

data

return(data)

}

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

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

dim(exp_FPKM) # 56380个基因 432个样本

## [1] 56380 432#分组信息

group_dat

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

print(dim(data))

library(org.Hs.eg.db)

ids1

ensembl_id=rownames(data))

ids2

toTable(org.Hs.egSYMBOL),by='gene_id')

ids

data

rownames(data)

print(dim(data))

return(data)

}

exp_FPKM

## [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

library(Seurat)

print('创建seurat对象...')

scRNA

meta.data=group)

print('质控...')

scRNA[['percent.mt']]

minGene=500;maxGene=4000;pctMT=15

scRNA minGene & nFeature_RNA < maxGene & percent.mt < pctMT)

print('归一化...')

scRNA

print('差异分析...')

diff_dat

group.by='group')

}

FPKM_diff

## [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 0.8,]

dim(FPKM_diff)## [1] 122 5

exp_FPKM_diff

2.4 差异结果可视化

热图my_heatmap

n

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

p.exp_FPKM_diff

箱图

my_boxplot

library(ggpubr)

if(!is.numeric(i)){

i=match(i,rownames(data))

}

df

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

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

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

####根据有基因长度的基因，筛选矩阵子集

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

head(lengths)

#### 计算样本文库大小，以及最后的fpkm计算

counts

counts[1:4,1:4]

total_count

head(total_count)

#根据counts、length、total_count计算fpkm

FPKM

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

g_l = g_l)

colnames(nm_my_FPKM)

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

g_l = g_l)

dim(tumor_FPKM)## [1] 56380 272

colnames(tumor_my_FPKM)

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

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

## [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

## [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 0.8,]

dim(my_FPKM_diff)## [1] 112 5

exp_my_FPKM_diff

3.4 可视化

热图my_heatmap(exp_my_FPKM_diff)

 ##### 箱图

my_boxplot(exp_my_FPKM_diff, 2)

my_boxplot(exp_my_FPKM_diff, 'PGK1')