📖 生信数据分析--分析流程，工具包等

1、关于WGCNA原理 1.1 建立共表达网络 在基因共表达网络中，节点node代表基因，边edge代表两个基因间共表达关系。 若一个基因同时与多个基因存在相关性，称为hub基因。 若一群基因存在高度互相相关，称为module。 基因共表达网络的展示形式一般为 n × n 邻接矩阵adjacency matrix(n个基因) ...

1、GSVA函数 1 2 3 4 5 6 7 # BiocManager::install("GSVA") library(GSVA) ?gsva gsva(expr = , #metrix格式表达矩阵(行--基因，列--样本) gset.idx.list = , #list格式基因集 method=c("gsva", "ssgsea", "zscore", "plage"), # defaul:gsva kcdf=c("Gaussian", "Poisson", "none")) # default:Gaussian 如上所示，gsva()函数有4个关键参数。 前两个参数分别提供表达矩阵与基因集，其中基因名需要为相同的格式。 method参数设置方法，默认为gsva kcdf参数设置分布模拟类型，默认为Gaussian，适合于array数据以及log标准化(CPM/RPKM/TPM等)的RNA-seq数据。Poisson分布适合RNA-seq的原始count矩阵。 2、示例数据 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 ## (1) 表达矩阵 # BiocManager::install("GSVAdata") library(GSVAdata) data(commonPickrellHuang) ### array矩阵 array_mt = exprs(huangArrayRMAnoBatchCommon_eset) dim(array_mt) # [1] 11508 36 array_mt[1:4,1:4] # NA19099 NA18523 NA19144 NA19137 # 8567 8.370526 8.544890 8.270395 8.585984 # 23139 7.536228 7.232781 7.469420 7.801986 # 7580 6.325102 6.405961 6.510249 6.285510 # 55619 9.134495 9.049938 9.704679 9.285202 ### count矩阵 count_mt = exprs(pickrellCountsArgonneCQNcommon_eset) dim(count_mt) # [1] 11508 36 count_mt[1:4,1:4] # NA19099 NA18523 NA19144 NA19137 # 8567 326 209 318 343 # 23139 255 169 245 361 # 7580 72 69 124 76 # 55619 487 590 678 540 ## (2) 基因集 library(tidyverse) library(msigdbr) hm_H = msigdbr(species = "Homo sapiens", category = "H", subcategory = NULL) %>% as.data.frame() %>% dplyr::select(gs_cat, gs_subcat, gs_name, entrez_gene) hm_H.list = split(hm_H$entrez_gene, hm_H$gs_name) str(head(hm_H.list)) # List of 6 # $ HALLMARK_ADIPOGENESIS : int [1:210] 19 11194 10449 33 34 35 47 50 51 112 ... # $ HALLMARK_ALLOGRAFT_REJECTION: int [1:335] 16 6059 10006 43 92 207 322 567 567 586 ... # $ HALLMARK_ANDROGEN_RESPONSE : int [1:102] 10257 11057 2181 87 9510 11047 9590 207 220 56172 ... # $ HALLMARK_ANGIOGENESIS : int [1:36] 350 351 894 1281 1290 6372 2260 11167 3685 182 ... # $ HALLMARK_APICAL_JUNCTION : int [1:231] 58 60 70 71 72 87 88 89 81 81 ... # $ HALLMARK_APICAL_SURFACE : int [1:46] 102 79602 79602 84632 9465 351 50617 5205 2683 672 ... 3、GSVA分析 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 ## (1) array + gsva res_gsva = gsva(array_mt, hm_H.list) class(res_gsva) # [1] "matrix" "array" dim(res_gsva) # [1] 50 36 res_gsva[1:4,1:4] # NA19099 NA18523 NA19144 NA19137 # HALLMARK_ADIPOGENESIS -0.21147644 0.08157479 0.004672524 -0.1656925 # HALLMARK_ALLOGRAFT_REJECTION -0.32620137 -0.44105423 -0.452996166 -0.4288909 # HALLMARK_ANDROGEN_RESPONSE 0.02056377 0.09123573 -0.135386552 -0.2198959 # HALLMARK_ANGIOGENESIS -0.02435541 -0.18282244 0.281235763 0.5114103 ## (2) count + ssgsea res_gsva = gsva(count_mt, hm_H.list, method = "ssgsea", kcdf = "Poisson") class(res_gsva) # [1] "matrix" "array" dim(res_gsva) # [1] 50 36 res_gsva[1:4,1:4] # NA19099 NA18523 NA19144 NA19137 # HALLMARK_ADIPOGENESIS 0.3288768 0.32899026 0.2893928 0.3297140 # HALLMARK_ALLOGRAFT_REJECTION 0.5018248 0.49183093 0.4672342 0.4531009 # HALLMARK_ANDROGEN_RESPONSE 0.3210905 0.28289473 0.2818011 0.2845855 # HALLMARK_ANGIOGENESIS 0.1794361 0.07705857 0.2076328 0.1276449 4、limma差异分析 以上面array + gsva的分析结果为例，随机模拟分组进行limma差异分析 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 ## 模拟分组 meta = data.frame(sample=colnames(array_mt), group=rep(c("G1","G2"), each=18)) head(meta) # sample group # 1 NA19099 G1 # 2 NA18523 G1 # 3 NA19144 G1 # 4 NA19137 G1 # 5 NA18861 G1 # 6 NA19116 G1 ## limma差异分析 library(limma) mod <- model.matrix(~ factor(meta$group)) colnames(mod) <- c("G1", "G2") fit <- lmFit(res_gsva, mod) fit <- eBayes(fit) tt <- topTable(fit, coef=2, n=Inf) head(tt) # logFC AveExpr t P.Value adj.P.Val B # HALLMARK_COMPLEMENT 0.08684907 -0.10731154 2.495909 0.01676849 0.3324198 -3.059104 # HALLMARK_HEME_METABOLISM 0.07942485 -0.04250442 2.238926 0.03076698 0.3324198 -3.486635 # HALLMARK_IL2_STAT5_SIGNALING 0.07839683 -0.09276742 2.193525 0.03411795 0.3324198 -3.558694 # HALLMARK_MITOTIC_SPINDLE -0.10721385 -0.07656961 -2.057416 0.04617883 0.3324198 -3.768076 # HALLMARK_IL6_JAK_STAT3_SIGNALING 0.09342746 -0.20104143 2.026107 0.04943145 0.3324198 -3.814790 # HALLMARK_INFLAMMATORY_RESPONSE 0.10211588 -0.10571870 1.996939 0.05263925 0.3324198 -3.857810

xCell xCell包 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 devtools::install_github('dviraran/xCell') library(xCell) data("xCell.data") ##查看支持的64种细胞类型，同下图 colnames(xCell.data$spill$K) ##预测函数的关键参数解释 ?xCellAnalysis() # expr = 交代表达矩阵； ##如果是array，不需要额外标准化；如果是RNAseq，需要TPM/FPKM/TPM。 ##对于基因ID格式需要是symbol格式。 # rnaseq = TRUE 数据是否为RNAseq数据，如果是芯片数据设置为FALSE # cell.types.use = NULL 提供一个字符串，说明想要预测64种细胞中的哪些细胞类型 # parallel.sz = 4 调用的线程数，默认为4 NOTE： ...

ClusterGVis包是中国药科大学Jun Zhang博士开发的系列可视化工具包之一，可以基因表达矩阵进行高级的热图可视化分析。如下根据其github以及微信教程简单整理一下自己感兴趣的用法。 ...

生存分析(survival analysis)的主要目的是发现与患者生存事件相关的指标因素，例如年龄性别、基因表达/突变等。如下学习相关基础知识及几种常见的生存分析方法。 ...

MuSiC(MUlti-Subject SIngle Cell deconvolution)是来自宾夕法尼亚大学Biostatistics, Epidemiology and Informatics系的Mingyao Li课题组于2019年发表于Nature Communication的一个工具R包，可根据单细胞转录组信息推测Bulk RNA-seq细胞组成。而后，该团队又于2022年在Briefing in bioinformatics发表了扩展版本MuSiC2，可以考虑更复杂的场景。 ...

RNA-seq数据比对流程主要分为三步（1）整理数据；（2）质控；（3）比对。其中每一步都涉及到若干软件的用法，如下简单整理基本的分析流程。 示例数据：GSE158623中6个样本的RNA-seq测序结果(human)，对应SRR12720999~SRR12721004 ...

1、HPA HPA：the Human Protein Atlas，由2003年来自瑞典的科研机构发起，旨在绘制综合性人类蛋白质图谱。 https://www.proteinatlas.org/ 蛋白质表达数据库，常见用途包括： （1）查看特定基因在不同组织、脑区，细胞类型，组织细胞类型，疾病（癌症），免疫细胞，肿瘤细胞系等表达情况。 （2）数据挖掘类文章常使用其进行比较基因在肿瘤部位与相应正常部位的蛋白水平表达差异。 2、ENCORI ENCORI，The Encyclopedia of RNA Interactomes，由中山大学生命科学学院屈良鹄团队开发，于2014年发表于Nucleic Acids Res。 https://starbase.sysu.edu.cn/index.php miRNA、lncRNA、RBP等多维相互作用网络，常见用途包括 miRNA/RBP的靶标查询(mRNA,lncRNA…) RNA interaction, ceRNA network TCGA肿瘤的差异表达，生存分析，相关性等 3、TCIA TCIA，The Cancer Immunome Database，由来自奥地利的因斯布鲁克大学医学院Zlatko Trajanoski团队开发，于2017年发表于Cell Reports https://tcia.at/home 20种solid cancer的免疫相关分析，例如 免疫基因表达、细胞浸润、肿瘤亚克隆等 亮点之一是提出Immunophenoscore指标用于预测免疫治疗响应 https://github.com/icbi-lab/Immunophenogram 样本临床信息也整理的较为完整 ...

1、TCGAbiolinks下载数据 使用TCGAbiolinks下载特定肿瘤类型的SNV数据 https://bioconductor.org/packages/release/bioc/vignettes/TCGAbiolinks/inst/doc/mutation.html 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 library(TCGAbiolinks) query <- GDCquery( project = "TCGA-CHOL", data.category = "Simple Nucleotide Variation", access = "open", legacy = FALSE, data.type = "Masked Somatic Mutation") GDCdownload(query) maf <- GDCprepare(query) dim(maf) # [1] 3764 141 ## (1) 因后续需要，修改Tumor_Sample_Barcode列 maf$long_Barcode = maf$Tumor_Sample_Barcode maf$Tumor_Sample_Barcode = substr(maf$Tumor_Sample_Barcode,1,12) length(unique(maf$Tumor_Sample_Barcode)) # 51 ## (2) 读取临床生存数据 clinical = readxl::read_xlsx("TCGA_Pan_Cancer_Clinical_Data_mmc1.xlsx") clinical_sle = clinical %>% dplyr::filter(type=="CHOL") %>% dplyr::select(bcr_patient_barcode, OS, OS.time, clinical_stage) %>% dplyr::rename(Tumor_Sample_Barcode=bcr_patient_barcode) dim(clinical_sle) # 45 2、maftools可视化 https://bioconductor.org/packages/release/bioc/vignettes/maftools/inst/doc/maftools.html https://www.jieandze1314.com/post/cnposts/237/ 1 2 3 4 5 6 7 8 9 library(maftools) maf_obj = read.maf(maf = maf, clinicalData = clinical_sle) #每个样本的突变情况统计 getSampleSummary(maf_obj) #每个基因的突变类型统计 getGeneSummary(maf_obj) 2.1概括图 1 2 3 plotmafSummary(maf = maf_obj, rmOutlier = TRUE, addStat = 'median', dashboard = TRUE, titvRaw = FALSE) 2.2 基因瀑布图 1 2 3 4 5 6 7 8 9 10 11 12 13 oncoplot(maf = maf_obj, top = 10) # 如下图 ## 添加临床注释 oncoplot(maf = maf_obj, top = 10, clinicalFeatures = c('clinical_stage',"OS"), draw_titv = TRUE) ## 选择特定基因集 set.seed(42) genes = sample(maf_obj@data$Hugo_Symbol,10) oncoplot(maf = maf_obj, genes = genes, clinicalFeatures = c('clinical_stage',"OS"), draw_titv = TRUE) 2.3 转换颠倒统计 1 2 3 4 5 # Transition 转换 ： 嘌呤(AG)或嘧啶(CT)内部之间转换 # Transversions 颠倒：嘌呤与嘧啶间互相转换 maf_obj.titv = titv(maf = maf_obj, plot = FALSE, useSyn = TRUE) plotTiTv(res = maf_obj.titv) 2.4 基因对突变统计 1 2 3 4 5 6 7 8 # green: co-occuring # yellow: mutually exclusive somaticInteractions(maf = maf_obj, top = 25, pvalue = c(0.05, 0.1)) set.seed(42) genes = sample(maf_obj@data$Hugo_Symbol,25) somaticInteractions(maf = maf_obj, genes = genes , pvalue = c(0.05, 0.1)) 2.5 生存分析 根据特定基因是否突变将病人分成WT与Mutant两组 1 2 3 4 5 6 7 8 9 10 11 mafSurvival(maf = maf_obj, genes = 'TP53', time = 'OS.time', Status = 'OS') # Group medianTime N # 1: Mutant 732 4 # 2: WT 650 41 ## 提取信息 # maf_obj@clinical.data %>% # dplyr::mutate(Group=ifelse(Tumor_Sample_Barcode %in% # subset(maf_obj@data, Hugo_Symbol=="TP53")$Barcode, # "Mutant","WT")) 2.6 基因对的生存相关性 1 2 3 4 5 6 7 8 9 10 prog_geneset = survGroup(maf = maf_obj, top = 200, geneSetSize = 2, time = "OS.time", Status = "OS", verbose = FALSE,minSamples = 3) prog_geneset # Gene_combination P_value hr WT Mutant # 1: PBRM1_PLXNA4 0.243 2.36e+00 42 3 # 2: PBRM1_PCLO 0.294 3.46e-01 42 3 # 3: PBRM1_TP53 0.320 3.64e-01 42 3 mafSurvGroup(maf = maf_obj, geneSet = c("PBRM1", "PLXNA4"), time = "OS.time", Status = "OS")

ConsensusClusterPlus包是肿瘤分型研究的常用工具，其于2010年发表于Bioinformatics。 Paper：https://academic.oup.com/bioinformatics/article/26/12/1572/281699 Tutorial：https://bioconductor.org/packages/release/bioc/vignettes/ConsensusClusterPlus/inst/doc/ConsensusClusterPlus.pdf 1 2 # BiocManager::install("ConsensusClusterPlus") library(ConsensusClusterPlus) 1、示例表达矩阵 行名为基因，列名为样本的基因表达矩阵；根据需要，进行标准化处理 基因集的选择是关键的一步，可根据统计学（高变基因）或者生物学（特定功能相关基因）进行选择 如下是参考教程的数据 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 library(ALL) data(ALL) d=exprs(ALL) dim(d) # [1] 12625 128 mads=apply(d,1,mad) d=d[rev(order(mads))[1:5000],] # 筛选高方差基因 dim(d) # [1] 5000 128 ## 标准化处理：每行（基因）减去每行的中位数 d = sweep(d,1, apply(d,1,median,na.rm=T)) d[1:4,1:4] # 01005 01010 03002 04006 # 36638_at 1.556121 0.9521271 -0.05018082 4.780378 # 39318_at 1.191353 2.5013225 -2.38793537 -1.199521 # 38514_at 1.020716 3.2785671 1.55949145 -3.345919 # 266_s_at 1.829260 0.3624327 1.54913247 -1.286294 2、样本亚型鉴定 实际就是一个函数即可：ConsensusClusterPlus() 其中涉及到较多参数的选择，具体如下 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 ## 默认参数 ConsensusClusterPlus(d, maxK=3, reps=10, pItem=0.8, pFeature=1, title="untitled_consensus_cluster", clusterAlg="hc", distance="pearson", plot=NULL, writeTable=FALSE, seed=42) ## 参数含义 # maxK：考虑的最大聚类数，建议取10~20 # reps：抽样次数，建议1000 # pItem与pFeature：分别表示对样本与基因的抽样比例 # title：图片或文件的保存路径 # clusterAlg：聚类方式 c("hc","pam","km") # distance: 距离计算方式 c("pearson","spearman","euclidean","binay","maximum","canberra","minkowski") # plot : 是否绘图，以及图形类型 c(NULL, "pdf", "png", "pngBMP") # writeTable: 是否保存文件 # seed：随机种子 results = ConsensusClusterPlus(d, maxK=10, reps=1000, pItem=0.8, pFeature=1, title="./tmp/", clusterAlg="hc", distance="pearson", plot="pdf", writeTable=FALSE, seed=42) R对象结果：list格式 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 ## 查看聚类数为2的结果信息 names(results[[3]]) # [1] "consensusMatrix" "consensusTree" "consensusClass" "ml" "clrs" table(results[[3]]$consensusClass) # 1 2 3 # 69 28 31 dim(results[[3]]$consensusMatrix) # [1] 128 128 results[[3]]$consensusMatrix[1:4,1:4] # [,1] [,2] [,3] [,4] # [1,] 1.0000000 0.3408360 0.7996870 0.3515249 # [2,] 0.3408360 1.0000000 0.1224806 1.0000000 # [3,] 0.7996870 0.1224806 1.0000000 0.1190108 # [4,] 0.3515249 1.0000000 0.1190108 1.0000000 ## consensus values 0 (never clustered together) ## consensus values 1 (always clustered together) 图形输出结果，主要两类图 （1）特定聚类结果的热图 ...