从GEO的soft文件中提取探针序列信息
今天是生信星球陪你的第334天
大神一句话，菜鸟跑半年。我不是大神，但我可以缩短你走弯路的半年~
就像歌儿唱的那样，如果你不知道该往哪儿走，就留在这学点生信好不好~
这里有豆豆和花花的学习历程，从新手到进阶，生信路上有你有我！
问题
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL16956
我需要注释这个平台，想要下载它的表格，可是直接表格点下面的view full table会显示在网页中，有五万多行，网页只能显示两万多行，根本不完整。
思路
所以我只能下载文件了，隐约记得应该下载soft文件，可我不记得怎么读取了。
解决方案
1.搜索：
找到了第一个链接：
https://warwick.ac.uk/fac/sci/moac/people/students/peter_cock/r/geo/
哦 原来是这个样子，getGEO可以。但是我的soft好大，接近三百兆，很容易下载不完整。
2.从Rstudio自带的terminal运行linux，下载soft文件到本地
用wget来下载到本地再读取，刚好Rstudio可以运行linux哦，之前发过一次：https://www.jianshu.com/p/c16c7095e4b2
wget -c ftp://ftp.ncbi.nlm.nih.gov/geo/platforms/GPL16nnn/GPL16956/soft/GPL16956_family.soft.gz
gunzip GPL16956_family.soft.gz #可选
好像还是挺慢的。。。。等的我花儿都谢了。
3.读取数据并处理
下载的是压缩包，解压不解压都行的。都可以读取
刚才的帖子里有用信息就是读取方法：
library(GEOquery)
x=getGEO(filename = "GPL16956_family.soft")
读取的过程非常之慢,等着就好了。
step1:从GPL中提取有用的部分，只需要数据框的’ID’,’SEQUENCE’两列。
y=x@dataTable@table[,c('ID','SEQUENCE')]
head(y)
## ID
## 1 ASHGA5P000001
## 2 ASHGA5P000002
## 3 ASHGA5P000003
## 4 ASHGA5P000004
## 5 ASHGA5P000005
## 6 ASHGA5P000006
## SEQUENCE
## 1 GAGATAACTAGAACAGTGTCCCTCCCCTTTTATAACCTGTGTTTTTAGATTTCAAAAAGA
## 2 CTAACTGTCAAACAAAGATGGGCTAATTAAACGGATGCCAACTAATCGGAAAACTTCTGG
## 3 CCAGAAGCTCCAGCATGGTCCACCTACCAGAAGCTCCAGCATGATTTACCCATGAGTAGC
## 4 TAAGTACGGTAATGCAAGAAGGTGCAGCAGGCGTAGTGCATTACAGCAACACTGAAAAAC
## 5 TATGAGAGATGAAGAGATAAGGTCTAACTGTTACCTGGGCTGAACCCTTGGACTTCTAGG
## 6 CAGTAGCGAGTATTTACTAAGTACTTTCTATTTGCGAGGCCCTGATAAAAGTACTGTCCT
#断线好几次 保存一下
save(y,file = "16956.Rdata")
step2.重复值统计
几千万行肯定不对，进行去重统计，发现每个都重复了八百多次！不明白为什么。
library(tidyverse)
x=count(y,ID,sort = T)
head(x)
## # A tibble: 6 x 2
## ID n
## <chr> <int>
## 1 !Sample_channel_count = 1 875
## 2 !Sample_molecule_ch1 = total RNA 875
## 3 !Sample_organism_ch1 = Homo sapiens 875
## 4 !Sample_platform_id = GPL16956 875
## 5 !sample_table_begin 875
## 6 !sample_table_end 875
step3.去重复，两种方法
test1=distinct(y,ID,SEQUENCE)
test2=y[!duplicated(y),]
发现有一些叹号开头的行没有被跳过，有NA。所以去除NA。与GEO网页上的行数进行对比，发现一致，说明正确了。保存一下
test3= na.omit(test2)
nrow(test3)
## [1] 58944
save(test3,file = "GPL16956.Rdata")
write.csv(test3,file = "GPL16956.csv",row.names = F)
test4=read.csv("GPL16956.csv")
head(test4)
## ID
## 1 ASHGA5P000001
## 2 ASHGA5P000002
## 3 ASHGA5P000003
## 4 ASHGA5P000004
## 5 ASHGA5P000005
## 6 ASHGA5P000006
## SEQUENCE
## 1 GAGATAACTAGAACAGTGTCCCTCCCCTTTTATAACCTGTGTTTTTAGATTTCAAAAAGA
## 2 CTAACTGTCAAACAAAGATGGGCTAATTAAACGGATGCCAACTAATCGGAAAACTTCTGG
## 3 CCAGAAGCTCCAGCATGGTCCACCTACCAGAAGCTCCAGCATGATTTACCCATGAGTAGC
## 4 TAAGTACGGTAATGCAAGAAGGTGCAGCAGGCGTAGTGCATTACAGCAACACTGAAAAAC
## 5 TATGAGAGATGAAGAGATAAGGTCTAACTGTTACCTGGGCTGAACCCTTGGACTTCTAGG
## 6 CAGTAGCGAGTATTTACTAAGTACTTTCTATTTGCGAGGCCCTGATAAAAGTACTGTCCT
附上sessionInfo
sessionInfo()
## R version 3.5.1 (2018-07-02)
## Platform: x86_64-conda_cos6-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.5 LTS
##
## Matrix products: default
## BLAS/LAPACK: /pub/anaconda3/lib/R/lib/libRblas.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] knitr_1.22 forcats_0.4.0 stringr_1.4.0
## [4] dplyr_0.8.0.1 purrr_0.3.2 readr_1.3.1
## [7] tidyr_0.8.3 tibble_2.1.1 ggplot2_3.1.0
## [10] tidyverse_1.2.1 GEOquery_2.50.5 Biobase_2.42.0
## [13] BiocGenerics_0.28.0
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.1 cellranger_1.1.0 pillar_1.3.1 compiler_3.5.1
## [5] plyr_1.8.4 tools_3.5.1 evaluate_0.13 lubridate_1.7.4
## [9] jsonlite_1.6 nlme_3.1-137 gtable_0.3.0 lattice_0.20-38
## [13] pkgconfig_2.0.2 rlang_0.3.3 cli_1.1.0 rstudioapi_0.10
## [17] yaml_2.2.0 xfun_0.5 haven_2.1.0 withr_2.1.2
## [21] httr_1.4.0 xml2_1.2.0 generics_0.0.2 hms_0.4.2
## [25] grid_3.5.1 tidyselect_0.2.5 glue_1.3.1 R6_2.4.0
## [29] fansi_0.4.0 readxl_1.3.1 limma_3.38.3 modelr_0.1.4
## [33] magrittr_1.5 backports_1.1.3 scales_1.0.0 rvest_0.3.2
## [37] assertthat_0.2.1 colorspace_1.4-1 utf8_1.1.4 stringi_1.4.3
## [41] lazyeval_0.2.2 munsell_0.5.0 broom_0.5.1 crayon_1.3.4
简书：小洁忘了怎么分身
隔壁生信技能树公益视频合辑（学习顺序是linux，r，软件安装，geo，小技巧，ngs组学！）
国内看B站，教学视频链接：https://m.bilibili.com/space/338686099
国外看YouTube，教学视频链接：https://m.youtube.com/channel/UC67sImqK7V8tSWHMG8azIVA/playlists友情链接：
点击底部的“阅读原文”，获得更好的阅读体验哦😻
初学生信，很荣幸带你迈出第一步。
我们是生信星球，一个不拽术语、通俗易懂的生信知识平台。由于是2018年新号，竟然没有留言功能。需要帮助或提出意见请后台留言、联系微信或发送邮件到jieandze1314@gmail.com，每一条都会看到的哦~
请关注“恒诺新知”微信公众号，感谢“R语言“，”数据那些事儿“，”老俊俊的生信笔记“，”冷🈚️思“，“珞珈R”，“生信星球”的支持！
