怎么样才能正确的学习生信分析呢？—从学徒做起

rm(list = ls ())
options(stringAsFactors = F)
f='GSE2513_eSet.Rdata'    #先设置一个Rdata赋值给对象f，让f成为一个Rdata对象。
library(GEOquery)
if(!file.exists(f)){
      gset <- getGEO('GSE2513',destdir='.',
                 AnnotGPL = F,    #注释文件
                 getGPL = F,        #平台文件 ）
         save(gset,file=f)   #保存文件到本地          
    }
a=gset[[1]]  ##因为GEO数据集只有一个GPL平台，所以下载到的是一个含有一个元素的list   
dat=exprs(a) ##a现在是一个对象，取a这个对象通过看说明书知道要用exprs这个函数
boxplot(dat,las=2)  ##原始数据做箱式图看数据分布                                         

BiocManager::install('oligo',ask = F,update = F)
library(oligo)
BiocManager::install('pd.hg.u133.plus.2',ask = F,update = F)
library(pd.hg.u133.plus.2)
dir='GSE2513_RAW/'  ##让这个文件赋值给dir这个变量
od=getwd()  ##当前所处路径
setwd(dir)  ##设置为文件路径，为了处理CEL文件
celFiles <- list.celfiles(listGzipped = T) ##为了列出路径下文件名中包含.cel或.CEL的文件，返回9值是一个包含有文件名的list，主要oligo包中的该函数只能列出被压缩的cel文件或是不被压缩的cel文件，不能同时列出，因此需要设定参数listGzipped。
celFiles
affyRaw <- read.celfiles(celFiles) ##读取CEL文件
setwd(od)
eSet <- rma(affyRaw)#rma函数是对表达矩阵进行normalization(mas5也有异曲同工之妙)
eSet#其实严格来说，这个芯片得到表达矩阵后，是需要过滤的，但是我想看看过滤和未过滤的对比，所以这里我就先不过滤了

a=eSet
dat=exprs(a)
dim(dat)
#22283features 12samples
dat[1:4,1:4]#每步都要查看数据结构
par(cex = 0.7)
n.sample=ncol(dat)
if(n.sample>40) par(cex=0.5)
cols <- rainbow(n.sample*1.2)
boxplot(dat,col=cols,main='new expression value',las=2,cex.axis=0.6)
#通过这些boxplot可以         看到各个芯片直接数据是否整齐，看了一下并不是很整齐，所以需要用normalizeBetweenArrays这个函数处理一  下，这步QC也很重要。

dat=normalizeBetweenArrays(dat)#运用quantile的进行组间归一化，如果上一步boxplot样本间不整齐，这   步肯定要执行的。只能在同一个数据集里面使用。
dat[1:4,1:4]
boxplot(dat,col=cols,main='new expression value',las=2,cex.axis=0.6)#可以看到不整齐的那些已经处理好了

dat[1:4,1:4]    #每步都要检查数据
dat=t(dat)    #转置矩阵，这步非常重要，如果之前的QC没有做好，这步容易出现NA值。转置的目的是符合PCA的输入数据格式：数值型数据，每一列是一个指标(例如基因)，每一行是一个观测(例如样本)
dat=as.data.frame(dat)    #转换成矩阵
dim(dat)    #查看一下维度
library(stringr)
data=data.table::fread('gsm_group.txt',data.table=F,header=F)    #多GSM排序不规则，利用stringr进行分组，做PCA要用到分组信息
names(data)=c('GSM','group')    #给数据的两列命变量名
head(data)
group_list=str_match(data$group,'(.*).*')    #正则表达式的哲学人生，请自行看我给出的网站解释
class(group_list)
head(group_list)
group=data.frame(GSM=data$GSM,group=group_list[,2])
head(group)    
#小白的人生就是每处理完一个数据格式，都要查看一下。
b=group[,2]
class(b)
b=as.character(b)    #一定要变因子
class(b)
colnames(b)='group'
table(b)
b$group[b$group=='Pterygium']='pterygium'    #在原始的txt中本来只有两个处理条件的，但是R是区分大小写的，所以如果你有一个样本是Pterygium，它也会认为这个样本和pterygium是不同的处理条件。
dat=cbind(dat,b)    #加上分组信息
library('FactoMineR')
library('factoextra')    #PCA必要用的两个函数
dim(dat)    #我是小白，每步我都要看看维度，才能放心
dat.pca=PCA(dat[,-ncol(dat)],graph=F)
fviz_pca_ind(dat.pca,
            geom.ind='point',
            col.ind=dat$group,
            palette=c('#00AFBB','#E7B800'),
            addEllipses=T,
            legend.title='Groups_PCA')    #画出PCA图

dat[1:4,1:4]    #我是小白，我检查
library(limma)    #利用limma包进行差异分析
design=model.matrix(~0+factor(group_list))    #这是可以通过limma包做到随心所欲的指定任意两组进行比较的代码
class(design)
dim(design)
colnames(design)=levels(factor(group_list))
exprSet=dat
exprSet[1:4,1:4]
rownames(design)=colnames(exprSet)
design    ##design就是我们制作好的分组矩阵，需要根据我们下载的芯片数据的实验设计方案来，此处例子是GSE2513数据集来探究，12个样本分成了两组，当然读者有自己的数据只需要按照同样方法制作即可！
contrast.matrix=makeContrasts('ptergium-conjunctiva',
                              levels=design）    ##contranst.matrix，这是要把pterygium组合conjunctiva组进行差异分析比较。                              
deg=function(exprSet,design,contrast.matrix){
                                  fit=lmFit(exprSet,design)
                                  fit2=contrasts.fit(fit,contrast.matrix)##这一步很重要
                                  fit2=eBayes(fit2)
                                  tempOutput=topTable(fit2,coef=1,n=Inf)
                                  nrDEG=na.omit(tempOutput)
                                  head(nrDEG)
                                  return(nrDEG)
 }
deg=deg(exprSet,design,contrast.matrix)
dim(deg)
head(deg)
save(deg,file = 'deg.Rdata')
load(file = 'deg.Rdata')
head(deg)
bp(dat[rownames(deg)[1],])
#为deg数据框添加1列
#1.加symbol列，把行名变成一列
library(dplyr)
deg <- mutate(deg,symbol=rownames(deg))
head(deg)
nrDEG=deg
head(nrDEG)
attach(nrDEG)#直接使用数据框中的元素，不用在用$
plot(logFC,-log10(P.Value))
library(ggpubr)
df=nrDEG
head(df)
df$v=-log10(P.Value)
dim(df)
ggscatter(df,x='logFC',y='v',size=0.5)
df$change=ifelse(df$P.Value > 0.05,'stable',
                 ifelse( df$logFC > 0.5849625,'up',
                        ifelse( df$logFC < -0.5849625,'down','stable') )
               )     ##logFC的条件照着文章处理
head(df)
table(df$change)
write.csv(df,file='DEGs.csv')

ggscatter(df,x='logFC',y='v',size = 0.5,color = 'change')
ggscatter(df,x='logFC',y='v',color = 'change',size = 0.5,
          label = 'symbol',repel = T,
          label.select = c('NDRG1','TNNI2','JUNB','MYC','CYP26A1','ATF3','UPK1B','RAB38','SERPINB1','ATP1B3','RBP1','TYRP1'),
          palette = c('#00AFBB','#D3D3D3','#FC4E07'))
ggsave('volcano.png')
dev.off()

dat[1:4,1:4]
head(df)
table(df$change)
df1=df[df$change=='up'|df$change=='down',]
dim(df1)
head(df1)
gene_up=df1[df1$change=='up',]
length(rownames(gene_up))
gene_down=df1[df1$change=='down',]
length(rownames(gene_down))
x=gene_up$logFC
head(x)
names(x)=rownames(gene_up)
head(x)
y=gene_down$logFC
head(y)
names(y)=rownames(gene_down)
head(y)
cg=c(names(head(sort(x),50)),names(head(sort(y),50)))#上调基因和下调基因各取50画图。
head(cg)
length(cg)
dat[1:4,1:4]
pheatmap(dat[cg,],show_rownames = F,show_colnames = F)
n=t(scale(t(dat[cg,])))
dim(n)
head(n)
max(n)
min(n)
n[n>2]=2
n[n< -2]=-2
n[1:4,1:4]
pheatmap(n,show_rownames = F,show_colnames = F)
ac=data.frame(g=group_list)
ac
rownames(ac)=colnames(n)
ac
pheatmap(n,
         show_colnames = T,
         show_rownames = T,
         cluster_cols = T,
         #cluster_rows = T,
         border_color = NA,
         annotation_col = ac,
         annotation_legend = T,
         cellwidth = 30,
         cellheight = 8,
         fontsize = 8,
         filename = 'heatmap_up50_down50_DEG.png')