芯片分析里面有一个比较困难的地方就是探针的注释，其实我们之前的帖子已经把各种不同的情况解决的差不多该了。
除了使用R包注释，没有R包，我们可以自己从平台文件中获取探针和基因的对应关系。
来完成你的生信作业，这是最有诚意的GEO数据库教程
skr！GEO芯片数据的探针ID转换
有些GEO平台的探针转换比较麻烦
正则表达式是我们认识世界的哲学
还有一种情况比较特殊，我们也解决了
GEO芯片中的NM_，NR_开头的识别号如何转换成基因名称？
其实看完了以上的帖子，90%的GEO表达谱芯片的探针ID转换都没有问题，剩下的就是平台中只有序列没有其他可用信息的芯片了，这种芯片常常是非编码芯片，比如GPL21827芯片
他的平台信息是这样的，只有一段序列，没有其他可用的信息。
这时候就需要把这个序列比对一下，才行，blast人人都会用，只要把这个序列复制一个，去NCBI比对一下就行了。不过一个个操作肯定不行，下载NCBI-blast批量操作也很慢。
在唐医生的帮助下，我们使用seqmap这个软件解决了问题。
1.转录本信息文件下载
比对这个是事情，如果暂时理解不了，就用NCBI的blast去理解，后面到了高通量测序的时候也会讲。比对是要比对到基因组上去或者比对到已知的转录本上去，所以我们还缺个记录所有转录本信息的文件。这个文件保存在genecode数据库中。
点进去之后是这个样子的
再点进去是这样的，我们这个芯片是人的，所以选择人
点进去，选择最新的版本，30
选择转录本信息下载
2.GEO平台文件下载
浏览器输入平台名称
找到平台信息下载
点击进去下载到本地
因为这个文件有点大，自己下载有可能出现断断续续，或者下载不完全的情况，请下载完之后，检查文件大小。
3.比对软件下载
浏览器输入seqmap
点进去拉到最后，选择自己电脑对应的版本
4.处理平台文件变成fasta
这时候用Rstudio新建一个project，再把下载的三个文件放在project文件夹中，解压成三个文件。大概是这样的。
在检索seqmap的时候，第二个条目是其用法
点开阅读后，发现seqmap的用法比较简单，输入fasta格式的文件，输入要比对的参考基因组，再输入输出文件的文字，最后加上一个自定义选项即可。
其中，最困难的地方是fasta文件获取，这个只能从平台文件来转换。
gpl <- data.table::fread("GPL21827_family.soft",skip = "ID",data.table = F)
fread的skip参数可以跳过数字，也可以跳过字符。如果觉得不够优雅，可以用下面的方法
library(GEOquery)
gpl <- Table(getGEO(filename = "GPL21827_family.soft"))
缺点就是慢，尤其在当前数据下，这个方法很慢。
此时数据是这样的：
现在就选择两列，一列是探针名称，一列是序列
## 选取想要的两列，一列是ID， 一列是序列
gpl <- gpl[,c(1,4)]
过滤掉没有序列的行
library(dplyr)
gpl <- gpl %>%
filter(nchar(SEQUENCE)!=0)
保存为fasta格式文件,十分巧妙！
gp <- paste0('>',gpl$ID,'\n', gpl$SEQUENCE)
写成文本文件的时候就可以看出fasta的格式了
write.table(gp,'GPL.fasta', quote = F, row.names = F, col.names = F)
打开后数据是这个样子的
此时，我们要的三个文件都齐了
5.批量比对
这个可以在终端里面完成，建议安装git for windows，更加方便
在刚才那个文件夹下，右击选择git bash
在里面输入以下命令
./seqmap-1.0.12-windows.exe 0 ./GPL.fasta ./gencode.v30.transcripts.fa seqmap_results.txt //output_all_matches
这个命令由6部分组成
-
./seqmap-1.0.12-windows.exe是软件名称， -
0表示匹配容错率为0 -
./GPL.fasta是平台fasta文件 -
./gencode.v30.transcripts.fa是参考基因组 -
seqmap_results.txt是生成文件的名字 -
//output_all_matchest输出所有匹配结果，其他系统一个斜杠即可。
运行后是这个样子的：
1分钟不到就完成了。
6.提取比对结果
## 读入数据
probe2ID <- data.table::fread("seqmap_results.txt",data.table = F)
## 重要的结果在第一列
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)
结果是这个样子的。
大功告成！！
重点在这！！
这个技能用的地方不多，我把常见的非编码芯片平台信息已经注释好了，
文件储存为Rdata格式，使用的时候直接
load即可。
load(file = "GPL21827_probe_ID.Rdata")
微信联系我获取。
我的微信是guotosky, 加的时候注明一下原因，简单自我介绍一下即可。
除此之外，如果还有其他平台信息不在那5个之内的，也欢迎留言，我可以尝试去解决。
网友评论