拖后腿学徒居然也完成作业，理解RNA-seq数据分析结果

前面我出了一个学徒作业，下载表达矩阵后绘制PCA图及热图，然后理解作者给出的RPM和raw_counts的差异，详见：很意外，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_r1PH_WT_tumors_r2
ENSMUSG00000000001_Gnai318272772
ENSMUSG00000000028_Cdc454488
ENSMUSG00000000031_H1919500
ENSMUSG00000000037_Scml2510
PH_WT_tumors_r3PH_WT_notreat_r1
ENSMUSG00000000001_Gnai321493264
ENSMUSG00000000028_Cdc457228
ENSMUSG00000000031_H19391
ENSMUSG00000000037_Scml2113

##RPM##
b<-read.table('GSE103788_filtered_RPM_genes.tsv.gz',header=T,row.names=1)
b[1:4,1:4]
##PH_WT_tumors_r1PH_WT_tumors_r2
ENSMUSG00000000001_Gnai3111.4598680149.6572460
ENSMUSG00000000028_Cdc452.68430994.7510237
ENSMUSG00000000031_H191.159133826.9944527
ENSMUSG00000000037_Scml20.30503520.5398891
PH_WT_tumors_r3PH_WT_notreat_r1
ENSMUSG00000000001_Gnai3120.3417347183.67178123
ENSMUSG00000000028_Cdc454.03192411.57561577
ENSMUSG00000000031_H192.18395890.05627199
ENSMUSG00000000037_Scml20.61598840.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='datacheck')
###去除低碱基质量的基因。
a=a[apply(a,1,function(x)sum(x>1)>6),]

ex<-log2(a+1)##转换
boxplot(ex,las=2,cex.axis=0.6,main='datacheck')

下图是转换前后的图片。

同上方法，处理一下RPM的数据。
##########boxplot可视化数据检查
boxplot(b,las=2,cex.axis=0.6,main='datacheck')
###转换
ex1<-log2(b+1)
boxplot(ex1,las=2,cex.axis=0.6,main='datacheck')

二，提取数据 ##提取数据。因为之前讲过，我们只看肿瘤周围的肝细胞跟正常肝周细胞对比。所以在此我们提取我们的目的数据。##
PHex<-ex[,1:6]
PHex[1:4,1:6]##查看一下数据。（rawcounts数据）
##PH_WT_tumors_r1PH_WT_tumors_r2
ENSMUSG00000000001_Gnai310.83605011.437232
ENSMUSG00000000028_Cdc455.4918536.475733
ENSMUSG00000000031_H194.3219288.968667
ENSMUSG00000000037_Scml22.5849633.459432
PH_WT_tumors_r3PH_WT_notreat_r1
ENSMUSG00000000001_Gnai311.07012111.672867
ENSMUSG00000000028_Cdc456.1898254.857981
ENSMUSG00000000031_H195.3219281.000000
ENSMUSG00000000037_Scml23.5849632.000000
PH_WT_notreat_r2PH_WT_notreat_r3
ENSMUSG00000000001_Gnai311.71295711.007027
ENSMUSG00000000028_Cdc456.3219283.584963
ENSMUSG00000000031_H190.0000002.000000
ENSMUSG00000000037_Scml23.1699250.000000
PHex1<-ex1[,1:6]
PHex1[1:4,1:6]##查看一下数据。（RPM数据）
##PH_WT_tumors_r1PH_WT_tumors_r2
ENSMUSG00000000001_Gnai36.81326647.2351263
ENSMUSG00000000028_Cdc451.88139442.5238188
ENSMUSG00000000031_H191.11045274.8070691
ENSMUSG00000000037_Scml20.38408870.6228264
PH_WT_tumors_r3PH_WT_notreat_r1
ENSMUSG00000000001_Gnai36.92293207.52881962
ENSMUSG00000000028_Cdc452.33111021.36491739
ENSMUSG00000000031_H191.67082170.07898138
ENSMUSG00000000037_Scml20.69241680.22504780
PH_WT_notreat_r2PH_WT_notreat_r3
ENSMUSG00000000001_Gnai37.64467146.9971839
ENSMUSG00000000028_Cdc452.50769490.7465790
ENSMUSG00000000031_H190.00000000.2447131
ENSMUSG00000000037_Scml20.56036640.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")
####adataframewithnrows(individuals)
####andpcolumns(numericvariables)
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)
###FirstlyrunDEseq2
###
###---------------
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)

#####baseMeanlog2FoldChangelfcSE
ENSMUSG00000026678_Rgs5975.61133.1520640.2122172
ENSMUSG00000026473_Glul7295.51973.2652480.2312006
ENSMUSG00000031375_Bgn2390.67603.6619760.2659914
ENSMUSG00000021268_Meg3306.87786.3705730.4675770
ENSMUSG00000002900_Lamb1242.49833.6792130.2704702
ENSMUSG00000069662_Marcks526.73892.8888590.2137369
statpvaluepadj
ENSMUSG00000026678_Rgs514.853016.651063e-501.003313e-45
ENSMUSG00000026473_Glul14.123012.740446e-452.066981e-41
ENSMUSG00000031375_Bgn13.767274.010678e-432.016702e-39
ENSMUSG00000021268_Meg313.624652.857755e-421.077731e-38
ENSMUSG00000002900_Lamb113.603033.841994e-421.159130e-38
ENSMUSG00000069662_Marcks13.515961.259102e-413.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)

##assumethat1stcolumnisID
##2ndcolumnisfoldchange

##feature1:numericvector(输入差异倍数列)
geneList<-DEG[,4]

##feature2:namedvector
names(geneList)<-as.character(DEG[,2])##(输入ID列)

##feature3:decreasingorder
geneList<-sort(geneList,decreasing=TRUE)
save(geneList,file='geneList.Rdata')
head(geneList)
>head(geneList)
216643539735743011758612323337924
13.40283511.65501611.09752710.8042769.9231259.889644

go注释

rm(list=ls())
options(stringsAsFactors=F)
library(clusterProfiler)
library(org.Mm.eg.db)
load('geneList.Rdata')

##GOclassification
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注释总体上跟文章中的通路还是类似的。毕竟如果不知道代码，很难一模一样复现文章的图的。只要趋势差不多就差不多了。