【原】居然可以把rpkm這樣的歸一化并且?guī)?shù)點(diǎn)的轉(zhuǎn)錄組表達(dá)量矩陣直接取整

健明 2021-08-25

展開全文

rm(list = ls())  ## 魔幻操作，一鍵清空~
options(stringsAsFactors = F)#在調(diào)用as.data.frame的時(shí)，將stringsAsFactors設(shè)置為FALSE可以避免character類型自動轉(zhuǎn)化為factor類型
# 注意查看下載文件的大小，檢查數(shù)據(jù) 
f='GSE103611_eSet.Rdata'
# https://www.ncbi.nlm./geo/query/acc.cgi?acc=GSE103611
library(GEOquery)
# 這個(gè)包需要注意兩個(gè)配置，一般來說自動化的配置是足夠的。
#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)       ## 平臺文件
  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./geo/query/acc.cgi?acc=GSE103611 里面的是 Affymetrix Human Gene 2.0 ST Array [probe set (exon) version] 的表達(dá)量芯片。

但是對對 https://www.ncbi.nlm./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ū)分開來。

rpkm 格式的表達(dá)量矩陣可以區(qū)分兩個(gè)分組

學(xué)徒做了差異分析，然后上下調(diào)各自選取100個(gè)差異基因，熱圖可視化如下所示：

上下調(diào)各自選取100個(gè)差異基因

看起來分析的有模有樣！簡單了問了，才知道，學(xué)徒僅僅是根據(jù)我的表達(dá)芯片的公共數(shù)據(jù)庫挖掘系列推文：

解讀GEO數(shù)據(jù)存放規(guī)律及下載，一文就夠
解讀SRA數(shù)據(jù)庫規(guī)律一文就夠
從GEO數(shù)據(jù)庫下載得到表達(dá)矩陣一文就夠
GSEA分析一文就夠（單機(jī)版+R語言版）
根據(jù)分組信息做差異分析- 這個(gè)一文不夠的
差異分析得到的結(jié)果注釋一文就夠

代碼如下所示：

> 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包對 counts矩陣的差異分析結(jié)果
# deg_cptac  是DEseq2包對 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)了，

image-20210824102502370

兩次差異分析的結(jié)果居然是出奇的吻合，至少是從變化倍數(shù)的角度來看！

有意思

文末友情推薦

做教學(xué)我們是認(rèn)真的，如果你對我們的馬拉松授課（直播一個(gè)月互動教學(xué)）有疑問，可以看完我們從2000多個(gè)提問互動交流里面精選的200個(gè)問答！ 2021第二期_生信入門班_微信群答疑整理，以及 2021第二期_數(shù)據(jù)挖掘班_微信群答疑筆記

與十萬人一起學(xué)生信，你值得擁有下面的學(xué)習(xí)班：