拖后腿学徒居然也完成作业，理解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注释总体上跟文章中的通路还是类似的。毕竟如果不知道代码，很难一模一样复现文章的图的。只要趋势差不多就差不多了。