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

挖掘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可以达到百兆的速度，建议首先尝试。但最近试了几次，容易报错，不稳定(报错内容如下)；有时候可以。 ...

refgenie：参考基因组（阿拉丁）商店 http://refgenie.databio.org/ Here we provide a web interface and a RESTful API to access genome assets for popular reference genome assemblies. 该平台由位于弗吉尼亚大学公共卫生基因组学中心的计算生物学和生物信息学研究小组(Sheffield lab of computational biology)建立。上次修改/更新是2021年11月。 ...

1 2 3 4 5 if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("STRINGdb") library(STRINGdb) 1、定义要使用的STRING版本、物种，以及PPI阈值分数 1 2 3 4 string_db <- STRINGdb$new(version="11", species=9606, score_threshold=200, input_directory="") 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 data(diff_exp_example1) genes = diff_exp_example1[1:50,] head(genes) # pvalue logFC gene # 1 0.0001018 3.333461 VSTM2L # 2 0.0001392 3.822383 TBC1D2 # 3 0.0001720 3.306056 LENG9 # 4 0.0001739 3.024605 TMEM27 # 5 0.0001990 3.854414 LOC100506014 # 6 0.0002393 3.082052 TSPAN1 ###基因名匹配protein ID #第一个参数是data.frame; 第二个参数是基因所在列的列名 genes_mapped <- string_db$map(genes, "gene" ) #Warning: we couldn't map to STRING 30% of your identifiers head(genes_mapped) # gene pvalue logFC STRING_id # 1 VSTM2L 0.0001018 3.333461 9606.ENSP00000362560 # 2 TBC1D2 0.0001392 3.822383 9606.ENSP00000481721 # 3 LENG9 0.0001720 3.306056 9606.ENSP00000479355 # 4 TMEM27 0.0001739 3.024605 9606.ENSP00000369699 # 40 LOC100506014 0.0001990 3.854414 <NA> # 5 TSPAN1 0.0002393 3.082052 9606.ENSP00000361072 #string_db$plot_network(genes_mapped$STRING_id) 3、下载这些基因间的互作关系 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 ppi = string_db$get_interactions(genes_mapped$STRING_id) %>% distinct() ppi = ppi %>% dplyr::left_join(genes_mapped[,c(1,4)], by=c('from'='STRING_id')) %>% #列标序号根据具体情况而定 dplyr::rename(Gene1=gene) %>% ##列名根据具体情况而定 dplyr::left_join(genes_mapped[,c(1,4)], by=c('to'='STRING_id')) %>% dplyr::rename(Gene2=gene) %>% dplyr::select(Gene1, Gene2, combined_score) head(ppi) # Gene1 Gene2 combined_score # 1 C3 TYROBP 240 # 2 ABCA12 GRHL3 308 # 3 FAM189A1 TM4SF20 400 # 4 ABCA12 NIPAL4 824 # 5 GRHL3 NIPAL4 275 # 6 GRHL3 IGDCC4 238

1、GO 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 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 GO_data <- clusterProfiler:::get_GO_data("org.Hs.eg.db", "ALL", "SYMBOL") names(GO_data) # [1] "PATHID2NAME" "EXTID2PATHID" "GO2ONT" "PATHID2EXTID" ###（1）GO term的组成基因 class(GO_data$PATHID2EXTID) #[1] "list" GO_data$PATHID2EXTID[1] # $`GO:0000002` # [1] "PARP1" "SLC25A4" "DNA2" "TYMP" "LIG3" "MEF2A" # [7] "MPV17" "OPA1" "TOP3A" "TP53" "LONP1" "AKT3" # [13] "PPARGC1A" "POLG2" "SLC25A36" "PIF1" "SESN2" "SLC25A33" # [19] "MGME1" "PRIMPOL" "STOX1" ###（2）基因所涉及的通路 class(GO_data$EXTID2PATHID) #[1] "list" GO_data$EXTID2PATHID[1] # $A1BG # [1] "GO:0001775" "GO:0002252" "GO:0002263" "GO:0002274" "GO:0002275" "GO:0002283" "GO:0002366" # [8] "GO:0002376" "GO:0002443" "GO:0002444" "GO:0002446" "GO:0002576" "GO:0003674" "GO:0005575" # ... ###（3）GO term的名字 class(GO_data$PATHID2NAME) #[1] "character" GO_data$PATHID2NAME[1] # GO:0000001 # "mitochondrion inheritance" ###（4）GO term的类别 class(GO_data$GO2ONT) #[1] "character" GO_data$GO2ONT[1] # GO:0000002 # "BP" table(GO_data$GO2ONT) # BP CC MF # 16013 1981 4755 library(tidyverse) go_name = reshape2::melt(GO_data$PATHID2NAME) %>% rownames_to_column("ID") %>% dplyr::rename("Name"="value") go_type = reshape2::melt(GO_data$GO2ONT) %>% rownames_to_column("ID") %>% dplyr::rename("Type"="value") go_info = inner_join(go_name, go_type) %>% dplyr::mutate(GSEA=toupper(gsub(" ","_",paste0("GO",Type," ",Name)))) dim(go_info) head(go_info) table(rownames(brca_enrich_kegg) %in% go_info$GSEA) # ID Name Type # 1 GO:0000002 mitochondrial genome maintenance BP # 2 GO:0000003 reproduction BP # 3 GO:0000009 alpha-1,6-mannosyltransferase activity MF # 4 GO:0000010 trans-hexaprenyltranstransferase activity MF # 5 GO:0000012 single strand break repair BP # 6 GO:0000014 single-stranded DNA endodeoxyribonuclease activity MF # GSEA # 1 GOBP_MITOCHONDRIAL_GENOME_MAINTENANCE # 2 GOBP_REPRODUCTION # 3 GOMF_ALPHA-1,6-MANNOSYLTRANSFERASE_ACTIVITY # 4 GOMF_TRANS-HEXAPRENYLTRANSTRANSFERASE_ACTIVITY # 5 GOBP_SINGLE_STRAND_BREAK_REPAIR # 6 GOMF_SINGLE-STRANDED_DNA_ENDODEOXYRIBONUCLEASE_ACTIVITY 此外 GO.db包也提供了除组成基因以外的GO注释信息 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 library(GO.db) keytypes(GO.db) # [1] "DEFINITION" "GOID" "ONTOLOGY" "TERM" goids = keys(GO.db, keytype = "GOID")[1:3] # [1] "GO:0000001" "GO:0000002" "GO:0000003" goids_anno = AnnotationDbi::select(GO.db, keys = goids, columns = c("TERM","ONTOLOGY","DEFINITION"), #其中DEFINITION为term的详细描述 keytype="GOID") #所有的BP term的GO id goBP = select(GO.db, keys = "BP", columns = c("GOID"), keytype="ONTOLOGY") 2、KEGG 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 hsa_kegg <- clusterProfiler::download_KEGG("hsa") names(hsa_kegg) # [1] "KEGGPATHID2EXTID" "KEGGPATHID2NAME" ### KEGG id与name head(z) # from to # 1 hsa00010 Glycolysis / Gluconeogenesis # 2 hsa00020 Citrate cycle (TCA cycle) # 3 hsa00030 Pentose phosphate pathway ### KEGG id的组成基因 head(hsa_kegg$KEGGPATHID2EXTID) # from to # 1 hsa00010 10327 # 2 hsa00010 124 # 3 hsa00010 125

一、基因组大小 （1）人类基因组主要由细胞核的23对染色体组成(核基因组)，还包括线粒体中的小DNA分子(线粒体基因组)。 （2）单倍体基因组大概有30亿个碱基对组成，具体到每个染色体的碱基对长度与基因数量如下所示(参照UCSC的hg38)。 ...

1、不同基因ID转换 1.1 org.Hs.eg.db包 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 library(dplyr) library(org.Hs.eg.db) keytypes(org.Hs.eg.db) # [1] "ACCNUM" "ALIAS" "ENSEMBL" "ENSEMBLPROT" "ENSEMBLTRANS" "ENTREZID" # [7] "ENZYME" "EVIDENCE" "EVIDENCEALL" "GENENAME" "GENETYPE" "GO" # [13] "GOALL" "IPI" "MAP" "OMIM" "ONTOLOGY" "ONTOLOGYALL" # [19] "PATH" "PFAM" "PMID" "PROSITE" "REFSEQ" "SYMBOL" # [25] "UCSCKG" "UNIPROT" gene_symbol=c("RHO","CALM1","MEG3","GNGT1","SAG","RPGRIP1","TRPM1","PCP2","PCP4","AP1B1") gene_ids<-AnnotationDbi::select(org.Hs.eg.db, keys=as.character(gene_symbol), columns=c("ENSEMBL","ENTREZID"), #目标格式 keytype="SYMBOL") #目前的格式 gene_ids ##去重 gene_ids %>% dplyr::distinct(ENTREZID, .keep_all = T) # SYMBOL ENSEMBL ENTREZID # 1 RHO ENSG00000163914 6010 # 2 CALM1 ENSG00000198668 801 # 3 MEG3 ENSG00000214548 55384 # 4 GNGT1 ENSG00000127928 2792 # 5 SAG ENSG00000130561 6295 # 6 RPGRIP1 ENSG00000092200 57096 # 7 TRPM1 ENSG00000134160 4308 # 8 PCP2 ENSG00000174788 126006 # 9 PCP4 ENSG00000183036 5121 # 10 AP1B1 ENSG00000100280 162 1.2 biomaRt包 1 2 3 4 5 6 7 8 9 10 11 12 library("biomaRt") ensembl = useMart("ensembl",dataset="hsapiens_gene_ensembl") attributes = listAttributes(ensembl) attributes[1:5,] # library(httr) # httr::set_config(config(ssl_verifypeer = 0L)) gene_symbol=c("RHO","CALM1","MEG3","GNGT1","SAG","RPGRIP1","TRPM1","PCP2","PCP4","AP1B1") gene_ids2 <- getBM(filters= "hgnc_symbol", attributes= c("hgnc_symbol","ensembl_gene_id","entrezgene_id"), values = gene_symbol, mart= ensembl) gene_ids2 2、鼠源基因转为人类基因ID 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 musGenes <- c("Hmmr", "Tlx3", "Cpeb4") ## 方式1：直接大小写转换 toupper(musGenes) # [1] "HMMR" "TLX3" "CPEB4" ## 方式2：通过biomaRt包（不稳定） require("biomaRt") # library(httr) # httr::set_config(config(ssl_verifypeer = 0L)) human = useMart("ensembl", dataset = "hsapiens_gene_ensembl",host = "dec2021.archive.ensembl.org") mouse = useMart("ensembl", dataset = "mmusculus_gene_ensembl",host = "dec2021.archive.ensembl.org") genes = getLDS(attributes = c("mgi_symbol"), filters = "mgi_symbol", values = musGenes, mart = mouse, attributesL = c("hgnc_symbol"), martL = human, uniqueRows=T) ## 方式3：MGI 数据库 # https://support.bioconductor.org/p/129636/ library(dplyr) mouse_human_genes = read.csv("http://www.informatics.jax.org/downloads/reports/HOM_MouseHumanSequence.rpt",sep="\t") convert_mouse_to_human <- function(gene_list){ output = c() for(gene in gene_list){ class_key = (mouse_human_genes %>% filter(Symbol == gene & Common.Organism.Name=="mouse, laboratory"))[['DB.Class.Key']] if(!identical(class_key, integer(0)) ){ human_genes = (mouse_human_genes %>% filter(DB.Class.Key == class_key & Common.Organism.Name=="human"))[,"Symbol"] for(human_gene in human_genes){ output = append(output,human_gene) } } } return (output) } convert_mouse_to_human(musGenes) # 1] "HMMR" "TLX3" "CPEB4" 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 # # https://github.com/lishensuo/utils # # library("biomaRt") # # library(httr) # # httr::set_config(config(ssl_verifypeer = 0L)) # human = useMart("ensembl", dataset = "hsapiens_gene_ensembl",host = "dec2021.archive.ensembl.org") # mouse = useMart("ensembl", dataset = "mmusculus_gene_ensembl",host = "dec2021.archive.ensembl.org") # # # https://www.gencodegenes.org/mouse/ # dat = data.table::fread("gencode.vM33.basic.annotation.gtf.gz") # dat = subset(dat, V3 == "gene") # dat_sub = dat[,"V9"] %>% # separate(V9, into = c("gene_id","gene_type","gene_name","mgi_id","havana_gene"), sep = "; ") # dat_sub$gene_name2 = gsub('gencode.vM33.basic.annotation.gtf.gz "','',dat_sub$gene_name) # dat_sub$gene_name2 = gsub('"','',dat_sub$gene_name2) # # genes = getLDS(attributes = c("mgi_symbol"), filters = "mgi_symbol", # values = dat_sub$gene_name2, # mart = mouse, # attributesL = c("hgnc_symbol"), # martL = human, uniqueRows=T) # write.csv(genes, file = "mgi2hgnc_biomart.csv", row.names = F, quote = F) # head(genes) 3、蛋白质与基因ID转换 https://www.uniprot.org/uploadlists/ ...

以前通路富集分析直接使用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 ...