GEO踩坑记录：subscript out of bounds

我是伪放假的小洁，回学校走毕业的各项流程，老板跟我说，我有一百多兆的代码要传授给你！我需要把踩坑过程和解决过程记录下来，emmm然后就有了今天的推文😄如果你也想学习，请找到生信菜鸟团数据挖掘系列：

GEO数据挖掘-第一期-胶质母细胞瘤（GBM）
跟着这个文章跑代码，突然报了一个错
（不需要示例数据了，看一看型推文😄）

AssayData <- newAssayData[ID2gene$V1,]
Error in newAssayData[ID2gene$V1, ] : subscript out of bounds

猜测一：索引范围大于向量范围

subscript out of bounds这个报错我见过，有可能是因为索引范围大于向量范围，例如一个向量只有三列，你索引取第四列，就报这个错。所以我来验证一下：
先看看这两个向量：

head(ID2gene$V1)
#[1] "1007_s_at" "1007_s_at" "1053_at"   "117_at"    "121_at"   
#[6] "1255_g_at"

head(rownames(newAssayData))
#[1] "1007_s_at" "1053_at"   "117_at"    "121_at"    "1255_g_at"
#[6] "1294_at"  

length(ID2gene$V1)
#[1] 53249
nrow(newAssayData)
#[1] 54613

看起来数据类型是一样的，并且索引没有比向量长，那么问题出在哪里？

all(rownames(ID2gene$V1 %in% newAssayData ))
[1] TRUE

猜测二：行名重复

众所周知，行名不能重复，会不会是因为重复值

length(unique(ID2gene$V1))
#[1] 6792

这个结果让我很震惊，五万多行，去重复完了剩六千？？？？
不太相信，换个函数再确认

library(dplyr)
nrow(count(ID2gene,V1))
#[1] 6792

啊？？？哦
那么用count看一下每个探针重复多少次。

head(count(ID2gene,V1))
# A tibble: 6 x 2
#  V1            n
#  <chr>     <int>
#1 NA        45857
#2 1007_s_at     2
#3 1053_at       1
#4 117_at        1
#5 121_at        1
#6 1255_g_at     1

原来是有缺失值？而且四万多个？
查看缺失值的另一种方法：

table(is.na(ID2gene$V1))
#FALSE  TRUE 
# 7392 45857 

谜底揭开，就是因为NA。在顺着代码往上看，发现上游的featureData数据框长这样

而这个数据框是这样来的：

featureData <- exprSet@featureData@data

再网上找，就会找到下载文件那一步，所以是数据下载或者读取不完整的锅巴。

GEOquery还有一个需要注意的问题，下载不完整，需要删除原来下载的文件，再重新下载。因为它可以自行判断，本地有原始文件，就从本地读取，本地没有才会下载。

需要重新下载，就删掉他俩

这个让我想联想到一件事，数据框下载和读取，是这样现有了行数列数，确定了框架再一行行往里填，如果中断，就会留下一大片NA。画图和写文件都是这样，如果错了，就会生成一个0K的空文件，存在但是打不开。