前面我出了一个学徒作业，下载表达矩阵后绘制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_r2ENSMUSG00000000001_Gnai318272772ENSMUSG00000000028_Cdc454488ENSMUSG00000000031_H1919500ENSMUSG00000000037_Scml2510PH_WT_tumors_r3PH_WT_notreat_r1ENSMUSG00000000001_Gnai321493264ENSMUSG00000000028_Cdc457228ENSMUSG00000000031_H19391ENSMUSG00000000037_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_r2ENSMUSG00000000001_Gnai3111.4598680149.6572460ENSMUSG00000000028_Cdc452.68430994.7510237ENSMUSG00000000031_H191.159133826.9944527ENSMUSG00000000037_Scml20.30503520.5398891PH_WT_tumors_r3PH_WT_notreat_r1ENSMUSG00000000001_Gnai3120.3417347183.67178123ENSMUSG00000000028_Cdc454.03192411.57561577ENSMUSG00000000031_H192.18395890.05627199ENSMUSG00000000037_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_r2ENSMUSG00000000001_Gnai310.83605011.437232ENSMUSG00000000028_Cdc455.4918536.475733ENSMUSG00000000031_H194.3219288.968667ENSMUSG00000000037_Scml22.5849633.459432PH_WT_tumors_r3PH_WT_notreat_r1ENSMUSG00000000001_Gnai311.07012111.672867ENSMUSG00000000028_Cdc456.1898254.857981ENSMUSG00000000031_H195.3219281.000000ENSMUSG00000000037_Scml23.5849632.000000PH_WT_notreat_r2PH_WT_notreat_r3ENSMUSG00000000001_Gnai311.71295711.007027ENSMUSG00000000028_Cdc456.3219283.584963ENSMUSG00000000031_H190.0000002.000000ENSMUSG00000000037_Scml23.1699250.000000PHex1<-ex1[,1:6]phex1[1:4,1:6]##查看一下数据。（rpm数据）##ph_wt_tumors_r1ph_wt_tumors_r2ensmusg00000000001_gnai36.81326647.2351263ensmusg00000000028_cdc451.88139442.5238188ensmusg00000000031_h191.11045274.8070691ensmusg00000000037_scml20.38408870.6228264ph_wt_tumors_r3ph_wt_notreat_r1ensmusg00000000001_gnai36.92293207.52881962ensmusg00000000028_cdc452.33111021.36491739ensmusg00000000031_h191.67082170.07898138ensmusg00000000037_scml20.69241680.22504780ph_wt_notreat_r2ph_wt_notreat_r3ensmusg00000000001_gnai37.64467146.9971839ensmusg00000000028_cdc452.50769490.7465790ensmusg00000000031_h190.00000000.2447131ensmusg00000000037_scml20.56036640.0000000 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,-2phex[phex>1.2]=1.2PHex[PHex<-1.2]=-1.2library(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)#####baseMeanlog2FoldChangelfcSEENSMUSG00000026678_Rgs5975.61133.1520640.2122172ENSMUSG00000026473_Glul7295.51973.2652480.2312006ENSMUSG00000031375_Bgn2390.67603.6619760.2659914ENSMUSG00000021268_Meg3306.87786.3705730.4675770ENSMUSG00000002900_Lamb1242.49833.6792130.2704702ENSMUSG00000069662_Marcks526.73892.8888590.2137369statpvaluepadjENSMUSG00000026678_Rgs514.853016.651063e-501.003313e-45ENSMUSG00000026473_Glul14.123012.740446e-452.066981e-41ENSMUSG00000031375_Bgn13.767274.010678e-432.016702e-39ENSMUSG00000021268_Meg313.624652.857755e-421.077731e-38ENSMUSG00000002900_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)<-bgene<-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:namedvectornames(geneList)<-as.character(DEG[,2])##(输入ID列)##feature3:decreasingordergeneList<-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")##GOclassificationgene<-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注释总体上跟文章中的通路还是类似的。毕竟如果不知道代码，很难一模一样复现文章的图的。只要趋势差不多就差不多了。
