本文是用主成分分析的方法去預測未知樣本的信息。
這裡有一批179個樣本的數據，正常組織27個，腫瘤組織129個，轉移組織18個，還有5個組織因為師妹不小心，標籤弄掉了 所以不知道信息。
那麼本次就是要通過主成分分析的方法，預測那5個樣本的信息。
第一，加載數據增加一些分組信息，比如已知的樣本分類，訓練組，測試組。
load(file = "taylor2010.Rdata")
index_normal <- grepl("N.P", colnames(taylor))
index_primary <- grepl("P.P", colnames(taylor))
index_met <- grepl("M.P", colnames(taylor))
n_normal <- sum(index_normal);
n_primary = sum(index_primary);
n_met = sum(index_met);
# class label (design vector)
taylor_classes = c(rep(0,n_normal), rep(1,n_primary), rep(2,n_met));
# train (known type samples), and test (unknown type samples)
train <- taylor[,1:174];
test <- taylor[,175:179];
# colors for plotting
cols = c(taylor_classes+2, 1,1,1,1,1)
tumortype_class <- factor(taylor_classes, levels = 0:2,
labels = c("Normal", "Primary", "Metastasized"))
train_samps <- 1:174
test_samps <- 175:179
這些數據已經被標準化處理過。
第二，批次處理在之前，我們需要看一下他有沒有批次效應
edata <- train
dist_mat <- dist(t(edata))
clustering <- hclust(dist_mat, method = "complete")
#plot(clustering, labels = tumortype_class)
plot(clustering, labels = taylor_classes)
這個明顯是有批次效應的，上次的那個帖子中
批次效應這樣矯正
是給出批次信息的，而這次沒有， 當沒有批次信息的時候，眼睛看也可以大致判定批次(如果沒有其他因素影響) 我們在這個帖子裡面就做過 批次效應有時候真麻煩！
但是，本次樣本量很多，我並不能判斷，所以這次不作批次處理，並不影響我們的下遊分析，
第三，差異分析。這裡的方法我們在這個帖子裡面也講過 Limma 求差異基因構建矩陣的兩種方式 這裡我想知道癌和癌旁中的差異基因，轉移和癌旁中的差異基因，轉移和癌的差異基因，所以選擇多組矩陣構建
library(limma)
f <- tumortype_class
# 分組矩陣
design <- model.matrix(~0+f)
colnames(design) <- levels(f)
fit <- lmFit(train, design)
# 比較矩陣
contrast.matrix <- makeContrasts(Primary-Normal,
Metastasized-Normal,
Metastasized-Primary,
levels=design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
把三個組的差異基因分別提取出來
PvsN <- topTable(fit2, coef=1,number = Inf,p.value=0.05, lfc=1.5)
MvsN <- topTable(fit2, coef=2,number = Inf,p.value=0.05, lfc=1.5)
MvsP <- topTable(fit2, coef=3,number = Inf,p.value=0.05, lfc=1.5)
length(rownames(PvsN))
length(rownames(MvsN))
length(rownames(MvsP))
分別獲得了，14, 546，189個差異基因，為什麼癌和癌旁的差異那麼大，第一，他們本身就有差異， 第二，因為我們設置了fdr過濾，所以樣本越大越有可能挑選出穩定差異的基因。
把三組差異基因合併起來，得到562個，說明，這幾個之間有重合。
set12 <-union(rownames(PvsN),rownames(MvsN))
set123 <-union(set12,rownames(MvsP))
length(set123)
第四，主成分分析PCAR中prcomp可以直接做主成分分析
edata <- taylor[set123,]
pca_raw <- prcomp(t(edata), center = TRUE, scale. = TRUE)
用biplot作圖
biplot(pca_raw)
可能是我年紀比較小吧，還欣賞不來這種藝術形式，所以我決定自己畫個圖
我不僅畫，還要把那些樣本的信息標註出來
edata_pc_df <- as.data.frame(pca_raw$x)
taylor_classes = c(rep("Normal",n_normal),
rep("Primary",n_primary),
rep("Metastasized",n_met),
rep("Unknown",5))
taylor_classes <- factor(taylor_classes,
levels = c("Normal", "Primary", "Metastasized","Unknown"))
edata_pc_df$sample <- taylor_classes
#head(edata_pc_df)
library(ggplot2)
library(ggrepel)
edata_pc_df$samplenames <- rownames(edata_pc_df)
ggplot(edata_pc_df, aes(x = PC1, y = PC2, color = taylor_classes)) +
geom_point() +
geom_text_repel(data=subset(edata_pc_df, edata_pc_df$sample=="Unknown"),
aes(label=samplenames),col="black",alpha = 0.8)
圖中紫色的點就代表位置樣本，我已經用箭頭指出來了，我們發現，依然三個樣本屬於轉移組，剩下的兩個 混跡在normal組中，肉眼可見的差異才叫差異嘛！
過了幾天，我又看到了另外一個包，可以直接對PCA結果可視化，我就試用了一下。
首先這個包的安裝需要devtools，然後要從github上安裝，說實話我沒有安裝成功。 如果你像我一樣，請你轉移到公眾號裡面的那篇，R包的無敵安裝，可以解決任意R包的安裝。 估計洲更看到這裡會跳起來打我，但是，沒有用的，有時候我就是這麼流氓。
install.packages("devtools")
library(devtools)
install_github("vqv/ggbiplot")
安裝的困難，用起來簡單
# PCA
edata <- taylor[set123,]
pca_raw <- prcomp(t(edata), center = TRUE, scale. = TRUE)
# 作圖
library(ggbioplot)
ggbiplot(pca_raw, obs.scale = 1, var.scale = 1,
groups = taylor_classes, ellipse = T, circle = F,var.axes = F)
這時候我們發現兩個normal要實錘了，但是那三個轉移的反倒是沒那麼可信了。
肯定有人會不理解，為什麼要先算差異基因呢，最起碼有兩個原因： 第一，也是最重要的，差異基因在本質上校正了批次效應，排除了很多背景。 第二，是我自己發現的，實際上不重要，SVM還有回歸這些算法不支持大量基因。 感興趣的翻看到公眾號的SVM那個帖子看看就知道了，他也是要接受差異基因的。
那麼在這裡，我們說幹就幹，把全基因帶入進來，看看PCA還能幹什麼
edata <- taylor
pca_raw <- prcomp(t(edata), center = TRUE, scale. = TRUE)
ggbiplot(pca_raw, obs.scale = 1, var.scale = 1,
groups = taylor_classes, ellipse = T, circle = F,var.axes = F)
主成分：我的特長就是抓重點，但是首先你得有重點。
總結：主成分分析也是能夠預測未知分類的，但是和K均值一樣需要有眼光。 不過，不要沮喪，沒關係啊，我們還有，支持向量機SVM，和K鄰近算法呢，
Kmeans 和 PCA 靠看，SVM 和 KNN 靠算，靠看容易看花眼，靠算長留天地間。