接着day49
一、看看整体表达情况
原始矩阵为变量exprSet,是从txt文件读入得到的数据框
归一化后矩阵为变量exprSet_new,用DESeq2进行归一化后计算而来,是数值型矩阵。从下面class()命令可以看出两个变量属性不同。
class(exprSet)
[1] "data.frame"
class(exprSet_new)
[1] "matrix" "array"
exprSet<-as.matrix(exprSet) #把data.frame格式转唯matix格式
hist(exprSet)
hist(exprSet_new)
说明
hist() 用来做直方图(柱形图),如果是一列数据,就会画成从小到大那样的柱子。是一个矩阵,就把全部数据按照大小和频率列成柱形图。可以看出归一化后的数据看着更均匀。


看一下这两个矩阵的样子:(这两个矩阵都是把变量用write.csv给导出的。)


二、看看不同样本之间表达情况:
par(cex = 0.7) #指定符号的大小
par(mar=c(10, 5, 5, 5))
n.sample=ncol(exprSet) #样本数
if(n.sample>40) par(cex = 0.5)
n.sample=ncol(exprSet)#样本数
cols <- rainbow(n.sample*1.2)
boxplot(exprSet, col = cols,main="expression value",las=2)
boxplot(exprSet_new, col = cols,main="expression value",las=2)
说明
par函数cex参数是用来控制文字和点的大小。cex参数大了表明散点图里的点和各种文字部分的字体都大,比如上面命令中设定为0.7,在样本很多的时候要缩小些字体,所以有了>40,改为0.5这个命令行。
cex对所有的文字进行统一设置外,针对不同的标题,还有对应的cex系列参数的用法很多,比如调节主标题,副标题字体大小,调节坐标轴的字体大小等,可以看下面链接。
参考教程:https://www.bbsmax.com/A/A7zg1Erl54/
而par函数里的mar参数则是设定图形四周的留空大小。
参考教程:https://www.yht7.com/news/136868


两个图都是箱型图,可归一化以后的更好看些,在有些测序结果里看到过,主要是看各个样本表达水平有没有大的不同。如果某个样本比别人差太多,那就是有问题,还有可能是这个细胞有什么特殊处理。
原始数据那个图太丑了。。
三、查看单一基因在组间差异
plotCounts(dds2, "CD38", intgroup = "group_list") #把前面dds2中的nomalization之后的某个特定基因表达水平展示出来。
说明
plotCounts也是DESeq2包中的一个命令,dds是DESeqDataSet., gene是一个特殊基因名称,intgroup:在colData(x)中,用于进行分组的名称。

如果有很多想看的基因,可以批量执行。
sigGene= read.csv('treat_vs_con_sig.csv',row.names=1) #读取文件
for(i in rownames(sigGene)){
genename=paste(i,sep="") #基因名
pdf(file=paste(i,'_counts.pdf',sep="")) #以pdf格式输出
plotCounts(dds2, genename, intgroup = "group_list")
dev.off() #不在Rstudio里输出,直接保存在当前工作路径下
}
很快就在工作目录下生成了这么多文件。

挑一个打开看看。这个数值是归一化后的counts值,找来表格里的数对照一下,还真没错。如果样本多一些,画的散点图应该挺好看的。


上面这个表里的数是log2之后的,所以原来的数据要用2(x)来计算一下。比如那个9.479736的,转成2(9.479736)=713.9781。
三、PCA聚类
tiff(filename = "PCA.tiff", width = 1500, height = 1500, units = "px", res = 150)
plotPCA(rld, intgroup="group_list")
dev.off() #不在Rstudio里输出,直接保存在当前工作路径下这个PCA.tiff文件了

网友评论