rm(list = ls()) ## 魔幻操作，一鍵清空~
options(stringsAsFactors = F)#在調(diào)用as.data.frame的時(shí)，將stringsAsFactors設(shè)置為FALSE可以避免character類型自動(dòng)轉(zhuǎn)化為factor類型
# 注意查看下載文件的大小，檢查數(shù)據(jù)
f='GSE103611_eSet.Rdata'
# https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE103611
library(GEOquery)
# 這個(gè)包需要注意兩個(gè)配置，一般來說自動(dòng)化的配置是足夠的。
#Setting options('download.file.method.GEOquery'='auto')
#Setting options('GEOquery.inmemory.gpl'=FALSE)
if(!file.exists(f)){
gset <- getGEO('GSE103611', destdir=".",
AnnotGPL = F, ## 注釋文件
getGPL = F) ## 平臺(tái)文件
save(gset,file=f) ## 保存到本地
}
load('GSE103611_eSet.Rdata') ## 載入數(shù)據(jù)
class(gset) #查看數(shù)據(jù)類型
length(gset) #
class(gset[[1]])
gset
# assayData: 352859 features, 48 samples
是可以拿到表達(dá)量矩陣，因?yàn)閔ttps://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE103611 里面的是 Affymetrix Human Gene 2.0 ST Array [probe set (exon) version] 的表達(dá)量芯片。
但是對(duì) 對(duì) https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE106292 上面的代碼就拿不到表達(dá)矩陣 ，因?yàn)?是轉(zhuǎn)錄組測序數(shù)據(jù)。其 提供的表達(dá)量矩陣是 GSE106292_Human_Matrix_final.xlsx這個(gè) 4.6 Mb，就很麻煩，因?yàn)樗⒉皇菢?biāo)準(zhǔn)的counts矩陣，不一定適合于 edgeR、DEseq2這樣的包！
但是最近給學(xué)徒安排了一個(gè)類似的僅僅是提供了rpkm這樣的歸一化并且?guī)?shù)點(diǎn)的轉(zhuǎn)錄組表達(dá)量矩陣項(xiàng)目做差異分析，發(fā)現(xiàn)了一個(gè)騷操作，我也是醉了。如下所示，rpkm矩陣仍然是可以把腫瘤和正常組織涇渭分明的區(qū)分開來。
學(xué)徒做了差異分析，然后 上下調(diào)各自選取100個(gè)差異基因，熱圖可視化如下所示：
看起來分析的有模有樣！簡單了問了，才知道，學(xué)徒僅僅是根據(jù)我的表達(dá)芯片的公共數(shù)據(jù)庫挖掘系列推文 ：
代碼如下所示：
> mat[1:4,1:3] # rpkm 格式的表達(dá)量矩陣
C3L-00418_Solid Tissue Normal C3L-00004_Primary Tumor C3N-00852_Primary Tumor
OR4F5 0.0000000 0.0000000 0.0000000
SAMD11 4.0715358 0.3719551 0.6190211
KLHL17 2.3878834 2.1623200 1.6117428
PLEKHN1 0.3629224 0.4835416 1.5632905
> exprSet=log2(mat+1)
> # 為了使用edgeR、DEseq2
> exprSet <- ceiling(exprSet)
> exprSet[1:4,1:3]
C3L-00418_Solid Tissue Normal C3L-00004_Primary Tumor C3N-00852_Primary Tumor
OR4F5 0 0 0
SAMD11 3 1 1
KLHL17 2 2 2
PLEKHN1 1 1 2
也就是說，學(xué)徒僅僅是把下載好了的 rpkm 格式的表達(dá)量矩陣進(jìn)行了log2和取整的操作， 然后就假裝它是一個(gè) 整數(shù)矩陣，去走 edgeR、DEseq2這樣的專門為轉(zhuǎn)錄組測序counts矩陣開發(fā)的差異分析流程！
我讓學(xué)徒重新走轉(zhuǎn)錄組測序數(shù)據(jù)分析流程，拿到真正的counts矩陣，再做差異分析，然后比較兩次差異分析結(jié)果。
#畫九宮格就是上下調(diào)基因畫在一起，這樣用logFC=1畫
#～～～～～～數(shù)據(jù)進(jìn)行比較～～～～～～～
# deg_gdc 是DEseq2包對(duì) counts矩陣的差異分析結(jié)果
# deg_cptac 是DEseq2包對(duì) rpkm 矩陣的差異分析結(jié)果
int_gene=intersect(deg_gdc$V1,deg_cptac$V1)
head(int_gene)
length(int_gene)# 17568
comp=cbind(deg_gdc[int_gene,],
deg_cptac[int_gene,])
head(comp)
dim(comp)
colnames(comp)
#～～～～～看上下調(diào)基因的交集～～～～～
table(comp[,c(8,16)])
# change.1
# change DOWN stable UP
# DOWN 611 1568 0
# stable 11 12647 1
# UP 0 2387 343
#～～～～～初步畫圖～～～～～
plot(comp[,c(3,11)])
dev.off()
#～～～～進(jìn)階畫圖～～～～
comp_logFC <- comp[,c(3,8,11,16)]
head(comp_logFC)
logFC_t = 1
#P.Value_t = 0.05
p <- ggplot(comp_logFC, aes(x=comp_logFC$log2FoldChange, y=comp_logFC$log2FoldChange.1))+
geom_point()+
labs(x = "GDC",
y = "CPTAC")+
scale_color_manual(values=c("blue", "grey","red"))+
geom_vline(xintercept = c(-logFC_t,logFC_t),lty=4,col="#ec183b",lwd=0.8) +
geom_hline(yintercept = c(-logFC_t,logFC_t),lty=4,col="#18a5ec",lwd=0.8) +
xlim(-5,5)+
ylim(-3,3)+
theme_ggstatsplot()+
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank())
p
ggsave(filename = 'sudoku.png',width = 10,height = 8)
dev.off()
有意思的地方出現(xiàn)了，
兩次差異分析的結(jié)果居然是出奇的吻合，至少是從變化倍數(shù)的角度來看！
有意思
做教學(xué)我們是認(rèn)真的，如果你對(duì)我們的馬拉松授課（直播一個(gè)月互動(dòng)教學(xué)）有疑問，可以看完我們從2000多個(gè)提問互動(dòng)交流里面精選的200個(gè)問答！ 2021第二期_生信入門班_微信群答疑整理，以及 2021第二期_數(shù)據(jù)挖掘班_微信群答疑筆記
與十萬人一起學(xué)生信，你值得擁有下面的學(xué)習(xí)班：
聯(lián)系客服
