15高通量测序-R与MDS/PCoA
R与MDS/PCoA
构造数据集
加载ggplot2包，构造数据集，数据将由一个矩阵组成，10列对应10个样本，100行对应100个基因的测量值。前5个样品为“wt”或“野生型”样品,“wt”样本是正常的,最后5个样品为“ko”或"敲除"样品，样本为我们敲掉一个基因
library(ggplot2)
data.matrix <- matrix(nrow=100, ncol=10)
colnames(data.matrix) <- c(
paste("wt", 1:5, sep=""),
paste("ko", 1:5, sep=""))
rownames(data.matrix) <- paste("gene", 1:100, sep="")
for (i in 1:100){
wt.values <- rpois(5, lambda = sample(x=10:1000, size = 1))
ko.values <- rpois(5, lambda = sample(x=10:1000, size = 1))
data.matrix[i,] <- c(wt.values, ko.values)
}
> head(data.matrix)
wt1 wt2 wt3 wt4 wt5 ko1 ko2 ko3 ko4 ko5
gene1 676 721 660 648 674 385 348 356 370 331
gene2 376 353 320 364 347 158 177 165 163 179
gene3 966 909 981 938 965 949 948 936 1016 1031
gene4 239 243 248 255 256 844 897 919 899 922
gene5 141 127 152 145 125 176 172 175 179 186
gene6 129 115 105 128 129 77 63 63 88 50
PCA分析
现在，为了比较，我们将对数据集进行PCA。我将轻松地通过这些步骤，因为我们已经在R StatQuest的PCA中讨论过它们。
pca <- prcomp(t(data.matrix),scale=TRUE)
pca$x
plot(pca$x[,1],pca$x[,2])
pca.var <- pca$sdev^2
pca.var.per <- round(pca.var/sum(pca.var)*100,1)
barplot(pca.var.per,
main="Scree Plot",
xlab="Principal Component",
ylab="Percent Variation")
pca.data <- data.frame(Sample=rownames(pca$x),
X=pca$x[,1],
Y=pca$x[,2])
pca.data
ggplot(data=pca.data, aes(x=X, y=Y, label=Sample))+
geom_text() +
xlab(paste("PC1 - ", pca.var.per[1], "%", sep=""))+
ylab(paste("PC2 - ", pca.var.per[2], "%", sep=""))+
theme_bw()+
ggtitle("My PCA Graph")
Rplot03.png
最终得到上图，其中PC1占据了变异的88.9%，PC2占据了变异的2.9%，WT和KO样本之间有巨大差异。
MDS/PCoA分析
步骤1:创建距离矩阵。我们用dist())函数来实现这一点。就像PCA一样，我们对矩阵进行转置使得样本是行。我们还对每个基因的测量值进行居中和缩放(现在是列)。最后，我们告诉dist()函数，我们希望它使用欧几里得距离创建矩阵。注意:dist()有6个不同的距离指标可供选择!
distance.matrix <- dist(scale(t(data.matrix), center=TRUE, scale=TRUE),
method="euclidean")
步骤2:使用cmdscale()对距离矩阵进行多维度缩放。Classical Muti-Dimensional Scaling=cmdscale。我们告诉cmdscale()我们希望eig=TRUE返回特征值,用这些特征来计算最终MDS图中的每个坐标占距离矩阵总变异的比例。x.ret=TRUE用于表示是否返回双中心对称的距离矩阵。如果使用eigen()函数替代cmdscale()函数来演示做MDS分析时，这就比较有用。
mds.stuff <- cmdscale(distance.matrix, eig=TRUE, x.ret=TRUE)
步骤3::计算特征值在MDS图中各轴所占的变化量
mds.var.per <- round(mds.stuff$eig/sum(mds.stuff$eig)*100,1)
mds.var.per
步骤4:用ggplot2绘图
mds.values <- mds.stuff$points
mds.data <- data.frame(Sample=rownames(mds.values),
X=mds.values[,1],
Y=mds.values[,2])
ggplot(data=mds.data, aes(x=X, y=Y, label=Sample))+
geom_text() +
xlab(paste("MDS1 - ", mds.var.per[1], "%", sep=""))+
ylab(paste("MDS2 - ", mds.var.per[2], "%", sep=""))+
theme_bw()+
ggtitle("MDS plot using Euclidean distance")
Rplot04.png
就像在PCA图中一样，野生型样本在图的左侧,而敲除样本就在右边。就像PCA图一样，MDS1占数据变化的88.9%，MDS2占据了变异的2.9%,实际上，PCA图和MDS图不仅仅看起来相似，它们是完全相同的!! 这是因为我们使用了欧几里得度量(euclidean metric)来计算距离矩阵。
log fold change的绝对值的平均值
现在让我们看看当我们使用其他距离计算距离矩阵时会发生什么。我们用log fold change的绝对值的平均值。对于基因表达的差异分析，通常使用edgeR来成，现在我们使用plotMDS()函数进行表达差异的log值计算。
计算这个矩阵的log2的值
log2.data.matrix <- log2(data.matrix)
因为log fold change的绝对值的平均值并不是dist()中内置的距离指标之一，所以我们将手动创建自己的距离矩阵。在这一步中，我们只是创建了一个空矩阵。
log2.data.matrix <- log2(data.matrix)
log2.distance.matrix<- matrix(0,
nrow=ncol(log2.data.matrix),
ncol=ncol(log2.data.matrix),
dimnames = list(colnames(log2.data.matrix),
colnames(log2.data.matrix)))
log2.distance.matrix
> log2.distance.matrix
wt1 wt2 wt3 wt4 wt5 ko1 ko2 ko3 ko4 ko5
wt1 0 0 0 0 0 0 0 0 0 0
wt2 0 0 0 0 0 0 0 0 0 0
wt3 0 0 0 0 0 0 0 0 0 0
wt4 0 0 0 0 0 0 0 0 0 0
wt5 0 0 0 0 0 0 0 0 0 0
ko1 0 0 0 0 0 0 0 0 0 0
ko2 0 0 0 0 0 0 0 0 0 0
ko3 0 0 0 0 0 0 0 0 0 0
ko4 0 0 0 0 0 0 0 0 0 0
ko5 0 0 0 0 0 0 0 0 0 0
向这个矩阵填log fold change的绝对值的平均值
for(i in 1:ncol(log2.distance.matrix)){ for(j in 1:i){
log2.distance.matrix[i,j]<-
mean(abs(log2.data.matrix[,i] - log2.data.matrix[,j]))
}
}
> log2.distance.matrix
wt1 wt2 wt3 wt4 wt5 ko1 ko2
wt1 0.00000000 0.00000000 0.0000000 0.00000000 0.000000 0.00000000 0.00000000
wt2 0.08670204 0.00000000 0.0000000 0.00000000 0.000000 0.00000000 0.00000000
wt3 0.09634131 0.09062371 0.0000000 0.00000000 0.000000 0.00000000 0.00000000
wt4 0.09939116 0.09688245 0.1066254 0.00000000 0.000000 0.00000000 0.00000000
wt5 0.09159530 0.09475646 0.1040647 0.08427843 0.000000 0.00000000 0.00000000
ko1 1.28509779 1.26807514 1.2659031 1.27554608 1.277869 0.00000000 0.00000000
ko2 1.26764742 1.25118819 1.2481033 1.25743510 1.260040 0.08848651 0.00000000
ko3 1.28493452 1.26855651 1.2628855 1.27616884 1.278074 0.07867476 0.08272827
ko4 1.28564850 1.27050213 1.2650763 1.27768246 1.277032 0.09501419 0.09984947
ko5 1.28946202 1.27372474 1.2706313 1.27994170 1.282890 0.08644142 0.09483029
ko3 ko4 ko5
wt1 0.00000000 0.0000000 0
wt2 0.00000000 0.0000000 0
wt3 0.00000000 0.0000000 0
wt4 0.00000000 0.0000000 0
wt5 0.00000000 0.0000000 0
ko1 0.00000000 0.0000000 0
ko2 0.00000000 0.0000000 0
ko3 0.00000000 0.0000000 0
ko4 0.08426407 0.0000000 0
ko5 0.09299244 0.1132473 0
现在我们对新的距离矩阵进行多维度缩放。as.dist(log2.distance.matrix)我只是把我们自制的矩阵转换成一个“真正的”距离矩阵，所以cmdscale())知道它是如何工作的。
mds.stuff <- cmdscale(as.dist(log2.distance.matrix),
eig=TRUE,
x.ret=TRUE)
接着计算MDS图中每个坐标所占总变异的比例
mds.var.per <- round(mds.stuff$eig/sum(mds.stuff$eig)*100,1)
mds.var.per
##转换成ggplot2能识别的数据格式
mds.values <- mds.stuff$points
mds.data <- data.frame(Sample=rownames(mds.values),
X=mds.values[,1],
Y=mds.values[,2])
mds.data
ggplot(data=mds.data, aes(x=X, y=Y, label=Sample)) +
geom_text()+
theme_bw()+
xlab(paste("MDS1 - ", mds.var.per[1],"%", sep=""))+
ylab(paste("MDS2 - ", mds.var.per[2],"%", sep=""))+
ggtitle("MDS plot using avg(logFC) as the distance")
Rplot05.png
两种不同的MDS图(一种使用欧几里得距离，另一种使用log fold change的绝对值的平均值是相似的，但不相同。在新的图表中，x轴占了更多的变化(99.1% vs 88.9%)。
