寫(xiě)在前面：今天學(xué)習(xí)了Rcis Target這個(gè)包闷哆，每一次在學(xué)習(xí)一個(gè)包之前噪生，苦惱和痛苦是不知道這個(gè)R包用來(lái)do for what苹享？why to do院领？今天看了http://www.reibang.com/p/6e1d71db4220周老師的簡(jiǎn)書(shū)弛矛，感覺(jué)講的就是我自己，下面我將用我的生物學(xué)知識(shí)來(lái)講關(guān)于Rcis Target的理解比然，假如有錯(cuò)的地方歡迎大家指正與批評(píng)Ｍ羲摺！谈秫！
前言：我們?cè)趩渭?xì)胞分析的時(shí)候，通過(guò)差異基因來(lái)判斷哪些clusters是什么類(lèi)型的細(xì)胞鱼鼓，所以每一種細(xì)胞類(lèi)型是一些基因（即基因集）選擇表達(dá)的結(jié)果拟烫，而基因的表達(dá)是受轉(zhuǎn)錄調(diào)控的（分子生物學(xué)知識(shí)），單從轉(zhuǎn)錄因子（TF）這一個(gè)方面來(lái)說(shuō)迄本，不同的轉(zhuǎn)錄因子調(diào)控不同的基因硕淑，進(jìn)而使細(xì)胞表現(xiàn)出不同的狀態(tài)或類(lèi)型。當(dāng)然轉(zhuǎn)錄調(diào)控不只有TF調(diào)控 還有很多其他的如：順式作用元件嘉赎。
1.do for what：用B細(xì)胞為例置媳，我們可以對(duì)B細(xì)胞中高度相關(guān)的基因集（即：讓這些細(xì)胞為B細(xì)胞的基因）進(jìn)行Rcis Target分析，就可以知道哪一些TF在B細(xì)胞是富集或特有的公条，籠統(tǒng)的說(shuō)Rcis Target分析就是用來(lái)預(yù)測(cè)基因集的轉(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)行比較寥袭，通過(guò)篩選過(guò)濾找出二者具有統(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對(duì)象
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ù)自己分析對(duì)象進(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ù)庫(kù)上下載，
#很大佳遣，每一個(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窗声，不是電腦的配置問(wèn)題
#方法二，讀取我自己提前下好的數(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開(kāi)始分析
可以直接使用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：基因集的名稱(chēng)
motif：motif的ID
NES：基因集中基序的標(biāo)準(zhǔn)化富集得分
AUC：曲線下的面積（用于計(jì)算NES）
TFinDB：指示突出顯示的TF是包含在高置信度注釋?zhuān)▋蓚€(gè)星號(hào)）還是低置信度注釋?zhuān)ㄒ粋€(gè)星號(hào)）中。
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)行凑耻。對(duì)感興趣的一個(gè)輸出進(jìn)行分析，所以建議分步法分析柠贤。
3.分步法分析
3.1 motif 富集分析（就是計(jì)算輸入的基因集富集到哪些motif香浩，并對(duì)富集的motif進(jìn)行打分）
# 1. Calculate AUC
motifs_AUC <- calcAUC(genelists, motifRankings)
3.2 motif-TF注釋 對(duì)富集到的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)控那些基因了镜盯，這就是說(shuō)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ò)圖都是自定義的，也可以畫(huà)TF與其互作基因的網(wǎng)絡(luò)圖
建議大家可以看看周老師的簡(jiǎn)書(shū)：
http://www.reibang.com/u/06ae70ef31bc
4.References:
https://cloud.tencent.com/developer/article/1734803
http://www.reibang.com/p/6e1d71db4220
http://www.reibang.com/p/fd82f17db2d8
