寫在前面：今天學(xué)習(xí)了Rcis Target這個(gè)包，每一次在學(xué)習(xí)一個(gè)包之前，苦惱和痛苦是不知道這個(gè)R包用來do for what？why to do？今天看了http://m.itdecent.cn/p/6e1d71db4220周老師的簡書，感覺講的就是我自己，下面我將用我的生物學(xué)知識來講關(guān)于Rcis Target的理解，假如有錯(cuò)的地方歡迎大家指正與批評?。?！
前言：我們在單細(xì)胞分析的時(shí)候，通過差異基因來判斷哪些clusters是什么類型的細(xì)胞，所以每一種細(xì)胞類型是一些基因（即基因集）選擇表達(dá)的結(jié)果，而基因的表達(dá)是受轉(zhuǎn)錄調(diào)控的（分子生物學(xué)知識），單從轉(zhuǎn)錄因子（TF）這一個(gè)方面來說，不同的轉(zhuǎn)錄因子調(diào)控不同的基因，進(jìn)而使細(xì)胞表現(xiàn)出不同的狀態(tài)或類型。當(dāng)然轉(zhuǎn)錄調(diào)控不只有TF調(diào)控 還有很多其他的如：順式作用元件。
1.do for what：用B細(xì)胞為例，我們可以對B細(xì)胞中高度相關(guān)的基因集（即：讓這些細(xì)胞為B細(xì)胞的基因）進(jìn)行Rcis Target分析，就可以知道哪一些TF在B細(xì)胞是富集或特有的，籠統(tǒng)的說Rcis Target分析就是用來預(yù)測基因集的轉(zhuǎn)錄因子。（還可以這樣運(yùn)用，即：當(dāng)發(fā)現(xiàn)tumor中的T細(xì)胞數(shù)比normal中的T細(xì)胞數(shù)明顯多，我們可以將tumor的T細(xì)胞與normal的T細(xì)胞進(jìn)行比較，通過篩選過濾找出二者具有統(tǒng)計(jì)學(xué)差異的基因集，進(jìn)行Rcis Target分析，就可以看出是什么TF導(dǎo)致了tumor和normal中T細(xì)胞的差異）
需要三個(gè)文件，基因集，Gene-motif rankings::提供每個(gè)motif的所有基因的排名(~得分)，轉(zhuǎn)錄因子的motif注釋
2.流程
2.1.得到genelist
library(RcisTarget)
library(Seurat)
library(SeuratData)
library(tidyverse)
#導(dǎo)入seruat對象
scRNA <- readRDS("scRNAsub.rds")
#獲取genelist
Idents(scRNA_harmony) <- "celltype"
geneLists <- FindAllMarkers(scRNA)
#篩選有意義的genelist
genelists <- geneLists %>% filter(cluster == 'endothelial-cells') %>% top_n(35,avg_log2FC)
head(geneLists)
genelists <- genelists$gene #這里也可以像周老師一樣弄成一個(gè)list
2.2使用RcisTarget內(nèi)置的轉(zhuǎn)錄因子的motif注釋數(shù)據(jù)
#這里有兩個(gè)motif注釋數(shù)據(jù)可以選 根據(jù)自己分析對象進(jìn)行選擇
# mouse:
# data(motifAnnotations_mgi)
# human:
data(motifAnnotations_hgnc) #選擇人的
2.3導(dǎo)入Gene-motif rankings:提供每個(gè)motif的所有基因的排名(~得分)
#這個(gè)數(shù)據(jù)是比較難的，這個(gè)數(shù)據(jù)要從一個(gè)數(shù)據(jù)庫上下載，
#很大，每一個(gè)文件大概1G，經(jīng)常會(huì)下載不完整，會(huì)導(dǎo)致讀入不行，會(huì)報(bào)錯(cuò)
#方法一，用R直接下載并讀取
featherURL <- "https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg19/refseq_r45/mc9nr/gene_based/hg19-tss-centered-10kb-7species.mc9nr.feather"
download.file(featherURL, destfile=basename(featherURL)) #這樣就會(huì)直接下載到當(dāng)前工作目錄下
#讀入
motifRankings <- importRankings("hg19-tss-centered-10kb-7species.mc9nr.feather") #這里假如數(shù)據(jù)下載不全讀取R會(huì)奔潰并restart，不是電腦的配置問題
#方法二，讀取我自己提前下好的數(shù)據(jù)，
motifRankings <- importRankings("cisTarget_databases/hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.feather") #我的數(shù)據(jù)存放在cisTarget_databases文件夾下，跟網(wǎng)上下載有啥不同呢 hg19與hg38就是數(shù)據(jù)的版本不同，10kb與500bp就是想要研究基因上下調(diào)控的范圍 10kb更大 所以讀入hg19或hg38都可以。
假如大家想要我的參考數(shù)據(jù)可以私信我
這樣3個(gè)讀入的文件已經(jīng)搞定了
2.4開始分析
可以直接使用cisTarget進(jìn)行一步法分析，cisTarget()運(yùn)行按順序執(zhí)行的步驟
(1)motif 富集分析,
(2)motif-TF注釋,和
(3)選擇的重要基因。
motifEnrichmentTable_wGenes <- cisTarget(genelists,
motifRankings,
motifAnnot=motifAnnotations_hgnc)#一步搞定分析
2.5,分析結(jié)果解讀
RcisTarget輸出：
RcisTarget的最終輸出的data.table（motifEnrichmentTable_wGenes）包含有關(guān)motif 富集的以下信息：
geneSet：基因集的名稱
motif：motif的ID
NES：基因集中基序的標(biāo)準(zhǔn)化富集得分
AUC：曲線下的面積（用于計(jì)算NES）
TFinDB：指示突出顯示的TF是包含在高置信度注釋（兩個(gè)星號）還是低置信度注釋（一個(gè)星號）中。
TF_highConf：根據(jù)'motifAnnot_highConfCat'注釋到基序的轉(zhuǎn)錄因子。 這個(gè)是最主要的 就可以從這里看出輸入的基因集富集到了哪一些轉(zhuǎn)錄因子，從結(jié)果中也可以看出一個(gè)TF可以與多個(gè)motif結(jié)合調(diào)節(jié)不同的基因。
TF_lowConf：根據(jù)'motifAnnot_lowConfCat'注釋到主題的轉(zhuǎn)錄因子。
erichedGenes：在給定motif上排名較高的基因。
nErnGenes：高度排名的基因數(shù)量
rankAtMax：在最大富集時(shí)的排名，用于確定富集的基因數(shù)。
2.6結(jié)果plot：plot前幾個(gè)富集到的motif和轉(zhuǎn)錄因子，蠻雞肋的
motifEnrichmentTable_wGenes_wLogo <- addLogo(motifEnrichmentTable_wGenes)
resultsSubset <- motifEnrichmentTable_wGenes_wLogo[1:10,] #
library(DT)
datatable(resultsSubset[,-c("enrichedGenes", "TF_lowConf"), with=FALSE],
escape = FALSE, # To show the logo
filter="top", options=list(pageLength=5))
關(guān)于motif的logo如何看：https://zhuanlan.zhihu.com/p/428416814
顯然一步分析法很難得到我們想要的結(jié)果，我們也可以將這些步驟作為單獨(dú)的命令分別運(yùn)行。對感興趣的一個(gè)輸出進(jìn)行分析，所以建議分步法分析。
3.分步法分析
3.1 motif 富集分析（就是計(jì)算輸入的基因集富集到哪些motif，并對富集的motif進(jìn)行打分）
# 1. Calculate AUC
motifs_AUC <- calcAUC(genelists, motifRankings)
3.2 motif-TF注釋 對富集到的motif進(jìn)行篩選，默認(rèn)篩選分?jǐn)?shù)為3的motifs，進(jìn)行TF富集。
# 2. Select significant motifs, add TF annotation & format as table
motifEnrichmentTable <- addMotifAnnotation(motifs_AUC,
nesThreshold=3,
motifAnnot=motifAnnotations_hgnc)
3.3 找出每個(gè)基序富集的基因 由于一種motif可能調(diào)控多個(gè)基因，看看每一個(gè)motif調(diào)控哪一些基因，那么就可以知道TF調(diào)控那些基因了，這就是說Rcis Target的分析是基于motif的
## 3. Identify significant genes for each motif
# (i.e. genes from the gene set in the top of the ranking)
# Note: Method 'iCisTarget' instead of 'aprox' is more accurate, but slower
motifEnrichmentTable_wGenes <- addSignificantGenes(motifEnrichmentTable,
geneSets=genelists,
rankings=motifRankings,
nCores=1,
method="aprox")
3.4 利用分步文件 構(gòu)建網(wǎng)絡(luò)圖（基因和motif）
signifMotifNames <- motifEnrichmentTable$motif[1:3]
incidenceMatrix <- getSignificantGenes(genelists,
motifRankings,
signifRankingNames=signifMotifNames,
plotCurve=TRUE, maxRank=5000,
genesFormat="incidMatrix",
method="aprox")$incidMatrix
library(reshape2)
edges <- melt(incidenceMatrix)
edges <- edges[which(edges[,3]==1),1:2]
colnames(edges) <- c("from","to")
library(visNetwork)
motifs <- unique(as.character(edges[,1]))
genes <- unique(as.character(edges[,2]))
nodes <- data.frame(id=c(motifs, genes),
label=c(motifs, genes),
title=c(motifs, genes), # tooltip
shape=c(rep("diamond", length(motifs)), rep("elypse", length(genes))),
color=c(rep("purple", length(motifs)), rep("skyblue", length(genes))))
visNetwork(nodes, edges) %>% visOptions(highlightNearest = TRUE,
nodesIdSelection = TRUE) #這種網(wǎng)絡(luò)圖都是自定義的，也可以畫TF與其互作基因的網(wǎng)絡(luò)圖
建議大家可以看看周老師的簡書：
http://m.itdecent.cn/u/06ae70ef31bc
4.References:
https://cloud.tencent.com/developer/article/1734803
http://m.itdecent.cn/p/6e1d71db4220
http://m.itdecent.cn/p/fd82f17db2d8
