转录组

在对转录组数据分析时，分组比较的差异分析前提时获得测序表达矩阵。 根据测序技术分为两种，对应的分析R包如下所示。 RNAseq的原始count表达矩阵 DESeq2 edgeR limma 芯片array标准化后的表达矩阵 limma 1、示例count表达矩阵 如下使用TCGAbiolinks包下载共9个癌症样本、35个正常样本的表达矩阵；其中涉及到8对配对样本（tumor与normal来自同一病人） ...

TCGAbiolinks包是一站式分析TCGA数据的R包工具，它集成了TCGA数据下载、分析、可视化的全部流程。此次系列笔记主要跟着 TCGAbiolinks帮助文档重新学习下TCGA数据挖掘流程。 ...

1、计算公式 sample 1 sample 2 …….. sample k Gene 1 10 12 30 Gene 2 20 25 60 …… … … … … Gene n 0 0 … 1 对于(n*k)表达矩阵中k个样本的n个基因的count表达数据。 在任意样本（列）中，基因i的长度为l，count表达值为q；则基因i的FPKM与TPM标准化方式计算如下 ...

挖掘GEO数据时，主要一方面是下载GEO的测序数据（包括基因芯片array与RNAseq两类）的表达矩阵。同时会涉及到一些细节问题，例如array芯片ID转换、样本meta信息等。 ...

1、准备conda环境与软件 1 2 3 4 5 6 7 8 # 准备download环境 conda create -n download conda activate download # 安装软件 conda install -c hcc aspera-cli conda install -c bioconda sra-tools conda install -c conda-forge pigz 以SRR13911909为例 2、aspera下载 https://www.ebi.ac.uk/ena/browser/view/SRR13911909 2.1 下载链接 1 2 3 4 5 # era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/SRR139/009/SRR13911909/SRR13911909_1.fastq.gz # era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/SRR139/009/SRR13911909/SRR13911909_2.fastq.gz # 观察上述链接规律后，可以自动生成下载链接 era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/${id:0:6}/00${id:0-1}/${id}/${id}_1.fastq.gz era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/${id:0:6}/00${id:0-1}/${id}/${id}_2.fastq.gz 2.2 aspera下载 1 2 3 4 5 6 7 8 9 10 11 12 13 id=SRR13911909 ascp -QT -l 300m -P33001 \ -i ~/miniconda3/envs/download/etc/asperaweb_id_dsa.openssh \ era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/${id:0:6}/00${id:0-1}/${id}/${id}_1.fastq.gz . ascp -QT -l 300m -P33001 \ -i ~/miniconda3/envs/download/etc/asperaweb_id_dsa.openssh \ era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/${id:0:6}/00${id:0-1}/${id}/${id}_2.fastq.gz . #最后一个点表示下载文件的储存路径 #如果是其它环境，将download替换为对应环境名即可 #如果是base环境： ~/miniconda3/etc/asperaweb_id_dsa.openssh 使用aspera可以达到百兆的速度，建议首先尝试。但最近试了几次，容易报错，不稳定(报错内容如下)；有时候可以。 ...

以前通路富集分析直接使用clusterprofiler包，阅读文献发现GSEA分析及可视化较多使用Broad团队研发的工具，现简要学习其(window版本)使用方法。 ...

1、背景知识 （1）两种富集分析 基于超几何检验的ORA(over representation analysis)富集分析 ① 假设对转录组分组测序的10000个基因表达数据进行差异分析； ② 按照规定cutoff阈值(logFC/Pvalue)筛选得到100个差异表达基因； ...

0、示例数据 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 library(clusterProfiler) library(org.Hs.eg.db) data(geneList, package="DOSE") str(geneList) # Named num [1:12495] 4.57 4.51 4.42 4.14 3.88 ... # - attr(*, "names")= chr [1:12495] "4312" "8318" "10874" "55143" ... gene_ids<-AnnotationDbi::select(org.Hs.eg.db, keys=as.character(names(gene_list)), columns="SYMBOL", #目标格式 keytype="ENTREZID") #目前的格式 # 'select()' returned 1:1 mapping between keys and columns ## ORA分析 gene_deg = names(geneList)[abs(geneList) > 2] ego <- enrichGO(gene = gene_deg, OrgDb = org.Hs.eg.db, ont = "CC", # "BP","MF","ALL" pAdjustMethod = "BH", pvalueCutoff = 0.01, qvalueCutoff = 0.05, readable = TRUE) ## GSEA分析 gseago <- gseGO(geneList = geneList, OrgDb = org.Hs.eg.db, ont = "CC", minGSSize = 100, maxGSSize = 500, pvalueCutoff = 0.05, verbose = FALSE) gseago = setReadable(gseago, OrgDb = org.Hs.eg.db) 1、ggplot2 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 ## ORA library(tidyverse) ego_df = ego@result %>% head(10) %>% dplyr::mutate(logp = -log10(p.adjust)) %>% dplyr::arrange(logp) %>% dplyr::mutate(Label=str_wrap(gsub("_"," ",Description), width = 50)) %>% dplyr::mutate(Label = factor(Label, levels=Label)) %>% tibble::remove_rownames() %>% dplyr::select(Label, logp) head(ego_df) # Label logp # 1 microtubule associated complex 5.419514 # 2 spindle midzone 5.681085 # 3 spindle microtubule 7.341138 ggplot(ego_df, aes(x=Label, y=logp)) + geom_bar(stat="identity", fill="#1b9e77") + xlab("Pathway names") + ylab("-log10(P.adj)") + scale_y_continuous(expand=c(0,0)) + coord_flip() + xlab("GO Pathway") + ylab("Adjusted P value(-log10)") + geom_text(aes(label = Label, y=0.2), hjust = 0, size=5) + theme_classic() + theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(), axis.title.y = element_text(size = 25), axis.title.x = element_text(size = 20), axis.text.x = element_text(size = 16)) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 ## GSEA library(tidyverse) gseago_df = gseago@result %>% head(20) %>% dplyr::mutate(logp = -log10(p.adjust)) %>% dplyr::arrange(logp) %>% dplyr::mutate(Label=str_wrap(gsub("_"," ",Description), width = 50)) %>% dplyr::mutate(Label = factor(Label, levels=Label)) %>% tibble::remove_rownames() %>% dplyr::select(Label, logp, NES) head(gseago_df) # Label logp NES # 1 polymeric cytoskeletal fiber 3.755730 1.487284 # 2 chromosome, telomeric region 4.311816 1.907511 # 3 microtubule associated complex 4.432327 1.944046 ggplot(gseago_df, aes(NES, fct_reorder(Label, NES), fill=qvalue)) + geom_bar(stat='identity') + scale_fill_continuous(low='red', high='blue', guide=guide_colorbar(reverse=TRUE)) + theme_minimal() + ylab(NULL) 2、GOplot 官方教程：https://wencke.github.io/ 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 # install.packages('GOplot') library(GOplot) # (1) 准备数据 ## 富集通路结果 res_enrich = ego@result %>% dplyr::select(ID, Description, geneID,p.adjust) %>% dplyr::mutate(Category = "GO", .before=1) %>% dplyr::mutate(geneID = gsub("/",", ", geneID)) %>% dplyr::rename("Term"="Description", "Genes"="geneID", "adj_pval"="p.adjust") %>% dplyr::select(Category,ID,Term,Genes,adj_pval) %>% tibble::remove_rownames() t(res_enrich[1,]) # 1 # Category "GO" # ID "GO:0005819" # Term "spindle" # Genes "CDCA8, CDC20, KIF23, CENPE, ASPM, DLGAP5, SKA1, NUSAP1, TPX2, TACC3, NEK2, CDK1, MAD2L1, KIF18A, BIRC5, KIF11, TRAT1, TTK, AURKB, PRC1, KIFC1, KIF18B, KIF20A, AURKA, CCNB1, KIF4A" # adj_pval "6.339976e-11" ## 如上，ID为optional，Category为required ## (2) 基因差异倍数 res_logfc = data.frame(ID=gene_ids$SYMBOL, logFC=geneList) head(res_logfc) # ID logFC # 1 MMP1 4.572613 # 2 CDC45 4.514594 # 3 NMU 4.418218 ## 如无相关数据，可随便设置为0 ## (3) 构建绘图对象 circ <- circle_dat(res_enrich, res_logfc) head(circ) Circular visualization ...

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