GEO探针ID转换教程(代码+注释)

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

### 完结撒花
### 完结撒花
### 完结撒花
生物信息学

批次矫正教程(全代码+注释)

2020-8-13 21:33:14

生物信息学

生存分析曲线作图代码

2020-8-13 21:34:12

加入Q群
0 条回复 A文章作者 M管理员
    暂无讨论,说说你的看法吧
个人中心
今日签到
有新私信 私信列表
搜索