R语言绘图-常用参数

先介绍一下原始数据：
该数据最开始是一套从NCBI下载的基因芯片数据，数据编号为GSE29272
发表该套数据的文章名字为：Affymetrix gene expression array data for cardia and non-cardia gastric cancer samples
该初始数据的下载网址为：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE29272
这套数据的基本信息都可以在上面的网址中查看到，小编在下载之后又对这一套数据进行了一定的预处理，然后又对处理过的数据进行了差异表达基因的筛选。从中筛选到的具有明显表达水平差异性的数据作为本次绘图的原始数据。这些数据来自于共168个样本的13个基因探针，在文章的最后会提供给需要的你。

再介绍一下背景：
在做完差异表达基因的筛选之后，小编对样本进行了谱系聚类，但是聚类结果出了一点小问题，

本来应该聚类成为两个大类（Normal、Tumor）的样本，在最右边又单独出现了一个小类。虽然这个小类只有六个样本的大小，而且也是和Tumor分在一支上，但是！但是！身为一个完美主义者，这个根本不能忍啊！我就开始了使用R语言的探索过程……

再介绍一下背景：
在做完差异表达基因的筛选之后，小编对样本进行了谱系聚类，但是聚类结果出了一点小问题，本来应该聚类成为两个大类（Normal、Tumor）的样本，在最右边又单独出现了一个小类。虽然这个小类只有六个样本的大小，而且也是和Tumor分在一支上，但是！但是！身为一个完美主义者，这个根本不能忍啊！我就开始了使用R语言的探索过程…

那一坨坨的是什么东西，那只是样本名重叠在一起了而已，虽然我室友都吐槽说像黑叔叔们的卷发 …

一、代码的读取和简单处理

首先读取原始数据，然后进行简单的数据提取，代码如下：

data<-read.csv(“PATH/差异表达基因t检验60.csv”,header = T)
test<-data[,2:169]

简单解释一下，read.csv()是R语言中读取CSV格式文件的一种方法，后面参数header指的是读入的数据是否带有表头。
我提供数据中只有第2列到第169列是我们需要的，第一列和最后两列并不是我们需要的数据，所以我将其去除。
小贴士：可使用nrow()和ncol()函数，来查看数据的总行数和总列数。当然如果你使用了Rstudio的话在右侧的数据栏中你可以轻易的查看行列数。

二、对数据进行主成分分析以及k-均值聚类
这一步骤不在我们绘图介绍的内容之中，主成分分析只是对数据的一种处理，所以不在此处过多赘述，你只需要知道我们在这一步中获得了一个新的数据pca_data。

而k-均值聚类，也是一种聚类的方法，是对先前的数据pca_data进行聚类分析使用的，它可以生成一个与pca_data中样本一一对应的分类结果。详情咨询百度。代码如下：

pca<-princomp(t(test),cor=T)
summary(pca)
pca_data <- predict(pca)
library(stats)
fit_km1 =kmeans(pca_data,center=2)

三、绘图正式开始
先拿出我们最先讲到的函数plot()，对pca_data进行绘图，代码如下：

plot(pca_data)

(一) kmeans聚类

#读取数据
data<-read.csv("F:/R/test/差异表达基因t检验60.csv",header = T)
#增加基因名称
test<-data[,2:169]

pca<-princomp(t(test),cor=T)
summary(pca)
pca_data <- predict(pca)
library(stats)
fit_km1 =kmeans(pca_data,center=2)

plot(pca_data)

运行结果如下：

这显然不是我们想要的直观分类结果。

我们在之前的函数中继续添加参数col，这个参数就是控制颜色的参数(color)，对于这个函数的赋值，你可以直接赋值为数字(1、2、3、4、5、6…)，也可以使用“red”、“green”、“blue” 等来赋值，但是注意这样做的时候，对应的颜色要用双引号括起来。你也可以使用一组对应的颜色向量来对其赋值，举例：col=1、col=“red”、col=1:3、col=c(“red”.“green”,“blue”)；修改代码如下：

plot(pca_data,col=(fit_km1$cluster)*2)

(二) fit_km1$cluster是我们聚类的结果

#读取数据
data<-read.csv("F:/R/test/差异表达基因t检验60.csv",header = T)
#增加基因名称
test<-data[,2:169]

pca<-princomp(t(test),cor=T)
summary(pca)

pca_data <- predict(pca)
library(stats)
fit_km1 =kmeans(pca_data,center=2)

plot(pca_data,col=(fit_km1$cluster)*2)

运行结果如下：

简单解释一下，这里赋值的fit_km1$cluster是我们聚类的结果，他的本质是一组数字向量，至于乘2，是因为默认1为黑色，黑色并不是很适合图像的展示，所以用简单的乘2来改变它的颜色。

根据赋予了色彩的图像基本上就可直接看的出来他被明显的分为了红色和蓝色两大类。

但是对于展示来说，我们不仅要看到聚类的结果，也要看到什么样的样本被聚在了一起，我们在尝试引入一个参数pch，这个参数是用来修改图中图形元素(plotting character)的，接受的赋值为数字或者数字向量，举例pch=1、pch=c(1,2)。因为我的样本本身就是Normal和Tumor交替出现的，所以修改代码:

plot(pca_data,col=(fit_km1$cluster)*2,pch=c(1,2))

(三)参数pch

#读取数据
data<-read.csv("F:/R/test/差异表达基因t检验60.csv",header = T)
#增加基因名称
test<-data[,2:169]

pca<-princomp(t(test),cor=T)
summary(pca)

pca_data <- predict(pca)
library(stats)
fit_km1 =kmeans(pca_data,center=2)

plot(pca_data,col=(fit_km1$cluster)*2,pch=c(1,2))

可以看出圆圈基本上被分到了红色聚集的地方，而三角则都聚集在另一边。

此时肯定会有人说不喜欢圆圈和三角，那好吧，我只能一抬手——甩给你25个其他选择，总有一款适合你：

此外，我们再次进行修改，在推出两个参数lwd和cex分别是线条宽度和图像元素的大小，只接受数字赋值，例如：lwd=2,cex=2；这些都是指默认参数的两倍。再次修改代码：

plot(pca_data,col=(fit_km1$cluster)*2,lwd=2,cex=1.5,pch=c(1,2))

(四)两个参数lwd和cex

#读取数据
data<-read.csv("F:/R/test/差异表达基因t检验60.csv",header = T)
#增加基因名称
test<-data[,2:169]

pca<-princomp(t(test),cor=T)
summary(pca)

pca_data <- predict(pca)
library(stats)
fit_km1 =kmeans(pca_data,center=2)

plot(pca_data,col=(fit_km1$cluster)*2,lwd=2,cex=1.5,pch=c(1,2))

运行结果如下：

看得出来，效果很明显。哦，对了，我最开始的目的是要看看那六个奇怪的样本在哪，那就再使用一个函数points()，这个函数是用来在已经绘制出来的图像上添加新的元素点的。使用方法和plot()几乎相同，我们尝试找出这几个样本，代码如下：

plot(pca_data,col=(fit_km1$cluster)*2,lwd=2,cex=1.5,pch=c(1,2))
points(pca_data[sp,],col=”blue”,lwd=2,pch=17)

(五)函数points()

#读取数据
data<-read.csv("F:/R/test/差异表达基因t检验60.csv",header = T)
#增加基因名称
test<-data[,2:169]

pca<-princomp(t(test),cor=T)
summary(pca)

pca_data <- predict(pca)
library(stats)
fit_km1 =kmeans(pca_data,center=2)

plot(pca_data,col=(fit_km1$cluster)*2,lwd=2,cex=1.5,pch=c(1,2))


points(pca_data[sp,],col="blue",lwd=2,pch=17)

结果如下：

嗯，果然这六个样本点离蓝色区域较远，同时又有靠近红色的趋势，所以这又代表了啥？
管他呢 ，今天就写到这了。

原始数据及所用代码下载：

原始数据及代码

#读取数据
data<-read.csv("F:/R/test/差异表达基因t检验60.csv",header = T)
#增加基因名称
test<-data[,2:169]

#谱系聚类
d<-dist(t(test), method = "euclidean")
hc<-hclust(d,"single")
plot(hc)
rect.hclust(hc,k=3)
result=cutree(hc,k=3)
result
which(result==1)
which(result==2)
sp<-which(result==3)
#PCA
pca<-princomp(t(test),cor=T)
summary(pca)
pca_data <- predict(pca)

#k-均值聚类
library(stats)
fit_km1 =kmeans(pca_data,center=2)
#绘图
plot(pca_data,col=(fit_km1$cluster)*2,lwd=2,cex=1.5,pch=c(1,2))
points(pca_data[sp,],col="blue",lwd=2,cex=1.5,pch=17)