前面已经给出了4个GEO芯片数据挖掘分析点，详见：
也是时候进阶到转录组测序了，它跟表达量芯片技术类似的，都是基于基因在不同样品的表达量矩阵做各式各样的分析，比如2022的文章：《Aberrant MYCN expression drives oncogenic hijacking of EZH2 as a transcriptional activator in peripheral T-cell lymphoma》，对应的数据集是：GSE198126
GSM5939531 m34-Thymus研究者们发现了在 Peripheral T cell lymphoma (PTCL)疾病里面有 recurrent overexpression of MYCN的现象，所以在做了一个小鼠模型：Inducible expression of MYCN in lymphoid cells in a mouse model ，做一个简单的转录组看看差异，如下所示：
简单的转录组看看差异 (A) Heat map showing the top differentially expressed genes in MYCN-driven TCL compared with normal CD4 T cells and depicting the transcription factor binding motifs of the top transcription factors identified by i-CisTarget in the up- and downregulated genes. NES = normalized enrichment score. (B) NESs for different hallmark gene sets in MYCN-driven TCL compared with normal CD4 T cells. (C) GSEA showing enrichment of an embryonic stem cell–like signature as described by Wong et al
值得注意是，作者不仅仅是做了差异分析和生物学功能数据库注释（gsea方法，癌症的50个hallmark基因集 ），而且还搭配了转录因子调控分析。
常规的差异分析 只需要读取作者在geo界面给出来了的表达量矩阵文件即可：
data"GSE198124_all.counts.csv.gz",F )# data$V1=str_split(data$V1,'_',simplify = T)[,2] # head(data) 2 ]1 :4 ,1 :4 ]1 ) > 1 ;table(keep_feature)1 :4 ,1 :4 ]1 :4 ,1 :4 ] library (AnnoProbe)'SYMBOL' ,'mouse' )'protein_coding' ,]in % ids$SYMBOL,]'CD4' , colnames(symbol_matrix) ),'control' ,'case' )'control' ,'case' )) 'symbol_matrix.Rdata' )然后使用转录组测序的金标准算法（DESeq2，edgeR，limma-voom）计算统计学显著的上下调基因即可，并且做如下所示的标准可视化：
统计学显著的上下调基因 精华当然是生物学功能数据库注释（gsea方法， 癌症的50个hallmark基因集 ），但是我一般来说会默认走一下kegg的gsea注释，代码如下所示：
rm(list = ls()) 'DEG_deseq2.Rdata' )'DEG_limma_voom.Rdata' )'DEG_edgeR.Rdata' ) # FoxO signaling pathway "log2FoldChange" , "padj" )]"logFC" , "adj.P.Val" )]'logFC' ,'P.Value' )library (org.Mm.eg.db)library (clusterProfiler)"ENTREZID" ,in % gene$SYMBOL,]# https://www.ncbi.nlm.nih.gov/gene/?term=SPARCL1 T )library (clusterProfiler) 'mmu' ,#按需替换 1000 ,10 ,0.99 ,FALSE )'org.Mm.eg.db' ,'ENTREZID' )#按需替换 ' - Mus musculus \\(house mouse\\)' ,'' ,kk@result$Description )'test' '_kegg.gsea.csv' )))'gsea_kk.Rdata' ))可以看到，虽然说没有选择癌症的50个hallmark基因集，但是结果是差不多的，上调的仍然是增殖相关通路，以及一个代谢通路。
kegg_top_gsea 常规转录组数据分析流程里面是没有转录因子调控分析的 我们上面的转录组测序的金标准算法（DESeq2，edgeR，limma-voom）能拿到统计学显著的上下调基因，但是没有这些基因集对应的转录因子，有很多计算工具可以实现，比如BART、i-cisTarget和Enrichr和cistrome的Lisa，其中这个文章选择了i-cisTarget。
但是按照我一直以来的代码，我会使用RcisTarget转录因子分析。它可以从一组基因中识别潜在的转录因子（Transcription Factors, TFs）调控网络，并通过 motif 分析 推测转录因子的作用机制。该包主要被应用于转录调控研究，特别是基因共表达网络或差异表达基因集的上游调控因素预测。需要看一些基本的文档：
SCENIC: single-cell regulatory network inference and clustering. Nat Methods. 2017 Nov;14(11):1083-1086. bioconductor: https://bioconductor.org/packages/release/bioc/vignettes/RcisTarget/inst/doc/RcisTarget_MainTutorial.R motif2annotaion: https://rdrr.io/github/aertslab/RcisTarget/man/motifAnnotations.html github: https://github.com/aertslab/RcisTarget RcisTarget database: https://resources.aertslab.org/cistarget/ 也可以直接使用我们的代码，只需要一句话代码cisTarget完成富集分析 ，如下所示：
library (RcisTarget)# Load gene sets to analyze. e.g.: 'RcisTarget' ), "hypoxiaGeneSet.txt" ), FALSE )[,1 ]## 重点在于拿到了上下调基因各自的列表哦： # Select motif database to use (i.e. organism and distance around TSS) # 这个RcisTarget包内置的motifAnnotations_hgnc是16万行 # 可以看到每个转录因子基因有多个motif，但是不到2000个转录因子 # 取决于R包“RcisTarget”的版本，确定是否运行下面的一行代码： if (F ){1 :4 ,1 :4 ]#motifRankings 1 :4 ,1 :4 ]function (x){in % colnames(motifRankings@rankings) ]# Motif enrichment analysis: # 输入数据可以是一个list，这里是geneLists，那就分别对list每一个元素，分别做TF的聚类 # 这个 cisTarget 函数等价于 下面的3个步骤： # 1. Motif enrichment analysis (calcAUC) # 2. Motif-TF annotation (addMotifAnnotation) # 3. Selection of significant genes (addSignificantGenes) 1 :4 ,1 :4 ]# 设置阈值，相应的TF对应富集到8个基因以上才被留下 8 '/' ,'motifEnrichmentTable_wGenes.csv' )) 上面的代码依赖于一个数据库文件，而且不同的物种不一样：
1.0G 11 19 2020 hg19-500bp-upstream-7species.mc9nr.feather
