以第一個基因為例，根據group_list來檢驗在分組之間是否存在差異
load(file='step1-output.RData')
dat[1:4,1:4]
table(group_list) #分組信息
boxplot(dat[1,]~group_list)
t.test(dat[1,]~group_list) #p=0.00998有顯著性差異
library(ggpubr)
df <- data.frame(gene=dat[1,], stage=group_list) #比較下一個基因可以改為dat[2,]
p <- ggboxplot(df, x = "stage", y = "gene",
color = "stage", palette = "jco",
add = "jitter")
p+stat_compare_means() #p值和之前不一樣，因為換了一種統計學檢驗方法
以第一個基因為例進行表達差異性分析.Rplot==Note== : 第一個基因是隨機挑選的，雖然在兩個類群中有差異性，但是從圖上可以看出，noTNBC 有一部分是被包含在TNBC中的，並不是完全獨立分離的關係，統計學顯著性也不好說。
使用limma來進行批量的全部的基因的差異分析
#將繪製箱圖的函數包裝成函數便於使用
pb <- function(g){
library(ggpubr)
df <- data.frame(gene=g, stage=group_list) #比較下一個基因可以改為dat[2,]
p <- ggboxplot(df, x = "stage", y = "gene",
color = "stage", palette = "jco",
add = "jitter")
p+stat_compare_means() #p值和之前不一樣，因為換了一種統計學檢驗方法
}
pb(dat[1,])
library(limma)
design <- model.matrix(~factor(group_list)) #design就是limma做差異分析需要的一個輸入值，了解它可以有group_list構建
fit <- lmFit(dat, design)
fit <- eBayes(fit)
options(digits = 4)
topTable(fit,coef = 2,adjust="BH")
pb(dat["241662_x_at",]) #
deg <- topTable(fit, coef=2, adjust="BH", number = Inf) #Inf就是無窮大,把所有的數值都拿出來limma篩選出來的在兩個分組中差異表達比較顯著的基因, logFC為負表示下調，以noTNBC作為對照
limma識別到的差異表達基因以上面的第一個基因241662_x_at為例繪製箱線圖，這個基因在兩個分組之間的表達差異非常顯著，而且沒有重疊部分，TNBC和noTNBC完全分開了。
241662_x_at.Rplot下面我們需要將這些差異基因的基因名找出來，也就是將探針ID轉換到基因symbol
##進行注釋
rm(list=ls())
load(file = "deg.Rdata")
deg$probe_id = rownames(deg)
library(hgu133plus2.db)
?hgu133plus2.db #學習包的內容
ids <- toTable(hgu133plus2ENSEMBL) #取出探針和基因名對應的數據集，並將其轉化為表格
#merge函數，根據兩個數據框的相同的列名進行合併
deg_1 <- merge(deg,ids,by="probe_id")繪製火山圖
將上調和下調的基因清楚地顯示出來
nrDEG = deg_1
head(nrDEG)
attach(nrDEG)
plot(logFC,-log10(P.Value))
library(ggpubr)
df = nrDEG
df$v = -log10(P.Value)
ggscatter(df,x="logFC", y ='v', size = 0.5)
df$g = ifelse(df$P.Value > 0.01, "stable",
ifelse(df$logFC>1.5,"up",
ifelse(df$logFC < -1.5,"down","stable")))
table(df$g)
df$name = rownames(df)
ggscatter(df,x = "logFC", y = "v", size=0.5, color = "g")
ggscatter(df,x ="logFC",y="v",color = "g",size=0.5,
label = "symbol", repel = T,
label.select = c("PROM1","GABRP","AGR3","SCGB2A2")) #標記比較顯著的基因
火山圖1繪製熱圖
火山圖不需要表達矩陣，只要差異分析結果的表格就可以##繪製熱圖
load(file='step1-output.RData')
dat[1:4,1:4]
x = deg$logFC
names(x) = deg$probe_id
cg= c(names(head(sort(x),100)), names(tail(sort(x),100))) #分別挑出上調和下調表達最多的100個基因
library(pheatmap)
n <- t(scale(t(dat[cg,])))
n[n>2]=2
n[n< -2,]=-2
pheatmap(n,show_rownames = F,show_colnames = F)
ac=data.frame(g=group_list) #添加分組信息
rownames(ac)=colnames(n)
pheatmap(n, show_colnames = F,
show_rownames = F,
cluster_cols = F,
cluster_rows = F,
annotation_col = ac)使用前表達差異性特別顯著的基因繪製熱圖就可以很明顯地看出基因表達的模式差異
top100_rail100熱圖.Rplot01