中文字幕理论片,69视频免费在线观看,亚洲成人app,国产1级毛片,刘涛最大尺度戏视频,欧美亚洲美女视频,2021韩国美女仙女屋vip视频

# 首先我下載了一個(gè)錯(cuò)誤的包
BiocManager::install('hgug4110b.db')
library('hgug4110b.db')
?hgug4110b
ls("package:hgug4110b.db")

#下面這段很眼熟，就是我上面在R語言20題里的截圖，我把視頻中的hgu95av2.db包換成我在谷歌搜索到的hgug4110b.db包
library('hgug4110b.db')
ids=toTable(hgug4110bSYMBOL)
length(unique(ids$symbol))
tail(sort(table(ids$symbol)))
table(sort(table(ids$symbol)))
plot(table(sort(table(ids$symbol))))


exprSet<-dat
table(rownames(exprSet) %in% ids$probe_id)
dim(exprSet)
exprSet=exprSet[rownames(exprSet) %in% ids$probe_id,]
dim(exprSet)


ids=ids[match(rownames(exprSet),ids$probe_id),]
head(ids)
exprSet[1:5,1:5]
tmp = by(exprSet,ids$symbol,
         function(x) rownames(x)[which.max(rowMeans(x))] )
probes = as.character(tmp)
exprSet=exprSet[rownames(exprSet) %in% probes ,]
dim(exprSet)
#寫一個(gè)函數(shù)，可以之間過濾掉其他探針并且進(jìn)行ID轉(zhuǎn)換，
#此時(shí)的id轉(zhuǎn)換是用bioconductor包來轉(zhuǎn)換的，如果下載的matrix，則不用轉(zhuǎn)換
jimmy <- function(exprSet,ids){
                             tmp = by(exprSet,
                             ids$symbol,
                              function(x) rownames(x)[which.max(rowMeans(x))] )
                probes = as.character(tmp)
                print(dim(exprSet))
                exprSet=exprSet[rownames(exprSet) %in% probes ,]
                print(dim(exprSet))
                return(exprSet)
}
new.exprSet <- jimmy(exprSet,ids)


#此時(shí)會報(bào)錯(cuò)，因?yàn)楝F(xiàn)在的exprSet和ids的長度不一樣

#代碼如下
gpl <- getGEO('GPL17077', destdir=".")
colnames(Table(gpl)) 
head(Table(gpl)[,c(1,6,7)])  ## you need to check this , which column do you need 
write.csv(Table(gpl)[,c(1,6,7)],"GPL17077.csv")
ids=read.csv('GPL17077.csv')
#這樣我們就同樣獲得了探針對應(yīng)基因的表達(dá)矩陣?yán)?/span>

a=read.table('GSE42872_series_matrix.txt.gz')
> class(a)
[1] "data.frame"
> str(gset)

gset <- getGEO("GSE42589")

library(GEOquery)
gset <- getGEO('GSE42872',destdir=".",AnnotGPL = F,getGPL = F) 

# 就這么來就對了，沒毛病
###step1-downloadR
rm(list = ls())  ## 魔幻操作，一鍵清空~
options(stringsAsFactors = F)
# 注意查看下載文件的大小，檢查數(shù)據(jù) 
f='GSE99507_eSet.Rdata'

library(GEOquery)
# 這個(gè)包需要注意兩個(gè)配置，一般來說自動化的配置是足夠的。
#Setting options('download.file.method.GEOquery'='auto')
#Setting options('GEOquery.inmemory.gpl'=FALSE)
if(!file.exists(f)){
  gset <- getGEO('GSE99507', destdir=".",
                 AnnotGPL = F,     ## 注釋文件
                 getGPL = F)       ## 平臺文件
  save(gset,file=f)   ## 保存到本地
}
load('GSE99507_eSet.Rdata')  ## 載入數(shù)據(jù)
class(gset)
length(gset)
class(gset[[1]])
# 因?yàn)檫@個(gè)GEO數(shù)據(jù)集只有一個(gè)GPL平臺，所以下載到的是一個(gè)含有一個(gè)元素的list
a=gset[[1]]
dat=exprs
dim(dat)
#  GPL3921  [HT_HG-U133A] Affymetrix HT Human Genome U133A Array
# Bioconductor - hgu133a.db
dat[1:4,1:4]


pd=pData(a) #獲得分組信息，就得到了前面一張截圖，即前面懸念的解答處-GSE號對應(yīng)的Samples分組信息，主要是根據(jù)文章里的那個(gè)截圖的分組信息來對字符串進(jìn)行拆分（str_split）和替換（str_replace）
library(stringr)
group_list = trimws(str_split(pd$title,' ',simplify = T)[,5])#拆分
group_list = str_replace(group_list,'LM2','Vector')
table(group_list)
save(dat,group_list,file = 'step1-output.Rdata')
#此時(shí)，要注意保存的這個(gè)文件step1-output.Rdata，我們要知道此時(shí)的dat表達(dá)矩陣是什么樣的

###step2-ids-filter.R

rm(list = ls())  ## 魔幻操作，一鍵清空~
options(stringsAsFactors = F)
load('step1-output.Rdata')
gpl <- getGEO('GPL17077', destdir=".")
colnames(Table(gpl)) ## [1] 41108    17
head(Table(gpl)[,c(1,7)])  ## you need to check this , which column do you need 
write.csv(Table(gpl)[,c(1,7)],"GPL17077.csv")
ids=read.csv('GPL17077.csv')
ids=ids[,-1]
dat[1:4,1:4] 
rownames(dat)

#如前所述，現(xiàn)在沒有這個(gè)平臺對應(yīng)的注視R包，所以下面是根據(jù)下載的dat和ids來進(jìn)行的ID轉(zhuǎn)換步驟
ids[ids$GENE_SYMBOL!='',]
ids=ids[ids$GENE_SYMBOL!='',]
match(rownames(dat),ids$ID)
ids[match(rownames(dat),ids$ID),]
ids=ids[match(rownames(dat),ids$ID),]
ids=ids[complete.cases(ids),]
dat=dat[ids$ID,]

#下面是針對那些有多個(gè)基因去重的步驟，注意理解為什么要取median（中位數(shù)），巧妙
ids$median=apply(dat,1,median)
ids=ids[order(ids$GENE_SYMBOL,ids$median,decreasing = T),]
ids=ids[!duplicated(ids$GENE_SYMBOL),]
dat=dat[ids$ID,]
rownames(dat)=ids$GENE_SYMBOL
dat[1:4,1:4]  
dim(dat)

save(dat,group_list,file = 'step2-output.Rdata')

#此時(shí)此時(shí)，要注意保存的這個(gè)文件step2-output.Rdata，我們要知道此時(shí)的dat表達(dá)矩陣是什么樣的

###step3-pca-heatmap

##PCA主成分分析

(list = ls())  ## 魔幻操作，一鍵清空~
options(stringsAsFactors = F)
load(file = 'step1-output.Rdata')
# 每次都要檢測數(shù)據(jù)
dat[1:4,1:4]
## 下面是畫PCA的必須操作，需要看說明書。
dat=t(dat)
dat=as.data.frame(dat)
dat=cbind(dat,group_list)
library("FactoMineR")
library("factoextra") 
# The variable group_list (index = ) is removed
# before PCA analysis
dat.pca <- PCA(dat[,-ncol(dat)], graph = FALSE)
fviz_pca_ind(dat.pca,repel =T,
             #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"
)



##熱圖：heatmap-sd

rm(list = ls())  ## 魔幻操作，一鍵清空~
options(stringsAsFactors = F)
load(file = 'step1-output.Rdata')
# 每次都要檢測數(shù)據(jù)
dat[1:4,1:4]
cg=names(tail(sort(apply(dat, 1, sd)),1000))
library(pheatmap)
pheatmap(dat[cg,],show_colnames =F,show_rownames = F)
rowMeans(dat[,1:3]) -  rowMeans(dat[,4:6]) 
cg1=names(tail(sort(rowMeans(dat[,1:3]) -  rowMeans(dat[,4:6]) ),500))
cg2=names(head(sort(rowMeans(dat[,1:3]) -  rowMeans(dat[,4:6]) ),500))
cg=c(cg1,cg2)
n=t(scale(t(dat[cg,])))
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)
pheatmap(n,show_colnames =F,
         show_rownames = F,
         cluster_cols = F,
         annotation_col=ac)

##step4-DEG
rm(list = ls())  ## 魔幻操作，一鍵清空~
options(stringsAsFactors = F)
load(file = 'step2-output.Rdata')######要知道此時(shí)是step2中的dat哦
# 每次都要檢測數(shù)據(jù)
dat[1:4,1:4]
table(group_list)

##boxplot圖
boxplot(dat[1,]~group_list)
t.test(dat[1,]~group_list)
bp=function(g){
  library(ggpubr)
  df=data.frame(gene=g,stage=group_list)
  p <- ggboxplot(df, x = "stage", y = "gene",
                 color = "stage", palette = "jco",
                 add = "jitter")
  #  Add p-value
  p + stat_compare_means()
}
bp(dat[1,])
bp(dat[2,])
bp(dat[3,])
bp(dat[4,])

##DEG-limma包

library(limma)
design=model.matrix(~factor( group_list ))
fit=lmFit(dat,design)
fit=eBayes(fit)
options(digits = 4)
#topTable(fit,coef=2,adjust='BH') 
topTable(fit,coef=2,adjust='BH')
bp(dat['KLK10',])
bp(dat['FXYD3',])
deg=topTable(fit,coef=2,adjust='BH',number = Inf)

head(deg) 
save(deg,file = 'deg.Rdata')


## for volcano 
if(T){
  nrDEG=deg
  head(nrDEG)
  attach(nrDEG)
  plot(logFC,-log10(P.Value))
  library(ggpubr)
  df=nrDEG
  df$v= -log10(P.Value)
  ggscatter(df, x = "logFC", y = "v",size=0.5)

  df$g=ifelse(df$P.Value>0.01,'stable',
              ifelse( df$logFC >1.5,'up',
                      ifelse( df$logFC < -1.5,'down','stable') )
  )
  table(df$g)#此時(shí)獲得的上調(diào)基因?yàn)?0，下調(diào)基因?yàn)?45
  df$symbol=rownames(df)
  ggscatter(df, x = "logFC", y = "v",size=0.5,color = 'g')
  ggscatter(df, x = "logFC", y = "v", color = "g",size = 0.5,
            label = "symbol", repel = T,
            #label.select = rownames(df)[df$g != 'stable'] ,
            label.select = rownames(head(head(deg))),
            palette = c("#00AFBB", "#E7B800", "#FC4E07") )

  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") )


}


## 熱圖：for heatmap 
rm(list = ls())  ## 魔幻操作，一鍵清空~
options(stringsAsFactors = F)
if(T){ 
  load(file = 'step2-output.Rdata')
  # 每次都要檢測數(shù)據(jù)
  dat[1:4,1:4]
  table(group_list)
  load(file = 'deg.Rdata')
  x=deg$logFC
  names(x)=rownames(deg)
  cg=c(names(head(sort(x),20)),#老大的代碼中是100、100，由于我獲得上下調(diào)基因少，就改了 
       names(tail(sort(x),145)))
  library(pheatmap)
  pheatmap(dat[cg,],show_colnames =F,show_rownames = F)
  n=t(scale(t(dat[cg,])))
  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)
  pheatmap(n,show_colnames =F,
           show_rownames = F,
           cluster_cols = F,
           annotation_col=ac)


}