一般來說我們獲得的芯片數(shù)據(jù)在進(jìn)行下游分析時(shí)，表達(dá)矩陣一般都應(yīng)該是log2后的結(jié)果，那么如果我手里有一份二代測(cè)序數(shù)據(jù)的代碼，全盤接收，也給log了，而我又沒細(xì)看這個(gè)地方還有一步log的代碼，再加上對(duì)這個(gè)log如果沒有足夠的清醒意識(shí)，不能嚴(yán)格區(qū)分log和不log的應(yīng)用場(chǎng)景，那么會(huì)得到怎樣的結(jié)果呢？數(shù)據(jù)集來自GSE21933，從下載數(shù)據(jù)到后面獲得火山圖的代碼，都在下面啦！

step1-getdata

rm(list = ls())  
options(stringsAsFactors = F)#在調(diào)用as.data.frame的時(shí)，將stringsAsFactors設(shè)置為FALSE可以避免character類型自動(dòng)轉(zhuǎn)化為factor類型
# 注意查看下載文件的大小，檢查數(shù)據(jù) 


if(T){
f='GSE21933_eSet.Rdata'
# https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE38959
library(GEOquery)
# 這個(gè)包需要注意兩個(gè)配置，一般來說自動(dòng)化的配置是足夠的。
#Setting options('download.file.method.GEOquery'='auto')
#Setting options('GEOquery.inmemory.gpl'=FALSE)
if(!file.exists(f)){
  gset <- getGEO('GSE21933', destdir='.',
                 AnnotGPL = F,     ## 注釋文件
                 getGPL = F)       ## 平臺(tái)文件
  save(gset,file=f)   ## 保存到本地
}
load('GSE21933_eSet.Rdata')  ## 載入數(shù)據(jù)
class(gset)  #查看數(shù)據(jù)類型  
length(gset)  #
class(gset[[1]]) 
# 因?yàn)檫@個(gè)GEO數(shù)據(jù)集只有一個(gè)GPL平臺(tái)，所以下載到的是一個(gè)含有一個(gè)元素的list
a=gset[[1]] #
dat=exprs(a) #a現(xiàn)在是一個(gè)對(duì)象，取a這個(gè)對(duì)象通過看說明書知道要用exprs這個(gè)函數(shù)
dim(dat)#看一下dat這個(gè)矩陣的維度，發(fā)現(xiàn) 有 53617 行 

dat[1:4,1:4] #查看dat這個(gè)矩陣的1至4行和1至4列，逗號(hào)前為行，逗號(hào)后為列
pd=pData(a)

tumor=rownames(pd[grepl('T',as.character(pd$title)),])
normal=rownames(pd[grepl('N',as.character(pd$title)),])

dat=dat[,c(tumor,normal)]
dim(dat)
dat[1:4,1:4] 

group_list=c(rep('tumor',length(tumor)),
             rep('normal',length(normal)))
table(group_list)
}

if(T){
  #你要的通過GPL下載
  library(GEOquery)
 # Download GPL file, put it in the current directory, and load it:
    gpl <- getGEO('GPL6254', destdir='.')
  colnames(Table(gpl))  
  head(Table(gpl)[,c(1,6)]) ## you need to check this , which column do you need
  probe2gene=Table(gpl)[,c(1,6)]
  dim(probe2gene)
  #也可以這樣
  # b=read.table('GPL6254.txt',
  #              sep = '\t',quote = '',
  #              fill=T,header = T)
  # colnames(b)
  # probe2gene=b[,c(1,10)]
  probe2gene=probe2gene[probe2gene$GENE_SYMBOL != '',]
  dim(probe2gene)
  head(probe2gene)
  colnames(probe2gene) <- c('probe_id','symbol')
  dim(dat)
  probe2gene <- probe2gene[probe2gene$probe_id %in% rownames(dat),]
  dim(probe2gene)
  dat <- dat[match(probe2gene$probe_id,rownames(dat)),]
  dim(dat)
  dat[1:4,1:4]

}

if(T){

ids <- probe2gene
dat <- as.data.frame(dat)
ids$median=apply(dat,1,median) #ids新建median這一列，列名為median，同時(shí)對(duì)dat這個(gè)矩陣按行操作，取每一行的中位數(shù)，將結(jié)果給到median這一列的每一行
ids=ids[order(ids$symbol,ids$median,decreasing = T),]#對(duì)ids$symbol按照ids$median中位數(shù)從大到小排列的順序排序，將對(duì)應(yīng)的行賦值為一個(gè)新的ids
ids=ids[!duplicated(ids$symbol),]#將symbol這一列取取出重復(fù)項(xiàng)，'!'為否，即取出不重復(fù)的項(xiàng)，去除重復(fù)的gene ，保留每個(gè)基因最大表達(dá)量結(jié)果s
dat=dat[ids$probe_id,] #新的ids取出probe_id這一列，將dat按照取出的這一列中的每一行組成一個(gè)新的dat
rownames(dat)=ids$symbol#把ids的symbol這一列中的每一行給dat作為dat的行名
dat[1:4,1:4]  #保留每個(gè)基因ID第一次出現(xiàn)的信息
dim(dat)
boxplot(dat)
dat=log(dat 1)
dat[1:4,1:4] 

save(dat,ids,group_list,file = 'step1-getdata.Rdata')
load(file = 'step1-getdata.Rdata')

}

上面得到的表達(dá)矩陣dat和ids如下

step2-check

rm(list = ls())  
options(stringsAsFactors = F)
load(file = 'step1-getdata.Rdata')
# 每次都要檢測(cè)數(shù)據(jù)
dat[1:4,1:4]
## 下面是畫PCA的必須操作，需要看說明書。
dat=t(dat)#畫PCA圖時(shí)要求是行名時(shí)樣本名，列名時(shí)探針名，因此此時(shí)需要轉(zhuǎn)換
dat=as.data.frame(dat)#將matrix轉(zhuǎn)換為data.frame
dat=cbind(dat,group_list) #cbind橫向追加，即將分組信息追加到最后一列
library('FactoMineR')#畫主成分分析圖需要加載這兩個(gè)包
library('factoextra') 
# The variable group_list (index = 54676) is removed
# before PCA analysis
dat.pca <- PCA(dat[,-ncol(dat)], graph = FALSE)#現(xiàn)在dat最后一列是group_list，需要重新賦值給一個(gè)dat.pca,這個(gè)矩陣是不含有分組信息的
fviz_pca_ind(dat.pca,
             geom.ind = 'point', # show points only (nbut not 'text')
             col.ind = dat$group_list, # color by groups
             palette = c('#00AFBB', '#E7B800'),
             addEllipses = TRUE, # Concentration ellipses
             legend.title = 'Groups'
)
ggsave('all_samples_PCA_by_pCR.png')




rm(list = ls())  ## 魔幻操作，一鍵清空~
load(file = 'step1-output.Rdata') #此步為一個(gè)小插曲，即計(jì)算一下從第一行開是計(jì)算每一行的sd值，知道最后一行所需要的時(shí)間
dat[1:4,1:4] 

cg=names(tail(sort(apply(dat,1,sd)),1000))#apply按行（'1'是按行取，'2'是按列?。┤∶恳恍械姆讲睿瑥男〉酱笈判?，取最大的1000個(gè)
library(pheatmap)
pheatmap(dat[cg,],show_colnames =F,show_rownames = F) #對(duì)那些提取出來的1000個(gè)基因所在的每一行取出，組合起來為一個(gè)新的表達(dá)矩陣
n=t(scale(t(dat[cg,]))) # 'scale'可以對(duì)log-ratio數(shù)值進(jìn)行歸一化
n[n>2]=2 
n[n< -2]= -2
n[1:4,1:4]
pheatmap(n,show_colnames =F,show_rownames = F)
ac=data.frame(g=group_list)
rownames(ac)=colnames(n) #把a(bǔ)c的行名給到n的列名，即對(duì)每一個(gè)探針標(biāo)記上分組信息（是'noTNBC'還是'TNBC'）
pheatmap(n,show_colnames =F,show_rownames = F,
         annotation_col=ac,filename = 'heatmap_top1000_sd.png')

得到的PCA圖如下，PCA圖顯示兩組分得很開，那么繼續(xù)往下做差異分析

step3-DEG

#差異分析
if(T){
design <- model.matrix(~0 factor(group_list))
colnames(design)=levels(factor(group_list))
head(design)
exprSet=dat
rownames(design)=colnames(exprSet)
design
contrast.matrix<-makeContrasts('tumor-normal',
                               levels = design)
contrast.matrix ##這個(gè)矩陣聲明，我們要把 Tumor 組跟 Normal 進(jìn)行差異分析比較


colnames(design)
deg = function(exprSet,design,contrast.matrix){
  ##step1
  fit <- lmFit(exprSet,design)
  ##step2
  fit2 <- contrasts.fit(fit, contrast.matrix) 
  ##這一步很重要，大家可以自行看看效果

  fit2 <- eBayes(fit2)  ## default no trend !!!
  ##eBayes() with trend=TRUE
  ##step3
  tempOutput = topTable(fit2, coef=1, n=Inf)
  nrDEG = na.omit(tempOutput) 
  #write.csv(nrDEG2,'limma_notrend.results.csv',quote = F)
  head(nrDEG)
  return(nrDEG)
}
deg = deg(exprSet,design,contrast.matrix)
head(deg)
save(deg,file = 'deg.Rdata')
}

看了一下，deg的結(jié)果，如下兩張圖

感覺不太對(duì)勁，順勢(shì)做下火山圖，

## for volcano 
if(T){
  nrDEG=deg
  head(nrDEG)
  attach(nrDEG)
  plot(logFC,-log10(P.Value))
  library(ggpubr)
  df=nrDEG
  df$v= -log10(P.Value) #df新增加一列'v',值為-log10(P.Value)
  ggscatter(df, x = 'logFC', y = 'v',size=0.5)

  df$g=ifelse(df$P.Value>0.05,'stable', #if 判斷：如果這一基因的P.Value>0.01，則為stable基因
              ifelse( df$logFC >1,'up', #接上句else 否則：接下來開始判斷那些P.Value<0.01的基因，再if 判斷：如果logFC >1.5,則為up（上調(diào)）基因
                      ifelse( df$logFC < -1,'down','stable') )#接上句else 否則：接下來開始判斷那些logFC <1.5 的基因，再if 判斷：如果logFC <1.5，則為down（下調(diào)）基因，否則為stable基因
  )
  table(df$g)
  df$name=rownames(df)
  head(df)
  ggscatter(df, x = 'logFC', y = 'v',size=0.5)
    ggscatter(df, x = 'logFC', y = 'v', color = 'g',size = 0.5,
            label = 'name', repel = T,
            #label.select = rownames(df)[df$g != 'stable'] ,
            label.select =head(rownames(deg)), #挑選一些基因在圖中顯示出來
            palette = c('#00AFBB', '#E7B800', '#FC4E07') )
  ggsave('volcano.png')

  ggscatter(df, x = 'AveExpr', y = 'logFC',size = 0.2)
  df$p_c = ifelse(df$P.Value<0.001,'p<0.001',
                  ifelse(df$P.Value<0.01,'0.001<p<0.01','p>0.01'))
  table(df$p_c )
  ggscatter(df,x = 'AveExpr', y = 'logFC', color = 'p_c',size=0.2, 
            palette = c('green', 'red', 'black') )
  ggsave('MA.png')


}

得到了火山圖，媽耶！

這也太好看了，如果上天告訴我實(shí)際上就只有這兩個(gè)差異基因，那我一定倍加珍惜它們！

可是，實(shí)際上并不是如此！我們?cè)倩剡^頭去看第一張dat的矩陣，范圍就是10左右，并沒有成敗上千和0的，而成百上千的數(shù)據(jù)，但未經(jīng)歸一化才出現(xiàn)既有成敗上千的數(shù)字又有表達(dá)量為0的。為了確定一下，我們?cè)偃PL平臺(tái)看看，到底是芯片還是二代測(cè)序數(shù)據(jù)，如下圖，所以我們不需要log了！
dat的矩陣在log如果boxplot以后看到，就這么規(guī)整，那么我們就不要再log一次了！所以如果是這樣的的矩陣，又不小心給log了，后果就是差異分析的結(jié)果變成logFC>1的有一個(gè)，logFC<1的也只有1個(gè)，就是上面的火山圖的樣子了。
最后，第一步step1-getdata中找到倒數(shù)第四行字dat=log(dat 1)去掉，再重新運(yùn)行一遍，就可以了！

下面是火山圖的分解代碼

plot(logFC,-log10(P.Value))

ggscatter(df, x = 'logFC', y = 'v',size=0.5)

ggscatter(df, x = 'logFC', y = 'v',size=0.5,color = 'g')

df$g=ifelse(df$P.Value>0.05,'stable', #if 判斷：如果這一基因的P.Value>0.01，則為stable基因
              ifelse( df$logFC >1,'up', #接上句else 否則：接下來開始判斷那些P.Value<0.01的基因，再if 判斷：如果logFC >1.5,則為up（上調(diào)）基因
                      ifelse( df$logFC < -1,'down','stable') )#接上句else 否則：接下來開始判斷那些logFC <1.5 的基因，再if 判斷：如果logFC <1.5，則為down（下調(diào)）基因，否則為stable基因
  )
ggscatter(df, x = 'logFC', y = 'v', color = 'g',size = 0.5,
            label = 'name', repel = T,
            #label.select = rownames(df)[df$g != 'stable'] ,
            label.select =head(rownames(deg)), #挑選一些基因在圖中顯示出來
            palette = c('#00AFBB', '#E7B800', '#FC4E07') )
  ggsave('volcano.png')

ggscatter(df, x = 'AveExpr', y = 'logFC',size = 0.2)

  df$p_c = ifelse(df$P.Value<0.001,'p<0.001',
                  ifelse(df$P.Value<0.01,'0.001<p<0.01','p>0.01'))
  table(df$p_c )
  ggscatter(df,x = 'AveExpr', y = 'logFC', color = 'p_c',size=0.2, 
            palette = c('green', 'red', 'black') )

很好的例子，可以理解了，什么時(shí)候該log，什么時(shí)候不log。

本站僅提供存儲(chǔ)服務(wù)，所有內(nèi)容均由用戶發(fā)布，如發(fā)現(xiàn)有害或侵權(quán)內(nèi)容，請(qǐng)點(diǎn)擊舉報(bào)。

免费视频淫片aa毛片_日韩高清在线亚洲专区vr_日韩大片免费观看视频播放_亚洲欧美国产精品完整版

step1-getdata

step2-check

step3-DEG