生信马拉松 Day9-10 GEO数据分析笔记
原创
今天正式开始教画图了，具体的代码其实挺多地方讲到了，上课的好处就是可以听到很多细节和经验，是自己零散地找资料不能相比的，收获很多，感觉要全部吞下来还要再复习几遍
贴具体的代码好像挺重复的，只贴一些过去自学没有了解到的有趣的TIPS：
1、转录组数据不能直接套芯片的流程，错套的第一个表现就是dim结果是 行 为0
2、non-coding和普通array可以统一处理，但不能做富集分析，富集分析需要用编码蛋白做，或者先靶基因预测然后再做富集分析
3、normalize不是一定要加，不normalize会相对不够整齐，成钉子形状的往往是没取log
不会因为千八百异常表达的使整体数据发生明显差别，因此整体box明显不同的时候一定是异常样本
如何处理异常样本：
a.删除异常样本
b.exp=limma::normalizeBetweenArrays(exp)，一般齐和非常齐之间拉平
有负值：
a.取过log，有少量负值——正常（取log没加1，不影响使用）
b.没取过log，有负值——错误数据
c.有一半负值，中位数为0——做了标准化（给了一个不可逆的半成品，是作者不希望我们使用）
对2和3一般弃用，非要用的话就处理原始数据
4、boxplot范围落在0-4之间可能是运行了两遍log
5、Bioconductor的注释包，用find_anno(gpl_number)提示信息来找，包括全部注释R包，部分平台文件、部分自主注释（idmap里整理过的）
6、GPL页面的表格文件解析
没有gene symbol有ENTREZ、GB_ACC也可以，ORF一列也是，对gene_assignment中，在"//"里的第二栏里，带“--”说明不对应任何symbol，需要删去
7、一个探针对应多个基因（非特异性探针），难以解释，这些行直接去掉
8、对于lnkRNA不能直接用页面上的TargetID，尽量寻找symbol列
9、官网下载对应产品的注释表格：往往是付费的，不能看到
10、不是所有GPL都能找到注释
11、为什么会有很多探针空白
探针设置的前瞻性，20年前的探针可以测到2年前的基因，以前设计的探针可以测到未来的基因
12、富集分析找不到校正p<0.05的通路的解决方法
a.调整logFC、pvalue阈值，改动差异基因数量
b.不使用默认的padj(富集的)，而是使用原始p值，在文章里说清楚即可
c.换富集方法，GSEA也可以做kegg富集
d.调参数maxGSSize = 500，默认参数，表示500个以上的通路不考虑，可以调大
13、探针去除的三种方式
a.随机去重
b.取行和/行平均值最大
c.取探针平均值作为数据
不同方式结果有细微的差异
14、设置分组信息group后要转换为因子，需要设置levels=c("Normal","Disease")，使对照组在前，疾病组在后，因为在因子这个数据类型中，默认第一个位置为参考水平
15、没有Bioconductor注释包，也不能用idmap注释的各种情况处理方法：
# 1.特点是GEO官网能找到这个GPL，且GPL文件里明显有gene symbol行
# GPL570-55999
a = read.delim('GPL570-55999.txt',comment.char = "#")
#看了一下发现前几行是以#开头，表示注释，用comment.char参数去掉
colnames(a)
b = a[,c('ID','Gene.Symbol')]
# 把列名改为标准列名probe_id,symbol
colnames(b) = c('probe_id','symbol')
#空行和非特异性探针不要
k1 = b$symbol!="";table(k1)
library(stringr)
k2 = !str_detect(b$symbol,'///');table(k2)
ids = b[k1&k2,]
# 2.特点是有GEO官网能找到GPL，没有symbol，但是有ENTREZ ID
# GPL23270
rm(list = ls())
a = read.delim('GPL23270.txt',skip = 6)
b = a[,1:2]
#把ENTREZ 转 symbol
library(clusterProfiler)
library(org.Hs.eg.db)
s2e = bitr(b$ENTREZ_GENE_ID,
fromType = "ENTREZID",
toType = "SYMBOL",
OrgDb = org.Hs.eg.db)#人类,注意物种
#一部分基因没匹配上是正常的。<30%的失败都没事。
#其他物种http://bioconductor.org/packages/release/BiocViews.html#___OrgDb
#merge和inner_join都可以
ids = merge(b,s2e,by.x='ENTREZ_GENE_ID',by.y='ENTREZID')
ids =ids[,-1]
colnames(ids) = c('probe_id','symbol')
#GBACC在bitr里是 ACCNUM
# 3.有些基因名中存在“.1"".2"，也就是版本号，需要去掉“.1"".2"再做匹配
# 下面用的是编的示例数据
f = c('a.1','b.11','c.6')
str_split(f,'\\.',simplify = T)[,1] #"."号不可以直接写
# 4.GEO官网上有GPL，没有symbol，但是有gene_assignment
rm(list = ls())
a = data.table::fread("GPL17586.txt",data.table = F)
# 有报错，末尾行丢了就丢了
colnames(a)
b = a[,c(1,8)]
library(stringr)
tmp = str_split(b$gene_assignment," // ",simplify = T)
b$gene_assignment = tmp[,2]
colnames(b) = c('probe_id','symbol')
k1 = b$symbol!="---";table(k1)
k2 = b$symbol!="";table(k2)
ids = b[k1&k2,]16、筛选下载数据中的部分样本进行数据分析
library(stringr)
# 方法1：按照行号，能数的时候可以自己数行号
keep = c(1,2,5,6)
exp = exp[,keep]
pd = pd[keep,]
# 方法2：按照逻辑值，根据自己的数据特点编写:
# 可以是提取要保留的数据有的，也可以是不要的数据有的
# 无论如何设置，只要满足自己的需求即可
keep = !str_detect(pd$title,"and");table(keep)
exp = exp[,keep]
pd = pd[keep,]总结：
小洁老师讲的真的非常好，细腻又形象，准备学生信的一定要来学一波让你丝滑入门，自己学过点又整不大明白的推荐再来回炉重造，切实解决了我很多之前迷糊的痛点，而且把容易犯的错误特别着重强调了
这两天的内容零零散散在网上多少能找到一些，但是小洁老师分享的编程思维是比较宝贵的
之前自己瞎摸索的时候会把和一个项目有关的所有代码都写在一个脚本里，几百行的代码阅读困难，前面做修改，后面要重新加载数据，或者是只画个别图的时候也很麻烦，其实没有想过把步骤拆解开来，子步骤用易读的英文命名，这样甚至后续可以把每一步修改成函数，会更方便
另外就是这套代码最有趣的地方是，除了01_start_GEO和02_group_ids文件需要对下载的数据稍微调整，以及取出分组信息之外，后面所有的几个脚本不需要几乎不需要做修改，需要做修改的地方也会安排一个变量
logFC_t = 1
p_t = 0.05
#思考，如何使用padj而非p值
k1 = (deg$P.Value < p_t)&(deg$logFC < -logFC_t)
k2 = (deg$P.Value < p_t)&(deg$logFC > logFC_t)
deg = mutate(deg,change = ifelse(k1,"down",ifelse(k2,"up","stable")))
table(deg$change)这样一个小小的改变，可以调整的地方就变得非常明显，不容易遗漏了~
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。