GEO生信Bioinformatics

GEO芯片分析的倒数第2个关卡没有了！

2019-06-20  本文已影响18人  9d760c7ce737

芯片分析里面有一个比较困难的地方就是探针的注释，其实我们之前的帖子已经把各种不同的情况解决的差不多该了。
除了使用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部分组成

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个之内的，也欢迎留言，我可以尝试去解决。

上一篇 下一篇

猜你喜欢

热点阅读