生成虛擬數(shù)據(jù)集
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 = "")
head(data.matrix)
以上代碼生成了100行基因，10列樣本的矩陣
前五列命名wt開頭+1-5，表示正?；?br> 后五列命名ko開頭+1-5，表示缺少基因的樣本（knock-out）
給每行基因都統(tǒng)一命名gene+1-100
head()函數(shù)默認(rèn)查看前6行
現(xiàn)在只是定義了矩陣的shape和name，還沒填充數(shù)值
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)
這段代碼的作用是生成一個(gè)大小為100x10的數(shù)據(jù)矩陣data.matrix，其中前5列是"wt"（wild-type）樣本的值，后5列是"ko"（knockout）樣本的值。
在循環(huán)中，對(duì)于每個(gè)i的取值（從1到100），首先使用sample(x=10:1000, size = 1)從10到1000之間的整數(shù)中隨機(jī)抽取一個(gè)數(shù)作為泊松分布的參數(shù)lambda。然后，使用rpois(5, lambda)函數(shù)生成一個(gè)具有泊松分布的隨機(jī)數(shù)向量，其中每個(gè)元素表示一個(gè)基因在"wt"樣本中的表達(dá)量。同樣的過程也用于生成"ko"樣本中的表達(dá)量。
最后，通過c(wt.values, ko.values)將"wt"和"ko"樣本的表達(dá)量合并為一個(gè)長度為10的向量，并將其賦值給data.matrix的第i行。
用for循環(huán)給依次給1-100行的前五列和后五列賦值，填充值介于10-1000之間。
初始虛擬數(shù)據(jù)集創(chuàng)建完畢，接下來用prcomp()函數(shù)分析各樣本之間關(guān)系。該函數(shù)默認(rèn)情況下以基因?yàn)榱校瑯颖緸樾?，和我們?chuàng)建的矩陣互為轉(zhuǎn)置，因此需要用到轉(zhuǎn)置函數(shù)t()
pca <- prcomp(t(data.matrix), scale = TRUE)
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 = "Screen Plot", xlab ="principal component", ylab = "percent variation")
prcomp(, scale = TRUE)表示對(duì)數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化
每一列表示一個(gè)基因所對(duì)應(yīng)的10個(gè)樣本，即一列只有十個(gè)數(shù)據(jù)
plot 生成一個(gè)2D的圖，前兩個(gè)主成分的散點(diǎn)圖
pca.var 表示 標(biāo)準(zhǔn)差的平方
pca.var.per 表示 每個(gè)變量所占的百分比，保留小數(shù)點(diǎn)后一位
可以看到前兩個(gè)成分所占比例最大，尤其是第一個(gè)成分
用 barplot 來直觀每個(gè)成分所占比例
ggplot2繪圖
pca.data <- data.frame(sample = rownames(pca$x),
X = pca$x[,1],
Y = pca$x[,2]
)
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 pac graph")
先按 ggplot2 需要的方式格式化數(shù)據(jù)，x軸用第一個(gè)成分，y軸用第二個(gè)成分
可以發(fā)現(xiàn)數(shù)據(jù)分布在了兩側(cè)，我們用prcomp()調(diào)用負(fù)載得分loading scores 的參數(shù)rotation
通過分析pca$rotation，可以了解該主成分與哪些基因相關(guān)性較高，哪些基因?qū)χ鞒煞值呢暙I(xiàn)較大。這對(duì)于解釋主成分分析的結(jié)果和理解數(shù)據(jù)的結(jié)構(gòu)和變化模式非常有幫助。
loading_scores <- pca$rotation[,1]
gene_score <- abs(loading_scores)
gene_score_ranked <- sort(gene_score,decreasing = TRUE)
top_10_genes <- names(gene_score_ranked[1:10])
pca$rotation[top_10_genes, 1]
這里我們查看 PC1 的負(fù)載得分，因?yàn)?PC1 解釋原始數(shù)據(jù)的 93.6%的方差
loading_scores <- pca$rotation[,1] 這行代碼的作用是將PCA分析結(jié)果中第一個(gè)主成分（即第一列）的負(fù)載得分（即100個(gè)基因數(shù)據(jù)）提取出來并賦值給loading_scores變量。
abs()取絕對(duì)值，再從小大到小排序，選取排名靠前的前十個(gè)基因top_10_genes
若想顯示每個(gè)基因?qū)?yīng)的 rotation，用以下代碼即可，展示的是帶正負(fù)值的
pca$rotation[top_10_genes, 1]
排名前10的基因在第一主成分上的負(fù)載得分文章來源：http://www.zghlxwxcb.cn/news/detail-507759.html
總的來說，本文對(duì)基因表達(dá)數(shù)據(jù)進(jìn)行了主成分分析，并可視化結(jié)果。
通過PCA主成分分析，可以降維并找到影響數(shù)據(jù)變化最大的主要因素，進(jìn)而進(jìn)行數(shù)據(jù)的可視化和分析。文章來源地址http://www.zghlxwxcb.cn/news/detail-507759.html
到了這里，關(guān)于R 語言 ggplot2 PCA 主成分分析（虛擬數(shù)據(jù)集）的文章就介紹完了。如果您還想了解更多內(nèi)容，請(qǐng)?jiān)谟疑辖撬阉鱐OY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！
