不過，我這里不是要批評(píng)誰，我就是感興趣這個(gè)文章是如何使用這個(gè)數(shù)據(jù)集的，后面也有幾個(gè)小問題，感興趣的粉絲可以來信跟我交流哈。

文章介紹

轉(zhuǎn)移是腫瘤致死的原因之一。在腫瘤轉(zhuǎn)移的過程中，原發(fā)部位腫瘤細(xì)胞可以刺激骨髓來源的細(xì)胞（稱為被腫瘤馴化的細(xì)胞），促使并招募其進(jìn)入遠(yuǎn)端擬轉(zhuǎn)移部位，形成轉(zhuǎn)移前微環(huán)境，為腫瘤細(xì)胞轉(zhuǎn)移至該部位提供合適“土壤”。因此識(shí)別被腫瘤馴化的細(xì)胞，明確其功能是揭露腫瘤發(fā)生發(fā)展的關(guān)鍵步驟。

2019年1月15日凌晨，曹雪濤院士團(tuán)隊(duì)在Nature Medicine雜志發(fā)表了題為Tumor-educated B cells selectively promote breast cancer lymph node metastasis by HSPA4-targeting IgG的最新研究成果，詳細(xì)揭示了B細(xì)胞能夠通過分泌靶向腫瘤抗原HSPA4（Heat shock protein 4）的病理性抗體促進(jìn)乳腺癌淋巴結(jié)轉(zhuǎn)移。

本研究闡明了B細(xì)胞及抗體介導(dǎo)的體液免疫在淋巴結(jié)轉(zhuǎn)移前微環(huán)境形成及腫瘤淋巴結(jié)轉(zhuǎn)移中的重要功能。首次發(fā)現(xiàn)除調(diào)節(jié)性B細(xì)胞介導(dǎo)的負(fù)向免疫調(diào)控功能外，B細(xì)胞能夠通過分泌靶向腫瘤抗原的病理性抗體直接促進(jìn)腫瘤轉(zhuǎn)移，同時(shí)尋找到糖基化的腫瘤膜抗原在病理性抗體的產(chǎn)生及促轉(zhuǎn)移中的重要功能，為深入認(rèn)識(shí)體B細(xì)胞介導(dǎo)的體液免疫功能及腫瘤轉(zhuǎn)移前微環(huán)境的形成提供了新的視角，為腫瘤治療尤其是腫瘤轉(zhuǎn)移的預(yù)測(cè)和防治提供了新的靶點(diǎn)。

芯片實(shí)驗(yàn)設(shè)計(jì)及文章數(shù)據(jù)處理結(jié)果

值得注意的是研究者僅僅是測(cè)了才4個(gè)樣本，就區(qū)分成兩種細(xì)胞，兩種處理，也就是說根本沒有生物學(xué)重復(fù)。

從文章里面拿到了芯片信息，再看看其熱圖的描述：

# lymph nodes from normal mice (NLNs)
# draining lymph nodes (DLNs)

# cell cycle genes (such as Cdc25c, Bub1, Ttk, and Cdk1) and 
# migration-related genes (such as Vcam1, Arhgap5, Cxcr3 and Ccr2) 
# were upregulated in DLN B cells compared with NLN B cells 

# Heat map of genes expressed in purified B cells 
# (CD19+) from NLNs and DLNs as determined by gene chip array 
# (n = 1 for gene chip analysis of d?f)

是B細(xì)胞的數(shù)據(jù)的指定的兩個(gè)基因集，奇怪的是作者的圖例是 log2P值，很詭異！??！

再看看GSEA結(jié)果：

# d, Gene set enrichment analysis (GSEA) analysis for gene signatures of 
# leukocyte chemotaxis and humoral immune response 
# in purified stromal cells (CD45?) from DLN compared with NLN control. 
# NES, normalized enrichment score.

這里是stromal細(xì)胞的數(shù)據(jù)了。

作者繪制的是兩個(gè)GO基因集：

# GO:0002690 Positive regulation of leukocyte chemotaxis.
# http://amigo.geneontology.org/amigo/term/GO:0030595
# Accession GO:0030595
# Name leukocyte chemotaxis
# Ontology biological_process
# Synonyms immune cell chemotaxis, leucocyte chemotaxis


# http://amigo.geneontology.org/amigo/term/GO:0006959
# Accession GO:0006959
# Name humoral immune response
# Ontology biological_process

后面我們要重復(fù)這種圖就需要有這樣的背景知識(shí)，注意看圖中的基因集，應(yīng)該是50個(gè)基因左右的數(shù)量。

最后，作者列了一下排名前20的差異基因，如下：'

可以看到作者選取的差異倍數(shù)，而不是logFC

文章結(jié)果復(fù)現(xiàn)

表達(dá)矩陣的下載很容易：

rm(list = ls())  ## 魔幻操作，一鍵清空~
options(stringsAsFactors = F)
library(GEOquery)
# stromal cells from tumor-draining lymph nodes

gset <- getGEO('GSE113249', destdir=".",
               AnnotGPL = F,     ## 注釋文件
               getGPL = F)       ## 平臺(tái)文件  
a=gset[[1]] #
dat=exprs(a) #a現(xiàn)在是一個(gè)對(duì)象，取a這個(gè)對(duì)象通過看說明書知道要用exprs這個(gè)函數(shù)
dim(dat)#看一下dat這個(gè)矩陣的維度
head(dat)
pd=pData(a) #通過查看說明書知道取對(duì)象a里的臨床信息用pData
colnames(dat)=pd$title
dat1=dat

gset <- getGEO('GSE113250', destdir=".",
               AnnotGPL = F,     ## 注釋文件
               getGPL = F)       ## 平臺(tái)文件  
a=gset[[1]] #
dat=exprs(a) #a現(xiàn)在是一個(gè)對(duì)象，取a這個(gè)對(duì)象通過看說明書知道要用exprs這個(gè)函數(shù)
dim(dat)#看一下dat這個(gè)矩陣的維度
head(dat)
pd=pData(a) #通過查看說明書知道取對(duì)象a里的臨床信息用pData
colnames(dat)=pd$title
dat2=dat

dat=cbind(dat1,dat2)

dat[1:4,1:4] 
boxplot(dat)
# GPL21163    Agilent-074809 SurePrint G3 Mouse GE v2 8x60K Microarray [Probe Name version]
if(F){
  library(GEOquery)
  #Download GPL file, put it in the current directory, and load it:
  gpl <- getGEO('GPL21163', destdir=".")
  colnames(Table(gpl))  
  head(Table(gpl)[,c(1,6)]) ## you need to check this , which column do you need
  probe2gene=Table(gpl)[,c(1,6)] 
  save(probe2gene,file='probe2gene.Rdata')
}

load(file='probe2gene.Rdata')
ids=probe2gene 
head(ids)
colnames(ids)=c('probe_id','symbol')  
ids=ids[ids$symbol != '',]
ids=ids[ids$probe_id %in%  rownames(dat),]

dat[1:4,1:4]   
dat=dat[ids$probe_id,] 

ids$median=apply(dat,1,median) #ids新建median這一列，列名為median，同時(shí)對(duì)dat這個(gè)矩陣按行操作，取每一行的中位數(shù)，將結(jié)果給到median這一列的每一行
ids=ids[order(ids$symbol,ids$median,decreasing = T),]#對(duì)ids$symbol按照ids$median中位數(shù)從大到小排列的順序排序，將對(duì)應(yīng)的行賦值為一個(gè)新的ids
ids=ids[!duplicated(ids$symbol),]#將symbol這一列取取出重復(fù)項(xiàng)，'!'為否，即取出不重復(fù)的項(xiàng)，去除重復(fù)的gene ，保留每個(gè)基因最大表達(dá)量結(jié)果s
dat=dat[ids$probe_id,] #新的ids取出probe_id這一列，將dat按照取出的這一列中的每一行組成一個(gè)新的dat
rownames(dat)=ids$symbol#把ids的symbol這一列中的每一行給dat作為dat的行名
dat[1:4,1:4]  #保留每個(gè)基因ID第一次出現(xiàn)的信息
dim(dat)

save(dat,file = 'step1-output.Rdata')

load(file = 'step1-output.Rdata')

可以看到樣本相關(guān)性與其分組是匹配的：

熱圖很容易重現(xiàn)，代碼如下：

rm(list = ls())  ## 魔幻操作，一鍵清空~
options(stringsAsFactors = F)
library(GEOquery)
load(file = 'step1-output.Rdata')
colnames(dat)
dat[1:4,1:4]

library(pheatmap)
ccg=trimws(strsplit('Cdc25c, Bub1, Ttk,  Cdk1',',')[[1]])
mrg=trimws(strsplit('Vcam1, Arhgap5, Cxcr3, Ccr2',',')[[1]])
mat=dat[c(ccg,mrg),3:4]
mat
pheatmap(mat)
pheatmap(dat[rownames(dat)[grepl('^Ig',rownames(dat))],])

主要是GSEA圖，我首先拿出基因排序：

rm(list = ls())  ## 魔幻操作，一鍵清空~
options(stringsAsFactors = F)
library(GEOquery)
load(file = 'step1-output.Rdata')
colnames(dat)
dat[1:4,1:4]

rownames(dat)=toupper(rownames(dat))

m=dat[,2]+dat[,1]/2
logFC1=log2(dat[,2]/dat[,1]);fivenum(logFC1)
logFC2=dat[,2]-dat[,1];fivenum(logFC2)
plot(logFC1,m)
plot(logFC2,m)
plot(logFC1,logFC2)

# 這里的 logFC1 和 logFC2 差異很小

geneList=scale( logFC2 )r
# geneList=scale( dat[,2]/dat[,1] )
names(geneList)=rownames(dat)
geneList=sort(geneList,decreasing = T)
head(geneList)

所有的 gsea 分析都應(yīng)該是基于 geneList這個(gè)對(duì)象，這個(gè)編程思想很重要。

然后創(chuàng)建基因集，我首先是要 gskb 這個(gè)R包

library(gskb)
#browseVignettes("gskb")  
library(gskb) 
data(mm_GO)
(g1=mm_GO$GO_BP_MM_HUMORAL_IMMUNE_RESPONSE)
(g2=mm_GO$GO_BP_MM_LEUKOCYTE_CHEMOTAXIS)
g1=g1[g1 %in% rownames(dat)]
g2=g2[g2 %in% rownames(dat)]
library(pheatmap)
pheatmap(dat[g1 ,1:2])
pheatmap(dat[g2 ,1:2])

然后進(jìn)行GSEA分析：

library(ggplot2)
library(clusterProfiler)
library(GSEABase) # BiocManager::install('GSEABase')
# geneset <- read.gmt('h.all.v6.2.symbols.gmt')
# Then we know geneset is a simple data.frame
# first column is name of geneset, second column is gene symbol
geneset=data.frame(ont=c(rep('HUMORAL_IMMUNE_RESPONSE',length(g1)),
                         rep('LEUKOCYTE_CHEMOTAXIS',length(g2))),
                   gene=c(g1,g2))
egmt <- GSEA(geneList, TERM2GENE=geneset, 
             verbose=T)
head(egmt)
gseaplot(egmt,'HUMORAL_IMMUNE_RESPONSE') 
gseaplot(egmt,'LEUKOCYTE_CHEMOTAXIS')

發(fā)現(xiàn)只有一個(gè)基因集是顯著的，另外一個(gè)慘不忍睹，我就不發(fā)圖了。

然后我嘗試使用 GO.db 包的數(shù)據(jù)集。

library(GO.db)  
ls("package:GO.db")  
library(org.Hs.eg.db)
eg2go=toTable(org.Hs.egGO2ALLEGS)
head(eg2go)
tmp=eg2go[eg2go$go_id=='GO:0030595',]
columns(org.Hs.eg.db)
g11=select(org.Hs.eg.db,tmp$gene_id,'SYMBOL', 'ENTREZID')[,2]
g11=g11[g11 %in% rownames(dat)]

tmp=eg2go[eg2go$go_id=='GO:0006959',]
columns(org.Hs.eg.db)
g22=select(org.Hs.eg.db,tmp$gene_id,'SYMBOL', 'ENTREZID')[,2]
g22=g22[g22 %in% rownames(dat)]
library(pheatmap)
pheatmap(dat[g11 ,1:2])
pheatmap(dat[g22 ,1:2])
geneset=data.frame(ont=c(rep('HUMORAL_IMMUNE_RESPONSE',length(g11)),
                         rep('LEUKOCYTE_CHEMOTAXIS',length(g22))),
                   gene=c(g11,g22))
egmt <- GSEA(geneList, TERM2GENE=geneset, 
             verbose=T)
head(egmt)
gseaplot(egmt,'HUMORAL_IMMUNE_RESPONSE') 
gseaplot(egmt,'LEUKOCYTE_CHEMOTAXIS')

發(fā)現(xiàn)也不是顯著富集，最后我再使用Y叔的R包，這里面有個(gè)R代碼技巧，很重要，值得大家花費(fèi)幾個(gè)小時(shí)慢慢理解。

rm(list = ls())  ## 魔幻操作，一鍵清空~
options(stringsAsFactors = F)
library(GEOquery)
load(file = 'step1-output.Rdata')
colnames(dat)
dat[1:4,1:4]

rownames(dat)=toupper(rownames(dat))

m=dat[,2]+dat[,1]/2
logFC1=log2(dat[,2]/dat[,1]);fivenum(logFC1)
logFC2=dat[,2]-dat[,1];fivenum(logFC2)
plot(logFC1,m)
plot(logFC2,m)
plot(logFC1,logFC2)

geneList=scale( logFC2 )
names(geneList)=rownames(dat)
geneList=sort(geneList,decreasing = T)
head(geneList)

s2g=select(org.Hs.eg.db,names(geneList),
           'ENTREZID','SYMBOL')

head(s2g)
s2g=s2g[!is.na(s2g$ENTREZID),]
s2g=s2g[s2g$SYMBOL %in%  names(geneList),]
geneList=geneList[s2g$SYMBOL]  
names(geneList)=s2g$ENTREZID  
head(geneList)
tail(geneList)

library(ggplot2)
library(clusterProfiler)
library(GSEABase) 
## 這里比較耗時(shí)：
go_bp_gsea <- gseGO(geneList     = geneList, 
             OrgDb='org.Hs.eg.db',
                    ont = 'BP',
                     nPerm        = 1000,
                     minGSSize    = 10,
                     pvalueCutoff = 0.9,
                     verbose      = FALSE)

gseaplot(go_bp_gsea, geneSetID = rownames(go_bp_gsea[1,]))
gseaplot(go_bp_gsea, geneSetID = "GO:0006959")
gseaplot(go_bp_gsea, geneSetID = "GO:0030595")

得到的GSEA圖，仍然是不顯著。

思考

為什么我工具芯片表達(dá)量算的兩種logFC的散點(diǎn)圖出現(xiàn)了染色體一樣的結(jié)構(gòu)：

為什么使用不同數(shù)據(jù)庫來源的同一個(gè)生物學(xué)功能基因集得到的結(jié)果可以千差萬別。

為什么我最后也無法重復(fù)出來一個(gè)基因集的顯著結(jié)果呢？是作者錯(cuò)了嗎，還是我分析方法不對(duì)？

提示

我默認(rèn)對(duì)同一個(gè)基因?qū)?yīng)的多個(gè)探針，選取了表達(dá)量最大的探針，但是我在檢查作者的top20基因的時(shí)候發(fā)現(xiàn)，特別詭異的是我看到的那些基因表達(dá)量并沒有差異，所以我懷疑是選取了表達(dá)量最大的探針來進(jìn)行過濾這個(gè)策略恐怕是有問題的。

比如 Anapc7 這個(gè)基因是有3個(gè)探針的，但是：

# 
# A_52_P384134 Anapc7 
# A_55_P2112330 Anapc7 
# A_55_P2731736 Anapc7

dat['A_52_P384134',]
dat['A_55_P2112330',]
dat['A_55_P2731736',]

我恰好挑的探針是沒有差異的，如下：

這樣問題就來了，作者只測(cè)一個(gè)樣本，這樣的數(shù)值可靠嗎？

本站僅提供存儲(chǔ)服務(wù)，所有內(nèi)容均由用戶發(fā)布，如發(fā)現(xiàn)有害或侵權(quán)內(nèi)容，請(qǐng)點(diǎn)擊舉報(bào)。

免费视频淫片aa毛片_日韩高清在线亚洲专区vr_日韩大片免费观看视频播放_亚洲欧美国产精品完整版

文章介紹

芯片實(shí)驗(yàn)設(shè)計(jì)及文章數(shù)據(jù)處理結(jié)果

文章結(jié)果復(fù)現(xiàn)

思考

提示