这个数据集是来自
GSE106292.关于PCA之前都是都是不是零或者是零但是肉眼不可见的矩阵，但是这次的测序数据一眼望去都是0，所以懵圈。但是老大给代码的PCA就是原始的矩阵+log处理做的，并没有做任何去除0的处理，所以我就想着那么不处理肯定没问题，那处理了会对PCA图有什么区别，然后分成了==exprSet1==、==exprSet2==、==exprSet3==
rm(list = ls())
load('step1-getdata.Rdata')
tmp<-pd$sample
rownames(pd)<-1:35
pd1<-pd[c(8:10,18:20,11:13,33:35,27:29),]
pd1$tissue<-rep(c('ligamentocyte','chondrocyte','tenocyte','osteoblast','myoblast'),each=3)
expr1<-expr[,match((pd1$sample),colnames(expr))]
save(expr1,pd1,file = 'expr1_pd1')
exprSet<-expr1
table(apply(exprSet,1, function(x) sum(x>1)))
table(apply(exprSet,1, function(x) sum(x>1)==0))
table(apply(exprSet,1, function(x) sum(x>1) > 5))
#下面是3个exprSet：exprSet1、exprSet2、exprSet3，后面分别用这三个做PCA
exprSet1=apply(exprSet, 2,as.numeric)
rownames(exprSet1)=rownames(exprSet)
exprSet2<-exprSet[apply(exprSet1,1,sum)!=0,]
exprSet3=exprSet[apply(exprSet2,1,function(x) sum(x>1) > 5),]
# exprSet3<-as.data.frame(exprSet2)
# exprSet3=2^exprSet3-1
# dat=log2(edgeR::cpm(exprSet3)+1)
group_list<-rep(c('ligamentocyte','chondrocyte','tenocyte','osteoblast','myoblast'),each=3)
# 过滤有表达的
save(exprSet3,group_list,file = 'step2-checkdata.Rdata')
###PCA
rm(list = ls())
load('step2-checkdata.Rdata')
dat<-exprSet1
dat<-log2(dat+1)
dat[1:4,1:4]
dat<-log2(dat+1)
## 下面是画PCA的必须操作，需要看说明书。
dat=t(dat)#画PCA图时要求是行名时样本名，列名时探针名，因此此时需要转换
dat=as.data.frame(dat)#将matrix转换为data.frame
dat=cbind(dat,group_list) #cbind横向追加，即将分组信息追加到最后一列
#dat<-as.data.frame(dat)
library("FactoMineR")#画主成分分析图需要加载这两个包
library("factoextra")
# The variable group_list (index = 54676) is removed
# before PCA analysis
dat.pca <- PCA(dat[,-ncol(dat)], graph = FALSE)#现在dat最后一列是group_list，需要重新赋值给一个dat.pca,这个矩阵是不含有分组信息的
fviz_pca_ind(dat.pca,
geom.ind = "point", # show points only (nbut not "text")
col.ind = dat$group_list, # color by groups
palette = c("#7ec0ee", "#ffaeb9", "#7d26d0", "#ff00ff", "#4876ff"),
addEllipses = TRUE, # Concentration ellipses
legend.title = "Groups"
)
1.PCA图
- 下面👇的dat是上面的
exprSet1
image-20191121210426802
- 下面👇的dat是上面的
exprSet3，是进行了sum(x>1) > 5)，可以看到去掉了一些基因（大概2000个左右）后，有两个样本离的特别近
image-20191121210450784
- 下面👇的dat是上面的
exprSet2，仅仅是做了exprSet2<-exprSet[apply(exprSet1,1,sum)!=0,]的处理，可以看到和上面的exprSet1得到的结果是一样的
image-20191121210459639
PCA得出的初步结论
- 思考如下
-
不能随便去除表达量为0的探针，因为这些基因是有意义的，有用，在主成分分析时，明明可以区分开的，但是用过滤条件apply(exprSet,1, function(x) sum(x>1) > 5))去掉了一部分以后，就会使本来能区分开的样本距离变近，说明这些基因是在作者进行了TPM处理后均有意义的基因。
2.热图
- 未进行log时的热图
scale前
image-20191121210511794
scale后
image-20191121210518787
加上分组
image-20191121210540581
-
进行==log了==的热图
scale前
image-20191121210549654
scale后
image-20191121210557557
加上分组
image-20191121210607390
原文中,和==原文比较==，感觉不是完全的一样，不过也没关系
image-20191121210616899
从上面的两张log和未log的矩阵热图比较，热图是一样的
heatmap得出的初步结论
- 由于可能是下一步有scale这个归一化的处理，所以前面即使有：有 上万也有0的矩阵，就是未log的矩阵，得到的热图展示是一样的。
网友评论