根据表达量筛选探针后，对主成分分析的PCA图有什么影响

2019-11-21 本文已影响0人小梦游仙境

这个数据集是来自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处理后均有意义的基因。