GEO探针ID转换教程(代码+注释)!!! #################################################################################### ### 01_使用R包注释 rm(list = ls()) platformMap <- data.table::fread("platformMap.txt") ### platformMap.txt文档下载:链接:https://pan.baidu.com/s/1R71PENgN-i5aSYKmFSVdXw 密码:0gsx ### platformMap.txt这个文档整理的是:平台名称与其对应的R包是啥 ### 那么如何知道平台的名称? ### 打开搜到的NCBI GEO网站:https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE42872 ### Platforms后面显示:GPL6244,就是平台 ### 怎么找到GPL6244对应的R包呢? ### 可以点开platformMap.txt文档自己慢慢找,在“bioc_package”这一栏里的就是R包,注意后面要加上“.db” ### 也可以用R语言找: index <- "GPL6244" paste0(platformMap$bioc_package[grep(index,platformMap$gpl)],".db") ### grep即抓取,在platformMap的gpl列找index在第几行 ### 看不懂的代码可以拆开运行,就能理解啦 ### 复制此时Console里出现的那行R包名称! ### 安装R包: options(BioC_mirror="https://mirrors.tuna.tsinghua.edu.cn/bioconductor") if(!require("hugene10sttranscriptcluster.db")) BiocManager::install("hugene10sttranscriptcluster.db",update = F,ask = F) ### 获取探针: probe2symbol_df <- toTable(get("hugene10sttranscriptclusterSYMBOL")) # R包名称+SYMBOL probe2symbol_df1 <- toTable(get("hugene10sttranscriptclusterENTREZID")) # R包名称+ENTREZID probe2symbol_df2 <- toTable(get("hugene10sttranscriptclusterENSEMBL")) # R包名称+ENSEMBL #################################################################################### ### 02_解析SOFT文件 ### 如果平台没有R包可用,怎么办呢?不管怎样,下面这招都有用! ### 最核心的技能:解析SOFT文件!!! ### 下载NCBI GEO网站的“SOFT formatted family file(s)” ### 示例:https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE42872 rm(list = ls()) GPL_anno <-data.table::fread("GSE42872_family.soft",skip ="ID",data.table = F) ### 跳过SOFT文档里前面不需要的内容,直接到字符“ID” ### 提示:Console里出现“Warning”不用管,只有“Error”才要管 ### 那么怎么提取基因名呢? library(dplyr) library(tidyr) probe2symbol_df <- GPL_anno %>% dplyr::select(ID,gene_assignment) %>% # 选择GPL_anno的ID列和gene_assignment列 filter(gene_assignment != "---") %>% # 删掉gene_assignment列是"---"的,因为啥也没有 separate(gene_assignment,c("drop","symbol"),sep="//") %>% # 根据"//"分割,分割得到的前两个命名为"drop"和"symbol" dplyr::select(-drop) # 把“drop”去掉 ### 这样就完成啦! #-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-# ### 下面是个小练习,再重复一下separate的使用:主要就是如何分出对照组和处理组 ### GEO的样本名称太多而且排序不规则,你们都是手动分组的么?:https://mp.weixin.qq.com/s/m3hc4slyV9arw70kKPOClg ### 加载示例数据:链接:https://pan.baidu.com/s/1XMJ8KvxEPRGeRwCV3pZ2eg 密码:7om3 data <- data.table::fread("gsm_group.txt",data.table = F,header = F) names(data) <- c("GSM","group") ### 如果不会正则表达式,我们可以用separate: library(dplyr) library(tidyr) metadata <- data %>% separate(group,into = c("drop","group"),sep = "_blood_") %>% # 根据“_blood_”把group这一列分成两列,“_blood_”前面的叫“drop”、后面叫“group” select(-drop) %>% # 去掉drop列 separate(group,into = c("group","drop"),sep = "_") %>% # 类似地,根据“_”把group这一列分成两列 select(-drop) # 去掉drop列 ### 这样分组信息就出来啦 #-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-# ### 再来一个小练习:示例二:https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL4381 ### 网站显示:基因名是在小括号里的 rm(list = ls()) ### 数据下载方式一: # 从网站的“Download full table”处下载数据:不推荐!!! GPL_anno <-data.table::fread("GPL4381-4306.txt",skip ="ID",data.table = F) ### 数据下载方式二: # 下载GEO平台注释信息SOFT文件:“SOFT formatted family file(s)”:推荐!!! GPL_anno2 <- data.table::fread("GPL4381_family.soft",skip ="ID",data.table = F) ### 开始处理 library(dplyr) library(tidyr) probe2symbol_df <- GPL_anno2 %>% select(ID,GB_DEFINITION) %>% # 选择需要的两列 filter(GB_DEFINITION != "") %>% # 过滤掉“GB_DEFINITION”为空的行 separate(GB_DEFINITION,into=c("gene_symbol","drop"),sep = "\)") %>% # \代表转义,\)其实就是),因为)是元字符才加了\来给它转义 separate(gene_symbol,into = c("drop","gene_symbol"),sep = "\(") %>% # 去掉两边的括号 select(ID,gene_symbol) # 选择最终需要呈现的列 ### 不过这时候有些基因名很奇怪,比如像“from clone DKFZp547P055”或者“A”这样的,原因是:原先的“GB_DEFINITION”列在那一行基因名的括号前面还有括号 ### 所以这种方法不够完美,那怎么办呢? ### 使用正则表达式!!!!!! #################################################################################### ### 03_学习正则表达式 ### 学习正则表达式,一定要看这一篇:什么是正则表达式(Regular expression)? ### 30分钟的教程写了13年,被名字耽误的正则表达式! ### 链接:https://mp.weixin.qq.com/s/ZHY5bM1NUHUSJrqPuKJ-uw ### 在这里先介绍正则表达式中的三个括号:[]{}() ## 中括号[]表示选项,代表内部数据任意选择:比如[ATCG],表示A、T、C、G四个字符随意选择 ## 如果是数字也是一个意思:[1356],表示有1、3、5、6这四个选项 ## 如果嫌麻烦,可用-连接起始代表范围: # 比如[A-Z],代表大写的26个字母 # 比如[a-z],代表小写的26个字母 # 比如[0-9],代表从0到9的10个数字 # 还可以混写:比如[0-9a-zA-Z],代表数字和字母 ## 大括号{}代表重复次数:比如[ATCG]{3}表示从A、T、C、G四个字符中选择3次 ## 如果是两个数,用逗号隔开,代表范围:比如[0-9]{4,10},表示产生4位数到10位数都可以 ## 如果两个数中的第二个缺失:比如[0-9]{4,},代表4位数及以上 ## 圆括号()代表成组,即分组:实际上()的用处很多,但是分组是核心功能 rm(list = ls()) ### 示例:电话号码 ### 创建字符串 strings <- c( "apple", "219 733 8965", "329-293-8753", "Work: 579-499-7527; Home: 543.355.3679" ) ### 再创建一个可以匹配电话号码的模式: ## [2-9][0-9]{2}:三个数字 ## [- .]:-或者 或者. ## [0-9]{3}:三个数字 ## [0-9]{4}:四个数字 pattern <- "([2-9][0-9]{2})[- .]([0-9]{3})[- .]([0-9]{4})" ## ()里是我们想要的东西! #-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-# library(stringr) ### 先掌握以下7个函数: ######### ######### ### 0.如何学习 ## 使用str_view()来学习stringr ## 我们可以运行这个 str_view(strings, pattern) ## str_view_all(strings, pattern) # 不带all的话,找到匹配的就不找了;带all则查找所有 ### 右侧显示的就是这个模式匹配的情况 ######### ######### ### 1.查找 ## str_detect:从strings中检测能否和pattern匹配,返回的是逻辑值,有点像grepl str_detect(strings, pattern) ######### ######### ### 2.定位 ## str_locate:从strings中检测能否和pattern匹配,返回匹配的起止位置 ## 返回的是个数据框 str_locate(strings, pattern) ## str_locate_all(strings, pattern) ######### ######### ### 3.取回 ## 取回原始数据 str_subset(strings, pattern) ## 取回匹配到的数据 str_extract(strings, pattern) # extract提取 ## 取回所有匹配到的数据 # str_extract_all(strings, pattern) ## 取回所有匹配到的数据(简洁) # str_extract_all(strings, pattern, simplify = T) ######### ######### ### 4.模式取回:重要!!! ## 返回的是数据框,第一列是str_extract的数据 ## 后面依次是括号()中的内容,模式中有多少个(),就返回多少列 str_match(strings, pattern) ## str_match_all(strings, pattern) ######### ######### ### 5.替换 ## 从strings中能够精确匹配pattern的内容,并替换为“¥¥¥-¥¥¥¥-¥¥¥¥” str_replace(strings, pattern, "¥¥¥-¥¥¥¥-¥¥¥¥") ## str_replace_all(strings, pattern, "¥¥¥-¥¥¥¥-¥¥¥¥") ######### ######### ### 6.其他 ## 分割、粘贴、计数什么的,见官网:https://cran.r-project.org/web/packages/stringr/vignettes/stringr.html ### 这边再介绍三个表示匹配次数的符号:?+ * ## 问号?表示0次或者1次,因为这是一个生存或是毁灭的问题 ## 加号+表示1次或者多次 ## 星号*表示0次或者多次 ### 这样我们就得到了六种表示重复次数的方法: # 代码/语法 说明 # * 重复零次或更多次 # + 重复一次或更多次 # ? 重复零次或一次 # {n} 重复n次 # {n,} 重复n次或更多次 # {n,m} 重复n次到m次 ### 正则表达式如何匹配编码区域CDS? ## ATG([ATCG]{3})+(TGA|TAG|TAA) ## 解释如下: # 首先以起始密码子ATG开头,接着出现三联密码子,用[ATCG]{3},出现的次数是1次或者多次,所以整体给个加号+ # 现在就是这个样子的:ATG([ATCG]{3})+ # 接着要出现终止密码子,有三种分别是TGA、TGA、TAA,正则表达式中竖线|表示或的意思 # 那么(TGA|TAG|TAA)就代表,三个终止密码子任意选择一个,所以最终CDS区就表示成这个样子啦 #################################################################################### ### 04_正则表达式实战 ### 正则表达式实战运用! rm(list = ls()) ### str_match实战 #-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-# ### 例子一:提取分组 ## 加载数据 data <- data.table::fread("gsm_group.txt",data.table = F,header = F) names(data) <- c("GSM","group") ### 如果不会正则表达式,我们可以用separate,如前所述 library(dplyr) library(tidyr) metadata1 <- data %>% separate(group,into = c("drop","group"),sep = "_blood_") %>% select(-drop) %>% separate(group,into = c("group","drop"),sep = "_") %>% select(-drop) ### 而现在! ### str_match可以轻松做到! ## 同样加载数据 data <- data.table::fread("gsm_group.txt",data.table = F,header = F) names(data) <- c("GSM","group") ## 使用str_match library(stringr) group <- str_match(data$group,"whole_blood_(.*)_.*") # 我们想要匹配的是data$group这一列 # 正则表达式中:用点号.代表任何字符,用*号代表出现0次或者多次,所以.*代表任意字符串 # 那么类似whole_blood_control_1就可以写成whole_blood_.*_.* # ()里是我们想要的东西 metadata2 <- data.frame(GSM=data$GSM,group=group[,2]) #-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-# #-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-# #-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-# ### 例子二:抽取基因!!! ### 例子二:抽取基因!!! ### 例子二:抽取基因!!! ### 正则表达式是我们认识世界的哲学:https://mp.weixin.qq.com/s/DlioHHXQd-W-96tXLWrQvA rm(list = ls()) ### 打开NCBI GEO网站:https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL4381 ### 下载SOFT formatted family file(s)文档: GPL_anno <-data.table::fread("GPL4381-4306.txt",skip ="ID",data.table = F) ### 如何把基因读取出来呢? ### 先观察总结一下GB_DEFINITION列的语句特征 GPL_anno$GB_DEFINITION ### 挑个例子展示一下: string = "Homo sapiens intraflagellar transport 80 homolog (Chlamydomonas) (IFT80), mRNA" ### 总结正则表达式: pattern = "\(([A-Za-z0-9]*)\)," library(stringr) str_view(string,pattern) str_match(string,pattern) str_match(string,pattern)[,2] ### 这里就验证了正则表达式的正确性! ### 前方高能! ### 前方高能! ### 前方高能! library(dplyr) library(tidyr) probeGenesymbol <- GPL_anno %>% ## 增加一个新的列gene_symbol,这个列获取了str_match括号中的内容 mutate(gene_symbol = str_match(.$GB_DEFINITION,pattern)[,2]) %>% # .代表当前数据,即管道符传下来的数据,选取GB_DEFINITION这一列 ## 过滤掉gene_symbol为空的行 filter(gene_symbol != "") %>% ## 选出需要的列 select(ID,gene_symbol,GB_DEFINITION) ### 从现在开始,我们从SOFT文件里解析想要的东西就几乎没有问题啦! #################################################################################### ### 05_使用bioMart来注释 ### GEO芯片中的NM_、NR_开头的识别号如何转换成基因名称? ### GEO芯片中的NM_、NR_开头的识别号如何转换成基因名称? ### GEO芯片中的NM_、NR_开头的识别号如何转换成基因名称? ### 参考帖:https://mp.weixin.qq.com/s/FdCcliMCYj4Yb4grzIQMaA 这边有讲各种开头对应的分子类型,亦可见下: ## RefSeq accession numbers and molecule types: # Accession prefix Molecule type Comment # AC_ Genomic Complete genomic molecule, usually alternate assembly # NC_ Genomic Complete genomic molecule, usually reference assembly # NG_ Genomic Incomplete genomic region # NT_ Genomic Contig or scaffold, clone-based or WGS # NW_ Genomic Contig or scaffold, primarily WGS # NZ_ Genomic Complete genomes and unfinished WGS data # NM_ mRNA Protein-coding transcripts (usually curated) # NR_ RNA Non-protein-coding transcripts # XM_ mRNA Predicted model protein-coding transcript # XR_ RNA Predicted model non-protein-coding transcript # AP_ Protein Annotated on AC_ alternate assembly # NP_ Protein Associated with an NM_ or NC_ accession # YP_ Protein Annotated on genomic molecules without an instantiated transcript record # XP_ Protein Predicted model, associated with an XM_ accession # WP_ Protein Non-redundant across multiple strains and species ### 下载SOFT formatted family file(s)文件:https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL16686 rm(list = ls()) GPL_anno <-data.table::fread("GPL16686_family.soft",skip ="ID",data.table = F) ### 点开看看,发现GB_ACC这一列的NM_、NR_似曾相识,实际上它们是RefSeq的基因识别号,别的啥信息也没有 ### 使用Ensemble旗下的biomaRt包处理 options(BioC_mirror="https://mirrors.tuna.tsinghua.edu.cn/bioconductor") if(!require("biomaRt")) BiocManager::install("biomaRt",update = F,ask = F) library(biomaRt) ### 测试一下 ensembl = useMart(biomart="ensembl", dataset="hsapiens_gene_ensembl") ## 这里dataset是人的!不是人要改!就像从购物车(Mart)里取出使用一样! ## 这里运行需要一点时间,大概一两分钟 searchDatasets(mart = ensembl, pattern = "hsapiens") refseqids = c("NM_005359","NM_000546") # 测试两个数据 ### 使用getBM函数来转换 ipro = getBM(attributes=c("refseq_mrna","hgnc_symbol"), # 返回什么样的数据类型 filters="refseq_mrna", # 现在我们的数据属于什么类型,下面紧接着会讲到 values=refseqids, # 我们测试的数据 mart=ensembl) # ensembl是上面刚获取的 ### 点开ipro,测试ID转换已经完成 ### 上面attributes参数输入的是要返回的列,可以是多个 ### 通过listAttributes函数可以获取可输入的attributes!!!!!! listAttributes <- listAttributes(mart = ensembl) index <- grep("refseq",listAttributes$name) # 查看refseq相关的attributes有哪些,grep抓取 listAttributes[index,] # 列举出来 ### 结果显示和refseq相关的有这些: # name description page # 86 refseq_mrna RefSeq mRNA ID feature_page # 87 refseq_mrna_predicted RefSeq mRNA predicted ID feature_page # 88 refseq_ncrna RefSeq ncRNA ID feature_page # 89 refseq_ncrna_predicted RefSeq ncRNA predicted ID feature_page # 90 refseq_peptide RefSeq peptide ID feature_page # 91 refseq_peptide_predicted RefSeq peptide predicted ID feature_page ### 再测试一下 ### 我们尝试把ncRNA的基因名提取出来: refseqids = c("NR_036215","NR_026927","NR_015434") ipro = getBM(attributes=c("refseq_ncrna","hgnc_symbol"), filters="refseq_ncrna", values=refseqids, mart=ensembl) ### 完成 #-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-# ### 实战!!! ### 实战!!! ### 实战!!! ### 实战!!! ### 实战!!! ### 实战!!! ### 下载SOFT formatted family file(s)文档: GPL_anno <- data.table::fread("GPL16686_family.soft",skip ="ID",data.table = F) ### 1.先提取NM开头的(编码RNA): library(dplyr) refseqids <- GPL_anno %>% filter(str_detect(GB_ACC, "NM.*")) %>% # str_detect函数的使用,*代表后面随便是什么都可以:在管道符前阻断#一下,查看意义 pull(GB_ACC) # 新的小函数pull ### head(refseqids)查看 ### 批量转换:代码和之前一样 ipro_NM = getBM(attributes=c("refseq_mrna","hgnc_symbol"), filters="refseq_mrna", values=refseqids, mart=ensembl) ### 2.再提取NR开头的(非编码RNA): library(dplyr) refseqids <- GPL_anno %>% filter(str_detect(GB_ACC, "NR.*")) %>% # 换成了NR.* pull(GB_ACC) ipro_NR = getBM(attributes=c("refseq_ncrna","hgnc_symbol"), filters="refseq_ncrna", # 换成了refseq_ncrna values=refseqids, mart=ensembl) ### 合并信息:把编码RNA和非编码RNA合并 ### 先处理一下,防止接下来不清楚哪个是哪个 meta_NM <- cbind(ipro_NM,meta="coding") # 给ipro_NM加一列:列名meta,内容是coding colnames(meta_NM) <- c("ID","symbol","meta") # 把列名重新命名一下,变成ID、symbol、meta meta_NR <- cbind(ipro_NR,meta="noncoding") # 给ipro_NR加一列:列名meta,内容是noncoding colnames(meta_NR) <- c("ID","symbol","meta") # 把列名重新命名一下,变成ID、symbol、meta ### 列名已经一致,用rbind把它们结合在一起 probeGenesymbol <- rbind(meta_NM,meta_NR) ### 点开probeGenesymbol,发现非编码RNA有一些是空白的 ### 去掉空白 probeGenesymbol <- probeGenesymbol %>% filter(symbol!="") ### 这样探针ID转换就完成啦!!! ### 强力推荐帖:如何让基因名称在多个数据库间随意转换? ### 链接:https://mp.weixin.qq.com/s/wsiceQmNVveoggiqeDSlmQ #################################################################################### ### 06_非编码RNA的芯片注释 ### 其实做到这里,90%的GEO表达谱芯片的探针ID转换都没有问题了 ### 剩下的就是平台中只有序列没有其它可用信息的芯片了,这种芯片常常是非编码芯片! ### 参考帖:GEO芯片分析的倒数第2个关卡被没有了 链接:https://mp.weixin.qq.com/s/X8rUnEasKy3Dk-EoUAvC2A ### 1.转录本信息文件下载 ## 我们要将序列比对到基因组或转录本上去,所以我们缺个记录所有转录本信息的文件 ## 这个文件保存在GENCODE数据库中,百度gencode,点击进入GENCODE官网:https://www.gencodegenes.org/ ## 点击How to access data ## 根据芯片是人或鼠,选择human或mouse: # The GENCODE release files can be found in this website or directly in our human and mouse FTP sites. ## 比如是人的,点击human ## 选择版本30(不用太新),点击release_30/(小鼠则选版本20) ## 下载转录本信息:gencode.v30.transcripts.fa.gz # human(30版本)ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_30/gencode.v30.transcripts.fa.gz # mouse(20版本)ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_mouse/release_M20/gencode.vM20.transcripts.fa.gz ### 2.GEO平台文件下载 ## Google输入平台名称,比如:GPL15314(不要用百度),点进去 # 网址这种格式:https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL15314 ## 下载这个平台的SOFT formatted family file(s)文档,也就是平台信息 # 这个文件有点大,自己下载可能出现断断续续或者下载不完全的情况,必须下载完成之后检查文件大小 ### 3.比对软件下载 ## 百度搜索SeqMap,点击进入SeqMap - Short Sequence Mapping Tool:http://www-personal.umich.edu/~jianghui/seqmap/ ## 点击进去拉到最后,下载全平台版本:1.0.13 source (for all platform) # SeqMap的用法:http://www-personal.umich.edu/~jianghui/seqmap/Docs.txt # 点开阅读后,发现SeqMap的用法比较简单:在第5部分会介绍 ### 4.处理平台文件变成fasta ## 用Rstudio新建一个R project:File - New Project - New Directory - New Project - 命名一下,地址选desktop,这样就出现了一个文件夹 ## 再把前面下载的三个文件放在project文件夹中并解压,包括: # gencode.v30.transcripts.fa.gz,并解压 # GPL15314.soft,即GEO平台的SOFT formatted family file(s)文件 # seqmap-1.0.13-src.zip,并解压 ## 新建一个R Script,命名为SeqMap,复制以下内容: #-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-# rm(list = ls()) ### 读取平台信息 gpl <- data.table::fread("GPL15314.soft",skip = "ID",data.table = F) ## 这个看起来不大对劲,可能skip = "ID"的读取方式不合适,需要换一种 gpl <- data.table::fread("GPL15314.soft",skip = 489,data.table = F) ## 这下对了,但是这个第489行是如何确定的呢?通过打开SOFT文件确定 ## 如何打开SOFT文件?Windows可以用Notepad++,Mac可以用Sublime Text ## 489这个数字是如何确定的? ## 它是“ID SPOT_ID CONTROL_TYPE ENSEMBL_ID ACCESSION_STRING CHROMOSOMAL_LOCATION SPOT_ID SEQUENCE”所在行的上一行的行数!!! ## skip = 489的意思是忽略我们所需数据列名上方的所有内容 library(dplyr) ### 下面开始选取想要的两列:一列是ID,一列是Sequence gpl <- gpl[,c(1,8)] # 因为ID和Sequence分别在第1和第8列 colnames(gpl) <- c("ID","SEQUENCE") # 换一下列名,供代码复用 ### 过滤掉没有序列的行以及阴性对照行: library(dplyr) gpl <- gpl %>% filter(nchar(SEQUENCE)!=0) ### 转换为fasta格式文件,十分巧妙! ### 转换为fasta格式文件,十分巧妙! ### 转换为fasta格式文件,十分巧妙! gp <- paste0('>',gpl$ID,'n',gpl$SEQUENCE) ### 保存fasta到文件夹: write.table(gp,'GPL.fasta',quote = F,row.names = F,col.names = F) ### 这是关键步骤,现在我们需要的3个文件(转录本信息文件、比对软件以及fasta文件)就全部准备好啦!!! #-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-# ### 5.批量比对 ### MAC可以在终端里面完成,Windows建议安装git for windows ### 下面介绍Mac版: ## 解压已下载的SeqMap全平台版本:1.0.13 source (for all platform) ## 查阅其中的Makefile文件,发现指令为: # seqmap:g++ -O3 -m64 -o seqmap match.cpp # clean:rm -f seqmap ## 右键点击解压得到的文件夹seqmap-1.0.13-src,点击“服务”-“新建位于文件夹位置的终端窗口” # 在终端里粘贴g++ -O3 -m64 -o seqmap match.cpp,点击Enter,得到一个叫seqmap的程序 # seqmap百度云链接:https://pan.baidu.com/s/1b5egNWEWZmUrC4fbM_UQtA 密码:5qsj ## 把得到的seqmap程序复制到我们的R project文件夹 ## 右键点击我们的R project文件夹,点击“服务”-“新建位于文件夹位置的终端窗口” # 在终端里粘贴: # ./seqmap 0 ./GPL.fasta ./gencode.v30.transcripts.fa seqmap_results.txt /output_all_matches ## 即可得到比对结果seqmap_results.txt文档 # 解释一下上面的命令,这个命令由6部分组成: # ./seqmap是软件名称 # 0表示匹配容错率为0 # ./GPL.fasta是平台fasta文件 # ./gencode.v30.transcripts.fa是参考基因组 # seqmap_results.txt是生成文件的名字 # /output_all_matchest输出所有匹配结果 #-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-# ### 读入比对结果 rm(list = ls()) probe2ID <- data.table::fread("seqmap_results.txt",data.table = F) ### 重要的数据在probe2ID的第一列 library(tidyr) library(dplyr) probe2ID <- probe2ID %>% select(probe_id,trans_id) %>% separate(trans_id,into = c("Ensembl", "drop1","drop2","drop3", "trans_Symble","gene_Symble","drop4","trans_biotype"),sep = "\|") %>% select(probe_id,Ensembl,trans_Symble,gene_Symble,trans_biotype) # 根据\|切割,不要的drop掉 ### 保存一下 save(probe2ID,file = "GPL15314_probe_ID.Rdata") ### 提取lncRNA: ncRNA <- c("processed_transcript", "non_coding", "3prime_overlapping_ncRNA", "antisense", "lincRNA", "sense_intronic", "sense_overlapping", "macro_lncRNA", "bidirectional_promoter_lncRNA") probe_lncRNA <- probe2ID %>% filter(trans_biotype %in% ncRNA) # trans_biotype这一列要在ncRNA范围里 ### 提取mRNA: probe_mRNA <- probe2ID %>% filter(trans_biotype =="protein_coding") #-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-# ### 这边我们已经把常见的非编码芯片平台信息注释好了! ### 文件储存为Rdata格式,使用的时候直接load即可:例如load(file = "GPL21827_probe_ID.Rdata") ### 平台包括:GPL13825、GPL15314、GPL16956、GPL19612、GPL21827!!! ### 非编码注释文件下载链接:https://pan.baidu.com/s/1z5Co_z6c7iV0EkVyUu2dzQ 密码:xevw #################################################################################### ### 07_环状RNA的注释 ### 环状RNA数据库:http://www.circbase.org/ ### circRNA ID(人)长这样:hsa_circ_0018207(举例) ### circRNA的序列是无法比对到基因组上的,而且实验上敲减、过表达也很难处理,不推荐做 ### 这里只提供两个关键的ID转换信息文档: ## 1.网址:https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL26925 # Title:Agilent-084217 CapitalBio Technology Human CircRNA Array v2 # 读取方式: GPL_anno1 <-data.table::fread("GPL26925_family.soft",skip ="ID",data.table = F) ## 2.与circbase的转换文件,来自于网络: # 读取方式: GPL_anno2 <-data.table::fread("CircRNA_ID_Arraystar_Human_CircRNA.csv",data.table = F) ### 有这两个文件就足够使用了 ### 下载链接:https://pan.baidu.com/s/1eivRg5zL5gQeI49GMnEibA 密码:w0hp #################################################################################### ### 08_完成探针ID转换 ### 最后一步!!! ### 最后一步!!! ### 最后一步!!! ### 整合表达量数据与探针ID转换文件 ### 示例:GPL6244平台的数据 rm(list = ls()) load(file = "exprSet_readGSE.Rdata") load(file = "probe2symbol_df.Rdata") ### 合并用的是merge或inner_join函数!两个函数功能是一样的! ### 这里最重要的一步是:保证两个数据中有名称一样的列,就是我们合并的轴 ### 换句话说:merge必须要有一样的列,而不能靠一样的行名 ### 经典方法:把行名变成第一列!:cbind代表按列合并 exprSet <- cbind(probe_id = rownames(exprSet),exprSet) # 第1列叫probe_id ### 那么两个数据里名称一样的列就是probe_id ### 1.方法一:使用merge(第1个数据,第2个数据,by = "共同的列") exprSet1 <- merge(probe2symbol_df,exprSet,by = "probe_id") ## 过河拆桥:去掉第一列 exprSet1 <- exprSet1[,-1] ### 2.方法二:使用inner_join(第1个数据,第2个数据,by = "共同的列") library(dplyr) exprSet2 <- inner_join(probe2symbol_df,exprSet,by = "probe_id") ## 过河拆桥:去掉第一列 exprSet2 <- exprSet2[,-1] ### 这一招,搞定所有探针ID转换!!! ### 这一招,搞定所有探针ID转换!!! ### 这一招,搞定所有探针ID转换!!! ### 因为存在多个探针对应一个基因的情况,所以要加上探针去重的步骤! ### 探针去重时几个探针识别同一基因,选取保留最大值(最大值说明检测到过这么大的表达量,其它的检测量较小可能是因为那种探针对应的东西降解了) library(tibble) exprSet3 <- exprSet %>% # 合并探针的信息 inner_join(probe2symbol_df,by="probe_id") %>% # 去掉多余信息 select(-probe_id) %>% # 重新排列 select(symbol,everything()) %>% # 求出平均数(这边的.代表上面传入的数据) # .[,-1]表示去掉输入数据的第一列,然后求行的平均值 mutate(rowMean = rowMeans(.[,-1])) %>% # 把表达量的平均值按从大到小排序 arrange(desc(rowMean)) %>% # 去重,symbol留下第一个 distinct(symbol,.keep_all = T) %>% # 反向选择去除rowMean这一列 select(-rowMean) ### 参考资料: ## GEO芯片中多个探针对应一个基因,是求平均值还是保留最大值?:https://mp.weixin.qq.com/s/sjmXP7z8jjUpQGD5mZNP5w ## 如果GEO中一个探针对应多个基因,如何把这个探针全部删掉?:https://mp.weixin.qq.com/s/iXlMiKY8RBLEkvCctKa6hA ## 凡是重复的,全部删掉,一个都不留!:https://mp.weixin.qq.com/s/OvZmMbnxqcC-5KABthb0AQ ## group_by和summrise连用后,分组计算就很方便:https://mp.weixin.qq.com/s/w6NdxsmetSitF19-vbbYXA ### 完结撒花 ### 完结撒花 ### 完结撒花
暂无讨论,说说你的看法吧