表达芯片数据分析2

> find_anno(gpl_number) #打出找注释的代码
`library(hgu133plus2.db);ids <- toTable(hgu133plus2SYMBOL)` and `ids <- AnnoProbe::idmap('GPL570')` are both avaliable
if you get error by idmap, please try different `type` parameters
> library(hgu133plus2.db);ids <- toTable(hgu133plus2SYMBOL)
########注释方法############
#方法1 BioconductorR包(最常用，已全部收入find_anno里面，不用看啦)
if(F){
  gpl_number #看看编号是多少
  #http://www.bio-info-trainee.com/1399.html #在这里搜索，找到对应的R包
  library(hgu133plus2.db)
  ls("package:hgu133plus2.db") #列出R包里都有啥
  ids <- toTable(hgu133plus2SYMBOL) #把R包里的注释表格变成数据框
}
# 方法2 读取GPL网页的表格文件，按列取子集
##https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL570
if(F){
  #注：表格读取参数、文件列名不统一，活学活用，有的表格里没有symbol列，也有的GPL平台没有提供注释表格，就不能用这种方法。
  b = read.delim("GPL570-55999.txt",
                 check.names = F,
                 comment.char = "#")
   b=data.table::fread("GPL570-55999.txt",data.table=F, skip=17)              
  colnames(b)
  #下一行代码里的列名是从colnames(b)的输出结果里复制的
  ids2 = b[,c("ID","Gene Symbol")]
  colnames(ids2) = c("probe_id","symbol") #和R包里的列名保持一致
  #下面三句代码是将没有对应到基因的探针和对应多个基因的探针去掉
  k1 = ids2$symbol!="";table(k1) 
  k2 = !str_detect(ids2$symbol,"///");table(k2)
  ids2 = ids2[ k1 & k2,]
  # ids = ids2 
  #使用方法二需要将42行F改为T，55行取消注释
}
# 方法3 官网下载注释文件并读取
# 方法4 自主注释，了解一下
#https://mp.weixin.qq.com/s/mrtjpN8yDKUdCSvSUuUwcA

library(GEOquery)
#打破下载时间的限制,改前60秒，改后10w秒
options(timeout = 100000) ##R默认设置，60s下载不完成就会停止
options(scipen = 20)#
a="GSE42872"
eSet = getGEO(a, destdir = '.', getGPL = F)
#网速太慢，下不下来怎么办
#1.从网页上下载/发链接让别人帮忙下，放在工作目录里
#2.试试geoChina,只能下载2019年前的表达芯片数据

class(eSet)
length(eSet)

eSet = eSet[[1]] 
class(eSet)


#(1)提取表达矩阵exp
exp <- exprs(eSet)
dim(exp)
range(exp)#看数据范围决定是否需要log，是否有负值，异常值
#exp = log2(exp+1) #需要log才log
boxplot(exp,las = 2)
#(2)临床信息
pd <- pData(eSet)

#(3)让exp列名与pd的行名顺序完全一致
p = identical(rownames(pd),colnames(exp));p
if(!p) {
  s = intersect(rownames(pd),colnames(exp))
  exp = exp[,s]
  pd = pd[s,]
}

#(4)提取芯片平台编号，后面要根据它来找探针注释
gpl_number <- eSet@annotation;gpl_number
save(pd,exp,gpl_number,file = "step1output.Rdata")

rm(list = ls())  
load(file = "step1output.Rdata")
# 1.Group----
library(stringr)
# 标准流程代码是二分组，多分组数据的分析后面另讲
# 生成Group向量的三种常规方法，三选一，选谁就把第几个逻辑值写成T，另外两个为F。如果三种办法都不适用，可以继续往后写else if
if(F){
  # 第一种方法，有现成的可以用来分组的列
  
}else if(F){
  # 第二种方法，眼睛数，自己生成
  Group = rep(c("Disease","Normal"),each = 10)
}else if(T){
  # 第三种方法，使用字符串处理的函数获取分组
  k = str_detect(pd$title,"Control");table(k)
  Group = ifelse(k,"control","vem")
}
###rep(c("T","N"),time=c(3,5))

# 需要把Group转换成因子，并设置参考水平，指定levels，对照组在前，处理组在后
Group = factor(Group,levels = c("control","vem"))
Group
##检查分组是否正确
data.frame(pd$title,Group)


#2.探针注释的获取-----------------
#捷径
library(tinyarray)
find_anno(gpl_number) #打出找注释的代码
if (!require("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

BiocManager::install("hugene10sttranscriptcluster.db")
library(hugene10sttranscriptcluster.db);ids <- toTable(hugene10sttranscriptclusterSYMBOL)
#四种方法，方法1里找不到就从方法2找，以此类推。
#捷径里面包含了全部的R包、一部分表格、一部分自主注释
#方法1 BioconductorR包(最常用，已全部收入find_anno里面，不用看啦)

save(exp,Group,ids,file = "step2output.Rdata")
###########################################
rm(list = ls()) 
load(file = "step2output.Rdata")
#差异分析
library(limma)
design = model.matrix(~Group)
fit = lmFit(exp,design)
fit = eBayes(fit)
deg = topTable(fit,coef = 2,number = Inf)

#为deg数据框添加几列
#1.加probe_id列，把行名变成一列
library(dplyr)
deg = mutate(deg,probe_id = rownames(deg))
#2.加上探针注释
ids = distinct(ids,symbol,.keep_all = T)
#其他去重方式在zz.去重方式.R
deg = inner_join(deg,ids,by="probe_id")
nrow(deg)

#3.加change列,标记上下调基因
logFC_t = 1
p_t = 0.05
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)
#火山图
library(ggplot2)
ggplot(data = deg, 
       aes(x = logFC, 
           y = -log10(P.Value))) +
  geom_point(alpha=0.4, size=3.5, 
             aes(color=change)) +
  scale_color_manual(values=c("blue", "grey","red"))+
  geom_vline(xintercept=c(-logFC_t,logFC_t),lty=4,col="black",linewidth=0.8) +
  geom_hline(yintercept = -log10(p_t),lty=4,col="black",linewidth=0.8) +
  theme_bw()
##结束