拖後腿學徒居然也完成作業，理解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注釋總體上跟文章中的通路還是類似的。畢竟如果不知道代碼，很難一模一樣複現文章的圖的。隻要趨勢差不多就差不多了。