拖后腿学徒居然也完成作业，理解RNA-seq数据分析结果
前面我出了一个学徒作业，下载表达矩阵后绘制PCA图及热图，然后理解作者给出的RPM和raw_counts的差异，详见：理解RNA-seq表达矩阵的两个形式 很意外，12月学徒肖一僧居然吭哧吭哧的完成了，吓我一跳！让我们看看他的表演 以下是正文 收到大佬的作业，第一次投稿。大佬的题目如下：通过一篇science文章，理解两种RNA-seq表达矩阵在数据分析的时候是否相同(大佬的意思是通过PCA和heatmap来看一下)。
稍微介绍 一下背景
Counts值
对给定的基因组参考区域，计算比对上的read数，又称为raw count（RC），也就我通常说的相对原始的数据，是没进行任何标准化操作的数据。
计数结果的差异的影响因素：落在参考区域上下限的read是否需要被统计，按照什么样的标准进行统计。
RPM (Reads per million mapped reads)
RPM方法：10^6标准化了测序深度的影响，以前说这个没有考虑转录本的长度的影响。表示RPM适合于产生的read读数不受基因长度影响的测序方法，比如miRNA-seq测序，miRNA的长度一般在20-24个碱基之间。但是现在的观点（Jimmy大佬说）认为，我们通常做差异表达，同于对比系统下转录本长度是一致（例如一个患者的癌跟癌旁），所以我们只要需要考虑测序深度的影响，所以这个RPM是比较好的一种数据标准化方式。
可能优于RPKM/FPKM (Reads/Fragments per kilo base per million mapped reads)。这个网页也说明了这个问题：https://www.jianshu.com/p/35e861b76486
下面给给出我的答案及代码
一，读取数据，并做一下常规的转换
###一些常规的设置
rm(list = ls())#清空环境变量
options(stringsAsFactors = F)##字符不作为因子读入
###读取数据。
##rawcounts##
a <- read.table('GSE103788_raw_counts_genes.tsv.gz',header = T,row.names = 1)
a[1:4,1:4]##大概看一下数据格式。
### PH_WT_tumors_r1 PH_WT_tumors_r2
ENSMUSG00000000001_Gnai3 1827 2772
ENSMUSG00000000028_Cdc45 44 88
ENSMUSG00000000031_H19 19 500
ENSMUSG00000000037_Scml2 5 10
PH_WT_tumors_r3 PH_WT_notreat_r1
ENSMUSG00000000001_Gnai3 2149 3264
ENSMUSG00000000028_Cdc45 72 28
ENSMUSG00000000031_H19 39 1
ENSMUSG00000000037_Scml2 11 3
##RPM##
b <- read.table('GSE103788_filtered_RPM_genes.tsv.gz',header = T,row.names = 1)
b[1:4,1:4]
## PH_WT_tumors_r1 PH_WT_tumors_r2
ENSMUSG00000000001_Gnai3 111.4598680 149.6572460
ENSMUSG00000000028_Cdc45 2.6843099 4.7510237
ENSMUSG00000000031_H19 1.1591338 26.9944527
ENSMUSG00000000037_Scml2 0.3050352 0.5398891
PH_WT_tumors_r3 PH_WT_notreat_r1
ENSMUSG00000000001_Gnai3 120.3417347 183.67178123
ENSMUSG00000000028_Cdc45 4.0319241 1.57561577
ENSMUSG00000000031_H19 2.1839589 0.05627199
ENSMUSG00000000037_Scml2 0.6159884 0.16881598
##查看数据是否需要做log转换。
qx <- as.numeric(quantile(a, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
LogC <- (qx[5] > 100) ||
(qx[6]-qx[1] > 50 && qx[2] > 0) ||
(qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2)
LogC
##[1] TRUE
##########boxplot可视化数据检查数据是否需要log等转换
boxplot(a,las=2,cex.axis=0.6,main='data check')
###去除低碱基质量的基因。
a=a[apply(a,1, function(x) sum(x>1) > 6),]
ex<- log2(a+1)##转换
boxplot(ex,las=2,cex.axis=0.6,main='data check')
下图是转换前后的图片。
同上方法，处理一下RPM的数据。
##########boxplot可视化数据检查
boxplot(b,las=2,cex.axis=0.6,main='data check')
###转换
ex1<- log2(b+1)
boxplot(ex1,las=2,cex.axis=0.6,main='data check')
二，提取数据
##提取数据。因为之前讲过，我们只看肿瘤周围的肝细胞跟正常肝周细胞对比。所以在此我们提取我们的目的数据。##
PHex <- ex[,1:6]
PHex[1:4,1:6]##查看一下数据。（rawcounts数据）
## PH_WT_tumors_r1 PH_WT_tumors_r2
ENSMUSG00000000001_Gnai3 10.836050 11.437232
ENSMUSG00000000028_Cdc45 5.491853 6.475733
ENSMUSG00000000031_H19 4.321928 8.968667
ENSMUSG00000000037_Scml2 2.584963 3.459432
PH_WT_tumors_r3 PH_WT_notreat_r1
ENSMUSG00000000001_Gnai3 11.070121 11.672867
ENSMUSG00000000028_Cdc45 6.189825 4.857981
ENSMUSG00000000031_H19 5.321928 1.000000
ENSMUSG00000000037_Scml2 3.584963 2.000000
PH_WT_notreat_r2 PH_WT_notreat_r3
ENSMUSG00000000001_Gnai3 11.712957 11.007027
ENSMUSG00000000028_Cdc45 6.321928 3.584963
ENSMUSG00000000031_H19 0.000000 2.000000
ENSMUSG00000000037_Scml2 3.169925 0.000000
PHex1<- ex1[,1:6]
PHex1[1:4,1:6]##查看一下数据。（RPM数据）
## PH_WT_tumors_r1 PH_WT_tumors_r2
ENSMUSG00000000001_Gnai3 6.8132664 7.2351263
ENSMUSG00000000028_Cdc45 1.8813944 2.5238188
ENSMUSG00000000031_H19 1.1104527 4.8070691
ENSMUSG00000000037_Scml2 0.3840887 0.6228264
PH_WT_tumors_r3 PH_WT_notreat_r1
ENSMUSG00000000001_Gnai3 6.9229320 7.52881962
ENSMUSG00000000028_Cdc45 2.3311102 1.36491739
ENSMUSG00000000031_H19 1.6708217 0.07898138
ENSMUSG00000000037_Scml2 0.6924168 0.22504780
PH_WT_notreat_r2 PH_WT_notreat_r3
ENSMUSG00000000001_Gnai3 7.6446714 6.9971839
ENSMUSG00000000028_Cdc45 2.5076949 0.7465790
ENSMUSG00000000031_H19 0.0000000 0.2447131
ENSMUSG00000000037_Scml2 0.5603664 0.0000000
三，画PCA图
library(stringr)
###分组
group_list1=str_split(colnames(PHex),'_',simplify = T)[,3]
group_list2=str_split(colnames(Lex),'_',simplify = T)[,3]
> group_list1
[1] "tumors" "tumors" "tumors" "notreat" "notreat" "notreat"
####PCA看分组情况
library("FactoMineR")
library("factoextra")
####a data frame with n rows (individuals)
####and p columns (numeric variables)
df.pca <- PCA(t(PHex), graph = FALSE)
fviz_pca_ind(df.pca,
geom.ind = "point",
col.ind = group_list1,
addEllipses = TRUE,
legend.title = "Groups"
)
df.pca1 <- PCA(t(PHex1), graph = FALSE)
fviz_pca_ind(df.pca1,
geom.ind = "point",
col.ind = group_list2,
addEllipses = TRUE,
legend.title = "Groups"
)
我觉得从PCA分群来看，这个两个数据格式，应该没有太多的区别，都是很好的区分这个两组。
四，画热图
##画热图。
PHex<-na.omit(PHex)##去一下个别缺失值。
table(is.na.data.frame(PHex))##检测一下缺失值是否去干净，因为热图十不可以有缺失值的。
cg=names(tail(sort(apply(PHex,1,sd)),200))##选两百个去做热图。
PHex <- PHex[cg,]
PHex=t(scale(t(PHex)))
#####查看scale处理后数据的范围
fivenum(PHex)
####目的是避免出现极大极小值影响可视化的效果
###2,-2
PHex[PHex>1.2]=1.2
PHex[PHex< -1.2] = -1.2
library(pheatmap)
pheatmap(PHex)##这个画的比较丑，下面调整一下颜色，加入分组之后会漂亮许多。
####调整下颜色，使正负值颜色的对比会更加鲜明
require(RColorBrewer)
bk = c(seq(-1.2,1.2, length=100))
annotation_col = data.frame(Group = rep(c("tumors", "notreat"), c(3,3)))
rownames(annotation_col)<-colnames(PHex)
####annotation_col和annotation_row的格式需为数据框
####breaks参数用于值匹配颜色
####看下，color和breaks的对应，进行理解
pheatmap(PHex,
breaks=bk,
show_rownames = F,
annotation_col = annotation_col,
color = colorRampPalette(c("navy", "white", "firebrick3"))(100))
####可以调整的内容有是否聚类、是否分组、是否显示行和列的内容等
save(PHex,PHex1,group_list1,group_list2,file = 'ex.Rdata')
##同上换成PHex1再画个图##
热图聚类都可以聚的很好，两个数据很相似。但是具体里面差异基因是不是一致的呢？我觉得通过这两个图只能看一个大概，啥也说明不了。好吧接下来试着复现一下文章中的go功能注释的图把。
五，差异分析（DESeq2+rawcounts）
load('ex.Rdata')
exprSet<- read.table('GSE103788_raw_counts_genes.tsv.gz',header = T,row.names = 1)
exprSet=exprSet[apply(exprSet,1, function(x) sum(x>1) > 6),]
exprSet <- exprSet[,1:6]
group_list <- factor(group_list1)
### Firstly run DEseq2
###
### ---------------
class(exprSet)
suppressMessages(library(DESeq2))
(colData <- data.frame(row.names=colnames(exprSet),
group_list=group_list) )
dds <- DESeqDataSetFromMatrix(countData = exprSet,
colData = colData,
design = ~ group_list)
dds <- DESeq(dds)
res <- results(dds,
contrast=c("group_list",'tumors','notreat'))
resOrdered <- res[order(res$padj),]
head(resOrdered)
DEG =as.data.frame(resOrdered)
DEG = na.omit(DEG)
head(DEG)
##### baseMean log2FoldChange lfcSE
ENSMUSG00000026678_Rgs5 975.6113 3.152064 0.2122172
ENSMUSG00000026473_Glul 7295.5197 3.265248 0.2312006
ENSMUSG00000031375_Bgn 2390.6760 3.661976 0.2659914
ENSMUSG00000021268_Meg3 306.8778 6.370573 0.4675770
ENSMUSG00000002900_Lamb1 242.4983 3.679213 0.2704702
ENSMUSG00000069662_Marcks 526.7389 2.888859 0.2137369
stat pvalue padj
ENSMUSG00000026678_Rgs5 14.85301 6.651063e-50 1.003313e-45
ENSMUSG00000026473_Glul 14.12301 2.740446e-45 2.066981e-41
ENSMUSG00000031375_Bgn 13.76727 4.010678e-43 2.016702e-39
ENSMUSG00000021268_Meg3 13.62465 2.857755e-42 1.077731e-38
ENSMUSG00000002900_Lamb1 13.60303 3.841994e-42 1.159130e-38
ENSMUSG00000069662_Marcks 13.51596 1.259102e-41 3.165593e-38
五，GO注释（这里主要用一下Y叔的优秀的通路分析包，图又漂亮，有方便使用）。
geneList准备，Y叔的注释包（clusterProfiler）只要这个做好后，后面代码基本上不用改了。直接运行了。所以这个很重要。
rm(list = ls())
options(stringsAsFactors = F)
library(ggplot2)
library(clusterProfiler)
library(stringr)
library(org.Mm.eg.db)
keytypes(org.Mm.eg.db)
load('DEG.Rdata')
b <- rownames(DEG)
b=str_split(rownames(DEG),'_',simplify = T)[,1]
rownames(DEG) <- b
gene<- bitr(b, fromType = "ENSEMBL", #fromType是指你的数据ID类型是属于哪一类的
toType = c('ENTREZID'), #toType是指你要转换成哪种ID类型，可以写多种，也可以只写一种
OrgDb = org.Mm.eg.db)#Orgdb是指对应的注释包是哪个
DEG$ENSEMBL<- rownames(DEG)
DEG<- merge(gene,DEG)
## assume that 1st column is ID
## 2nd column is fold change
## feature 1: numeric vector(输入差异倍数列)
geneList <- DEG[,4]
## feature 2: named vector
names(geneList) <- as.character(DEG[,2])##(输入ID列)
## feature 3: decreasing order
geneList <- sort(geneList, decreasing = TRUE)
save(geneList,file='geneList.Rdata')
head(geneList)
>head(geneList)
216643 53973 57430 117586 12323 337924
13.402835 11.655016 11.097527 10.804276 9.923125 9.889644
go注释
rm(list = ls())
options(stringsAsFactors = F)
library(clusterProfiler)
library(org.Mm.eg.db)
load('geneList.Rdata')
##GO classification
gene <- names(geneList)[abs(geneList) > 2]#筛选差异基因大于2的列表
gene.df <- bitr(gene, fromType = "ENTREZID",
toType = c("ENSEMBL", "SYMBOL"),##得到ENSEMBLID与基因名
OrgDb = org.Mm.eg.db)
head(gene.df)
##BP##
ego2 <- enrichGO(gene = gene.df$ENSEMBL,
OrgDb = org.Mm.eg.db,
keyType = 'ENSEMBL',
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05)
ego2 <- setReadable(ego2, OrgDb = org.Mm.eg.db)
##去掉一些重复的通路。
lineage1_ego <- simplify(ego2,
cutoff=0.5,
by="p.adjust",
select_fun=min)
##可视化
barplot(lineage1_ego,showCategory=25)
BP注释总体上跟文章中的通路还是类似的。毕竟如果不知道代码，很难一模一样复现文章的图的。只要趋势差不多就差不多了。
本文参与 腾讯云自媒体同步曝光计划，分享自微信公众号。
原始发表：2020-01-04，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读
目录
腾讯云开发者
