GEO探针ID转换教程（代码+注释）_转载，非原创
2021-02-12 本文已影响0人
whykm
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_matchlibrary(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### 完结撒花### 完结撒花### 完结撒花
