本文主要目的是確定篩選的差異基因主要富集在哪些細胞中，可查看官方protocol https://github.com/NathanSkene/EWCE

install.packages("devtools")
library(devtools)
install_github("nathanskene/ewce")
#若安裝ewce時報錯，請看[http://m.itdecent.cn/p/40676b44a1b2](http://m.itdecent.cn/p/40676b44a1b2)
#，親測有效
library(EWCE)
library(ggplot2)
library(cowplot)
library(limma)
library(readxl)

首先，確定ctd數據集

#.Loading datasets
data("cortex_mrna")  # a list that contain an expression matrix and an annotation data frame
View(cortex_mrna)
cortex_mrna[[1]][1:10,1:5]
cortex_mrna[[2]][1:10,1:10]

image.png

image.png

image.png

image.png

由此，我們可以根據上面的數據格式導入我們自己的數據，即exp就是一個數據表達矩陣，而annot就是一個注釋矩陣，這個矩陣不需要像文中example那么多，但須要納入"cell_id", "level1class" and "level2class"

其次，對于公開可用的轉錄組數據集來說，使用不正確的基因符號是非常常見的(通常一些基因名稱在Excel中打開時會出現(xiàn)錯誤)。因此這個包提供了一個函數來糾正這類錯誤。第一次運行它時，需要從MGI下載MRK_List2文件，該文件列出了MGI基因符號的已知同義詞。建議在所有輸入數據集上運行此操作。

if(!file.exists("MRK_List2.rpt")){
    download.file("http://www.informatics.jax.org/downloads/reports/MRK_List2.rpt", destfile="MRK_List2.rpt")
}
cortex_mrna$exp = fix.bad.mgi.symbols(cortex_mrna$exp,mrk_file_path="MRK_List2.rpt")

##若為人，則下載這個：
##若為人，
if(!file.exists("HGNC_homologene.rpt")){
  download.file("http://www.informatics.jax.org/downloads/reports/HGNC_homologene.rpt", destfile="HGNC_homologene.rpt")
}
cortex_mrna$exp = fix.bad.mgi.symbols(cortex_mrna$exp,mrk_file_path="HGNC_homologene.rpt")

第三、Calculate specificity matrices

# Generate celltype data for just the cortex/hippocampus data
exp_CortexOnly_DROPPED = drop.uninformative.genes(exp=cortex_mrna$exp,level2annot = cortex_mrna$annot$level2class)
annotLevels = list(level1class=cortex_mrna$annot$level1class,level2class=cortex_mrna$annot$level2class)
ctd = generate.celltype.data(exp=exp_CortexOnly_DROPPED,annotLevels=annotLevels,groupName="kiCortexOnly")
print(ctd)

##最后一步是必要的，如果你計劃比較人類數據，即由人類遺傳學產生的基因集；這里這么說我推測可能是因為提供的example數據集是動物的
#ctd = filter.genes.without.1to1.homolog(fNames_CortexOnly)
#print(ctd)

load(fNames_CortexOnly[1])

image.png

第四步. Application to genetic data

data("example_genelist")
print(example_genelist)

image.png

由此，這里就可以 替換成我們自己的差異基因了

###若為小鼠
data("mouse_to_human_homologs")
m2h = unique(mouse_to_human_homologs[,c("HGNC.symbol","MGI.symbol")])
mouse.hits = unique(m2h[m2h$HGNC.symbol %in% example_genelist,"MGI.symbol"])
#mouse.bg  = unique(setdiff(m2h$MGI.symbol,mouse.hits))
mouse.bg  = unique(m2h$MGI.symbol)

image.png

##若為人，則human.hit and human.bg的建立方法
data("mouse_to_human_homologs")
#mouse.bg  = unique(setdiff(m2h$MGI.symbol,mouse.hits))
human.hits<- example_genelist
print(human.hits)
human.bg = unique(c(human.hits,m2h$HGNC.symbol))
length(human.bg)
head(human.bg)

image.png

reps=1000 # <- Use 1000 bootstrap lists so it runs quickly, for publishable analysis use >10000
level=1 # <- Use level 1 annotations (i.e. Interneurons)

# Bootstrap significance testing, without controlling for transcript length and GC content
full_results = bootstrap.enrichment.test(sct_data=ctd,hits=mouse.hits,bg=mouse.bg,
                                reps=reps,annotLevel=level)
print(full_results$results[order(full_results$results$p),3:5][1:6,])
print(ewce.plot(full_results$results,mtc_method="BH"))

image.png

image.png

如果您想查看列表中每個基因的富集特性，那么generate.bootstrap。應該使用plot函數。這將這些繪圖保存到BootstrapPlots文件夾中。

generate.bootstrap.plots(sct_data=ctd,hits=mouse.hits,bg=mouse.bg,reps=100,annotLevel=1,full_results=full_results,listFileName="VignetteGraphs")

image.png

第五、若數據為轉錄組數據，那么

data(tt_alzh)
tt_results = ewce_expression_data(sct_data=ctd,tt=tt_alzh,annotLevel=1,ttSpecies="human",sctSpecies="mouse")
ewce.plot(tt_results$joint_results)

image.png

image.png

這主要是閱讀官方protocol后自己的一點總結，若有不同意見，歡迎大家提出來一起探討。