進行基因差異分析
1.讀取數據
rm(list = ls()) ## 魔幻操作醋火，一鍵清空~
options(stringsAsFactors = F)#在調用as.data.frame的時益老，將stringsAsFactors設置為FALSE可以避免character（字符）類型自動轉化為factor（因子，如0那先，1）類型
# 注意查看下載文件的大小，檢查數據
load(file = 'step1-output.Rdata')#加載第一步保存的文件赚爵，此文件包括dat和group_list庶柿，加載后得到此兩數據
dat
group_list
2.檢測數據
# 每次都要檢測數據
dat[1:4,1:4] #查看1到4行和1到4列的dat文件
table(group_list) #table函數村怪，查看group_list中的分組個數
dat[1:4,1:4]
table(group_list)
3.為每個數據集繪制箱形圖，比較數據集中的數據分布（單一某個基因分析）
boxplot(dat[1,]~group_list) #按照group_list分組畫箱線圖,定義分組澳泵，（此步驟可以省略）
bp=function(g){ #定義一個函數g实愚，函數為{}里的內容
library(ggpubr)
df=data.frame(gene=g,stage=group_list)#把分組和數據糅合在一起
p <- ggboxplot(df, x = "stage", y = "gene",#定義圖形的橫縱坐標
color = "stage", palette = "jco",
add = "jitter")#
p + stat_compare_means()# Add p-value值
}
bp(dat[1,]) ## 調用上面定義好的函數，避免同樣的繪圖代碼重復多次敲兔辅。
bp(dat[2,])
bp(dat[3,])
bp(dat[4,])
dim(dat)
4.對某種基因在樣本中的表達情況根據分組統(tǒng)一分析（所有基因一起分析）
4.1數據表達形式
用limma包，這里注意击喂，limma包是對基因芯片表達矩陣的分析维苔，不能對逆轉錄RNAseq表達矩陣進行分析（因為數據特征不同），RNAseq需要用另一種方法
library(limma)
design=model.matrix(~factor( group_list ))#把分組信息變化成因子
design
fit=lmFit(dat,design)#在給定一系列陣列的情況下懂昂，擬合每個基因的線性模型
QQ截圖20190326210520.jpg
fit=eBayes(fit)#eBayes函數指在微陣列線性模型擬合下介时，通過經驗Bayes對標準誤差向一個共同值的緩和，計算緩和t-統(tǒng)計凌彬、緩和f-統(tǒng)計和微分表達式的對數概率沸柔。 用法舉列：EBAYES（配合，比例=0.01铲敛，標準偏差褐澎，系數，極限值=C（0.1,4）
## 上面是limma包用法的一種方式
#把上面計算出的差異再進行處理
options(digits = 4) #設置全局的數字有效位數為4
#topTable(fit,coef=2,adjust='BH')
topTable(fit,coef=2,adjust='BH') #提取一個表的列上的基因從一個線性模型的擬合
#例子：topTable(fit, coef=NULL, number=10, genelist=fit$genes, adjust.method="BH",
#sort.by="B", resort.by=NULL, p.value=1, lfc=0, confint=FALSE)
所有基因差異總數據表
解讀此表
表中數據意義
但是上面的用法做不到隨心所欲的指定任意兩組進行比較伐蒋，所有還有下一種方法
design <- model.matrix(~0+factor(group_list))#把分組變成因子形式
以分組信息為因子得到design
colnames(design)=levels(factor(group_list))#改design列名為分組信息
改之后列名變了
head(design)#查看分組文件
exprSet=dat#exprSet 為表達矩陣
rownames(design)=colnames(exprSet)#更改分組信息design行名為分組名工三，也就是樣本名字
design
#構造了design矩陣迁酸，得出每一個樣本屬于哪一個組，屬于為1俭正，不屬于為0
QQ截圖20190326213059.jpg
處理好了分組信息奸鬓，再自定義比較元素
#自定義比較元素
contrast.matrix<-makeContrasts("noTNBC-TNBC",
levels = design)
#也可以用下面的，和上面等同意思
#contrast.matrix<-makeContrasts("paste0(unique(group_list),collapse = "-")掸读，levels = design)
contrast.matrix ##這個矩陣聲明串远，我們要把 Tumor 組跟 Normal 進行差異分析比較
自定義函數進行比較
deg = function(exprSet,design,contrast.matrix){#定義一個叫deg的函數
##step1#logFC后是前的多少倍,進行差異化比較
fit <- lmFit(exprSet,design)
##step2
fit2 <- contrasts.fit(fit, contrast.matrix)
##這一步很重要，大家可以自行看看效果
fit2 <- eBayes(fit2) ## default no trend !!!
##eBayes() with trend=TRUE
##step3
tempOutput = topTable(fit2, coef=1, n=Inf)
nrDEG = na.omit(tempOutput)
#write.csv(nrDEG2,"limma_notrend.results.csv",quote = F)
head(nrDEG)
return(nrDEG)
}
deg = deg(exprSet,design,contrast.matrix)
head(deg)#查看前幾行文件
save(deg,file = 'deg.Rdata')#存儲文件
#得出了兩兩差異化比較的圖
得出所有基因差異化數據表
4.2分析基因表達差異的可視化表示：火山圖儿惫，熱圖
for volcano 火山圖
if(T){
nrDEG=deg
head(nrDEG)
attach(nrDEG)
plot(logFC,-log10(P.Value))
library(ggpubr)
df=nrDEG
df$v= -log10(P.Value) #df新增加一列'v',值為-log10(P.Value)
ggscatter(df, x = "logFC", y = "v",size=0.5)
df$g=ifelse(df$P.Value>0.01,'stable
', #if 判斷：如果這一基因的P.Value>0.01抑淫，則為stable基因
ifelse( df$logFC >2,'up', #接上句else 否則：接下來開始判斷那些P.Value<0.01的基因，再if 判斷：如果logFC >1.5,則為up（上調）基因
ifelse( df$logFC < -2,'down','stable') )#接上句else 否則：接下來開始判斷那些logFC <1.5 的基因姥闪，再if 判斷：如果logFC <1.5始苇，則為down（下調）基因，否則為stable基因
)
table(df$g)
df$name=rownames(df)
head(df)
ggscatter(df, x = "logFC", y = "v",size=0.5,color = 'g')
ggscatter(df, x = "logFC", y = "v", color = "g",size = 0.5,
label = "name", repel = T,
#label.select = rownames(df)[df$g != 'stable'] ,
label.select = c('TTC9', 'AQP3', 'CXCL11','PTGS2'), #挑選一些基因在圖中顯示出來
palette = c("#00AFBB", "#E7B800", "#FC4E07") )
ggsave('volcano.png')
ggscatter(df, x = "AveExpr", y = "logFC",size = 0.2)
df$p_c = ifelse(df$P.Value<0.001,'p<0.001',
ifelse(df$P.Value<0.01,'0.001<p<0.01','p>0.01'))
table(df$p_c )
ggscatter(df,x = "AveExpr", y = "logFC", color = "p_c",size=0.2,
palette = c("green", "red", "black") )
ggsave('MA.png')
}
for heatmap
if(T){
load(file = 'step1-output.Rdata')
# 每次都要檢測數據
dat[1:4,1:4]
table(group_list)
x=deg$logFC #deg取logFC這列并將其重新賦值給x
names(x)=rownames(deg) #deg取probe_id這列筐喳，并將其作為名字給x
cg=c(names(head(sort(x),100)),#對x進行從小到大排列催式，取前100及后100，并取其對應的探針名避归，作為向量賦值給cg
names(tail(sort(x),100)))
library(pheatmap)
pheatmap(dat[cg,],show_colnames =F,show_rownames = F) #對dat按照cg取行荣月，所得到的矩陣來畫熱圖
n=t(scale(t(dat[cg,])))#通過“scale”對log-ratio數值進行歸一化，現在的dat是行名為探針梳毙，列名為樣本名哺窄，由于scale這個函數應用在不同組數據間存在差異時，需要行名為樣本账锹，因此需要用t(dat[cg,])來轉換萌业，最后再轉換回來
n[n>2]=2
n[n< -2]= -2
n[1:4,1:4]
pheatmap(n,show_colnames =F,show_rownames = F)
ac=data.frame(g=group_list)
rownames(ac)=colnames(n) #將ac的行名也就分組信息（是‘no TNBC’還是‘TNBC’）給到n的列名，即熱圖中位于上方的分組信息
pheatmap(n,show_colnames =F,
show_rownames = F,
cluster_cols = F,
annotation_col=ac,filename = 'heatmap_top200_DEG.png') #列名注釋信息為ac即分組信息
}
write.csv(deg,file = 'deg.csv')
dev.new()
熱土和火山圖都是傻瓜式的奸柬，只要的前面得出的deg數據(也就是基因差異表達數據）是正確的
