R_test

題目：
初級：http://www.bio-info-trainee.com/3793.html
中級：http://www.bio-info-trainee.com/3750.html
高級： http://www.bio-info-trainee.com/3415.html

basic

#1打開 `Rstudio` 告訴我它的工作目錄野哭。
getwd()
#2新建6個(gè)向量姻乓，基于不同的`原子類型`胰默。（重點(diǎn)是字符串，數(shù)值话浇，邏輯值）
R中的原子型向量可以存儲以下6種數(shù)據(jù)類型：雙整型(double)，整型(integer)闹究，字符型(character)幔崖，邏輯型(logical)，復(fù)數(shù)類型(complex)渣淤，原始類型(raw)
dou <- c(1,2,3,4,5,6)
int <- c(1:5)
text <- c("Hello", "World")
logic <- c(TRUE, FALSE, TRUE)
comp <- c(1+1i, 2+2i, 1+3i)
raw(3)

#3告訴我在你打開的rstudio里面 `getwd()` 代碼運(yùn)行后返回的是什么赏寇？

#4新建一些數(shù)據(jù)結(jié)構(gòu)，比如矩陣砂代，數(shù)組蹋订，數(shù)據(jù)框，列表等重點(diǎn)是數(shù)據(jù)框刻伊，矩陣）
矩陣：test <- matrix(c(1:9),nrow = 3)
數(shù)組：test <- array(1:27,c(3,3,3))
數(shù)據(jù)框：test <- data.frame(c(1,2,3),c(4,5,6),c(7,8,9))
列表：list <- list(one=c(1,2,3),two=c(4,5,6),three=c(7,8,9))

#5在你新建的數(shù)據(jù)框進(jìn)行切片操作露戒，比如首先取第1，3行捶箱， 然后取第4智什，6列
取第一行：test[c(1),]
取第二列：test[,c(2)]

#6使用data函數(shù)來加載R內(nèi)置數(shù)據(jù)集 `rivers` 描述它。并且可以查看更多的R語言內(nèi)置的數(shù)據(jù)集：[https://mp.weixin.qq.com/s/dZPbCXccTzuj0KkOL7R31g]
data('rivers')

#7下載 [https://www.ncbi.nlm.nih.gov/sra?term=SRP133642](https://www.ncbi.nlm.nih.gov/sra?term=SRP133642) 里面的 `RunInfo Table` 文件讀入到R里面丁屎，了解這個(gè)數(shù)據(jù)框荠锭，多少列，每一列都是什么屬性的元素晨川。
SRT <- read.csv("D:/RStudio/working_path/SraRunTable.txt")
dim(SRT)

#8下載 [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE111229](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE111229) 里面的`樣本信息sample.csv`讀入到R里面证九，了解這個(gè)數(shù)據(jù)框删豺，多少列，每一列都是什么屬性的元素愧怜。（參考 [https://mp.weixin.qq.com/s/fbHMNXOdwiQX5BAlci8brA]) 獲取樣本信息sample.csv）)
library(GEOquery)
GSE_name = 'GSE111229'
options( 'download.file.method.GEOquery' = 'libcurl' ) 
gset <- getGEO( GSE_name, getGPL = F )
save( gset, file = 'gset.Rdata' )
load("gset.Rdata")
Gset <- gset[[1]]
GEO<-pData(Gset)

#9把前面兩個(gè)步驟的兩個(gè)表（`RunInfo Table` 文件呀页，樣本信息sample.csv）關(guān)聯(lián)起來，使用merge函數(shù)拥坛。
identical(SRT$Sample.Name , GEO$geo_accession)  
m=merge(SRT,GEO,by.x = 'Sample.Name',by.y = 'geo_accession')

#10對前面讀取的 RunInfo Table 文件在R里面探索其MBases列蓬蝶，包括 箱線圖(boxplot)和五分位數(shù)(fivenum)，還有頻數(shù)圖(hist)猜惋，以及密度圖(density) 丸氛。
e=m[,c("Bases","title")]
boxplot(e[[1]])
fivenum(e[[1]])
hist(e[[1]])
plot(density(e[[1]]))

#11把前面讀取的樣本信息表格的樣本名字根據(jù)下劃線分割看第3列元素的統(tǒng)計(jì)情況(第三列代表該樣本所在的plate)
plate=unlist(lapply(e[,2],function(x){
 x
 strsplit(x,'_')[[1]][3]
}))
table(plate)

#12根據(jù)plate把關(guān)聯(lián)到的 RunInfo Table 信息的MBases列分組檢驗(yàn)是否有統(tǒng)計(jì)學(xué)顯著的差異。
t.test(e[,1]~plate)
e$plate=plate
#分組繪制箱線圖(boxplot)著摔，頻數(shù)圖(hist)缓窜，以及密度圖(density) 
boxplot(e[,1]~plate)

#13使用ggplot2把上面的圖進(jìn)行重新繪制。
library(ggplot2)
colnames(e)
ggplot(e,aes(x=plate,y=Bases))+geom_boxplot()

#14使用ggpubr把上面的圖進(jìn)行重新繪制谍咆。
library(ggpubr)
p <- ggboxplot(e, x = "plate", y = "MBases",
 color = "plate", palette = "jco",
 add = "jitter")
# Add p-value
p + stat_compare_means(method = 't.test')

medium

1. 根據(jù)R包org.Hs.eg.db找到下面ensembl 基因ID 對應(yīng)的基因名(symbol)

library(org.Hs.eg.db)
g2s=toTable(org.Hs.egSYMBOL)
g2e=toTable(org.Hs.egENSEMBL)
#taget.txt
ensembl_id
ENSG00000000003.13
ENSG00000000005.5
ENSG00000000419.11
ENSG00000000457.12
ENSG00000000460.15
ENSG00000000938.11
t <- read.csv("D:/RStudio/working_path/target.txt", sep="")
n=unlist(lapply(t[,1],function(x){
 x
 strsplit(x,'[.]')[[1]][1]}))
m=data.frame(ensembl_id=n)
a=merge(g2e,m,by="ensembl_id")
b=merge(a,g2s,by="gene_id")
圖片.png

2. 根據(jù)R包hgu133a.db找到下面探針對應(yīng)的基因名(symbol)

3. 找到R包CLL內(nèi)置的數(shù)據(jù)集的表達(dá)矩陣?yán)锩娴腡P53基因的表達(dá)量雹洗，并且繪制在 progres.-stable分組的boxplot圖

library(CLL)
data(sCLLex)
exprSet=exprs(sCLLex) 
CLL=pData(sCLLex)
library(hgu95av2.db)
ids=toTable(hgu95av2SYMBOL)
probe_tp53 <- c("1939_at","1974_s_at","31618_at")
boxplot(exprSet[probe_tp53[i],] ~ CLL$Disease)
圖片.png

4. 找到BRCA1基因在TCGA數(shù)據(jù)庫的乳腺癌數(shù)據(jù)集([Breast Invasive Carcinoma (TCGA, PanCancer Atlas)]的表達(dá)情況

(http://www.cbioportal.org/datasets)

rm(list=ls())
options(stringsAsFactors = F)
f <- read.csv("plot.txt", sep = "\t")
colnames(f) <- c("id", "subtype", "expression", "mut")
library(ggstatsplot)
ggbetweenstats(data = f, 
               x = subtype,
               y = expression)
library(ggplot2)
ggsave("e4-BRCA1-TCGA.png")
圖片.png

圖片.png

5找到TP53基因在TCGA數(shù)據(jù)庫的乳腺癌數(shù)據(jù)集的表達(dá)量分組看其是否影響生存

(http://www.oncolnc.org/)

rm(list=ls())
options(stringsAsFactors = F)
f <- read.csv('BRCA_7157_50_50.csv')
cp <- f
library(ggstatsplot)
ggbetweenstats(data = cp,
               x = Group,
               y = Expression)
library(ggplot2)
library(survival)
library(survminer)
cp$Status <- ifelse(cp$Status == "Dead", 1, 0)
survf <- survfit(Surv(Days,Status)~Group, data=cp)
ggsurvplot(survf, conf.int = F, pval = T)
# complex survplot
ggsurvplot(survf,palette = c("orange", "navyblue"),
           risk.table = T, pval = T,
           conf.int = T, xlab = "Time of days",
           ggtheme = theme_light(),
           ncensor.plot = T)
ggsave("survival_TP53_in_BRCA_taga.png")
圖片.png

6. 下載數(shù)據(jù)集GSE17215的表達(dá)矩陣并且提取下面的基因畫熱圖

ACTR3B ANLN BAG1 BCL2 BIRC5 BLVRA CCNB1 CCNE1 CDC20 CDC6 CDCA1 CDH3 CENPF CEP55 CXXC5 EGFR ERBB2 ESR1 EXO1 FGFR4 FOXA1 FOXC1 GPR160 GRB7 KIF2C KNTC2 KRT14 KRT17 KRT5 MAPT MDM2 MELK MIA MKI67 MLPH MMP11 MYBL2 MYC NAT1 ORC6L PGR PHGDH PTTG1 RRM2 SFRP1 SLC39A6 TMEM45B TYMS UBE2C UBE2T

rm(list=ls())
options(stringsAsFactors = F)
#下載樣本
library(GEOquery)
GSE_name = 'GSE17215'
options( 'download.file.method.GEOquery' = 'libcurl' ) 
geo <- getGEO( GSE_name, getGPL = F )
save( geo, file = paste0(GSE_name,'.set.Rdata'))
load(paste0(GSE_name,'.set.Rdata'))
expr <- exprs(geo[[1]])
#提取目標(biāo)基因
library(hgu133a.db)
p2s=toTable(hgu133aSYMBOL)
#剔除沒有注釋的id
#expr <- expr[p2s$probe_id,] 
target <- "ACTR3B ANLN BAG1 BCL2 BIRC5 BLVRA CCNB1 CCNE1 CDC20 CDC6 CDCA1 CDH3 CENPF CEP55 CXXC5 EGFR ERBB2 ESR1 EXO1 FGFR4 FOXA1 FOXC1 GPR160 GRB7 KIF2C KNTC2 KRT14 KRT17 KRT5 MAPT MDM2 MELK MIA MKI67 MLPH MMP11 MYBL2 MYC NAT1 ORC6L PGR PHGDH PTTG1 RRM2 SFRP1 SLC39A6 TMEM45B TYMS UBE2C UBE2T"
target <- strsplit(target, ' ')[[1]]
tmp <- dplyr::filter(p2s, p2s$symbol %in% target)
tmp2 <- tibble::rownames_to_column(data.frame(expr),"probe_id")
tmp3 <- merge(tmp,tmp2,by="probe_id")
#計(jì)算每個(gè)probe的表達(dá)均值并畫熱圖
##apply函數(shù)經(jīng)常用來計(jì)算矩陣中行或列的均值、和值的函數(shù)
##第一個(gè)參數(shù)是指要參與計(jì)算的矩陣卧波；
##第二個(gè)參數(shù)是指按行計(jì)算還是按列計(jì)算时肿，1——表示按行計(jì)算，2——按列計(jì)算港粱；
##第三個(gè)參數(shù)是指具體的運(yùn)算參數(shù)
tmp3$median <- apply(tmp3[,3:ncol(tmp3)], 1, median)
#只保留平均表達(dá)量最大的探針
tmp4 <- tmp3[order(tmp3$symbol, tmp3$median, decreasing = F),]
#去重復(fù)螃成，只保留表達(dá)量最大的一列
tmp5 <- tmp4[!duplicated(tmp4$symbol),]
#將每一行以基因名命名
rownames(tmp5) <- tmp5$symbol
#除去“probe_id" "symbol" "median" 三列
tmp6 <- tmp5[,-c(1,2,ncol(tmp5))]
#tmp6 <- tmp5[,c(3:8)]
log2 <- log2(tmp6)
pheatmap::pheatmap(log2, scale = "row")
圖片.png

7. 下載數(shù)據(jù)集GSE24673的表達(dá)矩陣計(jì)算樣本的相關(guān)性并且繪制熱圖，需要標(biāo)記上樣本分組信息

rm(list=ls())
options(stringsAsFactors = F)
#下載樣本
library(GEOquery)
GSE_name = 'GSE24673'
options( 'download.file.method.GEOquery' = 'libcurl' ) 
geo <- getGEO( GSE_name, getGPL = F )
save( geo, file = paste0(GSE_name,'.set.Rdata'))
load(paste0(GSE_name,'.set.Rdata'))
expr <- exprs(geo[[1]])
pdata <- pData(geo[[1]]) 
# 根據(jù)pdata第八列做一個(gè)分組信息矩陣
grp <- pdata[,8]
grp_df <- data.frame(group=grp)
rownames(grp_df) <- colnames(expr)
cor <- cor(expr)
pheatmap::pheatmap(cor, annotation_col = grp_df)
圖片.png

8. 找到 GPL6244 platform of Affymetrix Human Gene 1.0 ST Array 對應(yīng)的R的bioconductor注釋包查坪，并且安裝它寸宏！

通過GSE序列號下載數(shù)據(jù)→從得到的list中獲取GPL序列號→找到對應(yīng)的bioconda注釋包
https://mp.weixin.qq.com/s/sVSsV2fWZOQmNd72Vb3SmQ

9. 下載數(shù)據(jù)集GSE42872的表達(dá)矩陣，并且分別挑選出 所有樣本的(平均表達(dá)量/sd/mad/)最大的探針偿曙，并且找到它們對應(yīng)的基因氮凝。

SD:Standard Deviation
MAD：Median absolute deviation，就是先求出給定數(shù)據(jù)的中位數(shù)(注意并非均值)然后原數(shù)列的每個(gè)值與這個(gè)中位數(shù)求出絕對差,然后新數(shù)列的中位值就是MAD望忆。平均絕對誤差可以避免誤差相互抵消的問題罩阵，因而可以準(zhǔn)確反映實(shí)際預(yù)測誤差的大小。

rm(list=ls())
options(stringsAsFactors = F)
#下載樣本
library(GEOquery)
GSE_name = 'GSE42872'
options( 'download.file.method.GEOquery' = 'libcurl' ) 
geo <- getGEO( GSE_name, getGPL = F )
save( geo, file = paste0(GSE_name,'.set.Rdata'))
load(paste0(GSE_name,'.set.Rdata'))
expr <- exprs(geo[[1]])
library(hugene10sttranscriptcluster.db)
p2s=toTable(hugene10sttranscriptclusterSYMBOL)
#剔除沒有注釋的id
expr <- expr[p2s$probe_id,] 
mean <- sort(apply(expr,1,mean),decreasing = T)[1]
p2s[which(p2s == names(mean)),]
sd <- sort(apply(expr,1,sd),decreasing = T)[1]
p2s[which(p2s == names(sd)),]
mad <- sort(apply(expr,1,mad),decreasing = T)[1]
p2s[which(p2s == names(mad)),]
圖片.png

10. 下載數(shù)據(jù)集GSE42872的表達(dá)矩陣启摄，并且根據(jù)分組使用limma做差異分析稿壁，得到差異結(jié)果矩陣

rm(list=ls())
options(stringsAsFactors = F)
#下載樣本
library(GEOquery)
GSE_name = 'GSE42872'
options( 'download.file.method.GEOquery' = 'libcurl' ) 
geo <- getGEO( GSE_name, getGPL = F )
save( geo, file = paste0(GSE_name,'.set.Rdata'))
load(paste0(GSE_name,'.set.Rdata'))
expr <- exprs(geo[[1]])
pdata <- pData(geo[[1]]) 

grp <- unlist(lapply(pdata$title, function(x){
  strsplit(x, ' ')[[1]][4]
}))

library(limma)
#limma needs：表達(dá)矩陣（expr：需要用log后的矩陣）、分組矩陣（design）歉备、比較矩陣（contrast）
#表達(dá)矩陣:expr
logexpr <- log2(expr)
#分組矩陣:design
design <- model.matrix(~0+factor(grp))
colnames(design) <- levels(factor(grp))
rownames(design) <- colnames(expr)
#比較矩陣:contrast
contrast<-makeContrasts(paste0(unique(grp),collapse = "-"),levels = design)
contrast
#limma差異分析流程
fit <- lmFit(logexpr,design)
fit2 <- contrasts.fit(fit, contrast) 
fit3 <- eBayes(fit2)  
output = topTable(fit3, coef=1, n=Inf)
DEG = na.omit(output) 
# 火山圖
if(T){
  logFC_cutoff <- with(DEG,mean(abs(logFC)) + 2*sd(abs(logFC)) )
  
  DEG$change = as.factor(ifelse(DEG$P.Value < 0.05 & abs(DEG$logFC) > logFC_cutoff,
                                ifelse(DEG$logFC > logFC_cutoff ,'UP','DOWN'),'NOT'))
  title <- paste0('log2FoldChange cutoff: ',round(logFC_cutoff,3),
                  '\nUp-regulated genes: ',nrow(DEG[DEG$change =='UP',]) ,
                  '\nDown-regulated genes: ',nrow(DEG[DEG$change =='DOWN',]))
}
library(ggplot2)
plot = ggplot(data=DEG, aes(x=logFC, y=-log10(P.Value), color=change)) +
  geom_point(alpha=0.4, size=1.75) +
  theme_set(theme_set(theme_bw(base_size=20)))+
  xlab("log2FoldChange") + ylab("-log10 p-value") +
  ggtitle( title ) + theme(plot.title = element_text(size=10,hjust = 0.8))+
  scale_colour_manual(values = c('blue','black','red')) # according to the levels(DEG$change)
print(plot)
圖片.png

advanced

library(CLL)
data(sCLLex)
exprSet=exprs(sCLLex)   
samples=sampleNames(sCLLex)
pdata=pData(sCLLex)
group_list=as.character(pdata[,2])
library(hgu95av2.db)
ids=toTable(hgu95av2SYMBOL)
#6.共有多少個(gè)gene
length(unique(ids$symbol))
#一個(gè)基因?qū)?yīng)多少個(gè)probe
table(sort(table(ids$symbol)))
plot(table(sort(table(ids$symbol))))
#7. expr中有多少個(gè)probe沒在注釋包中  
table(rownames(exprSet) %in% ids$probe_id)
#8. 剔除expr中注釋包沒有的探針
exprSet=exprSet[rownames(exprSet) %in% ids$probe_id,]
#沒懂,再一次對expr和ids匹配傅是？
ids=ids[match(rownames(exprSet),ids$probe_id),]
#9. 只保留expr中所有樣品平均表達(dá)量最大的探針
##第一種方法
###by(data, INDICES, FUN)函數(shù)的典型用法： 是將data數(shù)據(jù)框或矩陣按照INDICES因子水平進(jìn)行分組，然后對每組應(yīng)用FUN函數(shù)
tmp = by(exprSet,ids$symbol,
function(x) rownames(x)[which.max(rowMeans(x))] )
probes <-  as.character(tmp)
exprSet_filtered=exprSet[rownames(exprSet) %in% probes ,]
###match返回值為匹配上的probe在ids中的行數(shù)
rownames(exprSet_filtered)=ids[match(rownames(exprSet_filtered),ids$probe_id),2]
exprSet <- exprSet_filtered
##第二種方法
###ids新建median這一列，列名為median喧笔，同時(shí)對dat這個(gè)矩陣按行操作帽驯，取每一行的中位數(shù)，將結(jié)果給到median這一列的每一行
ids$median=apply(exprSet,1,median) 
###按symbol分組后
ids=ids[order(ids$symbol,ids$median,decreasing = T),]
###將symbol這一列取取出重復(fù)項(xiàng)书闸，'!'為否界拦，即取出不重復(fù)的項(xiàng)，去除重復(fù)的gene 梗劫，保留每個(gè)基因最大表達(dá)量結(jié)果
ids=ids[!duplicated(ids$symbol),]
###按照新probe_id，取出exprSet中的各樣品表達(dá)量
exprSet=exprSet[ids$probe_id,]
rownames(exprSet)=ids$symbol
#對表達(dá)量作圖
boxplot(exprSet[,1])
hist(exprSet[,1])
density(exprSet[,1])
boxplot(exprSet['GAPDH',])
boxplot

hist
library(reshape2)
exprSet_L=melt(exprSet)
colnames(exprSet_L)=c('probe','sample','value')
exprSet_L$group=rep(group_list,each=nrow(exprSet))
library(ggplot2) 
#rep(x, time = , length = , each = ,)
#each：代表的是對向量中的每個(gè)元素進(jìn)行復(fù)制的次數(shù)
p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_boxplot()
print(p)
p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_violin()
print(p)
p=ggplot(exprSet_L,aes(value,fill=group))+geom_histogram(bins = 200)+facet_wrap(~sample, nrow = 4)
print(p)
p=ggplot(exprSet_L,aes(value,col=group))+geom_density()+facet_wrap(~sample, nrow = 4)
print(p)
p=ggplot(exprSet_L,aes(value,col=group))+geom_density() 
print(p)
p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_boxplot()
p=p+stat_summary(fun.y="mean",geom="point",shape=23,size=3,fill="red")
p=p+theme_set(theme_set(theme_bw(base_size=20)))
p=p+theme(text=element_text(face='bold'),axis.text.x=element_text(angle=30,hjust=1),axis.title=element_blank())
print(p)
boxplot

violin

histogram

density

density

圖片.png
##12 mean,median,max,min,sd,var,mad
g_mean <- tail(sort(apply(exprSet,1,mean)),50)
g_median <- tail(sort(apply(exprSet,1,median)),50)
g_max <- tail(sort(apply(exprSet,1,max)),50)
g_min <- tail(sort(apply(exprSet,1,min)),50)
g_sd <- tail(sort(apply(exprSet,1,sd)),50)
g_var <- tail(sort(apply(exprSet,1,var)),50)
g_mad <- tail(sort(apply(exprSet,1,mad)),50)
g_mad
names(g_mad)
##13 heatmap 
library(pheatmap)
choose_gene=names(g_mad)
choose_matrix=exprSet[choose_gene,]
#scale()函數(shù)按列計(jì)算截碴，所以先用t()進(jìn)行矩陣倒置
choose_matrix=t(scale(t(choose_matrix)))
pheatmap(choose_matrix)
圖片.png
##14 使用UpSetR包來看他們之間的overlap情況梳侨。
library(UpSetR)
g_all <- unique(c(names(g_mean),names(g_median),names(g_max),names(g_min),
                  names(g_sd),names(g_var),names(g_mad) ))
dat=data.frame(g_all=g_all,
               g_mean=ifelse(g_all %in%  names(g_mean) ,1,0),
               g_median=ifelse(g_all %in%  names(g_median) ,1,0),
               g_max=ifelse(g_all %in%  names(g_max) ,1,0),
               g_min=ifelse(g_all %in%  names(g_min) ,1,0),
               g_sd=ifelse(g_all %in%  names(g_sd) ,1,0),
               g_var=ifelse(g_all %in%  names(g_var) ,1,0),
               g_mad=ifelse(g_all %in%  names(g_mad) ,1,0)
)
upset(dat,nsets = 7)
overlap
#16 對樣品聚類并且繪圖，然后添加樣本的臨床表型數(shù)據(jù)信息（更改樣本名）
colnames(exprSet)=paste(group_list,1:22,sep='')
# Define nodePar
nodePar <- list(lab.cex = 0.6, pch = c(NA, 19), 
                cex = 0.7, col = "blue")
hc=hclust(dist(t(exprSet)))
#par()修改plot邊距
par(mar=c(5,5,5,10)) 
plot(as.dendrogram(hc), nodePar = nodePar,  horiz = TRUE)
cluster
#17 PCA 
library(ggfortify)
df=as.data.frame(t(exprSet))
df$group=group_list 
autoplot(prcomp( df[,1:(ncol(df)-1)] ), data=df,colour = 'group')

library("FactoMineR")#畫主成分分析圖需要加載這兩個(gè)包
library("factoextra") 
#現(xiàn)在dat最后一列是group_list日丹，需要重新賦值給一個(gè)dat.pca,這個(gè)矩陣是不含有分組信息的
df=as.data.frame(t(exprSet))
dat.pca <- PCA(df, graph = FALSE)
fviz_pca_ind(dat.pca,
             geom.ind = "point", # show points only (nbut not "text")
             col.ind = group_list, # color by groups
             # palette = c("#00AFBB", "#E7B800"),
             addEllipses = TRUE, # Concentration ellipses
             legend.title = "Groups"
)
PCA1

PCA2
#18 t.test  
dat = exprSet
group_list=as.factor(group_list)
group1 = which(group_list == levels(group_list)[1])
group2 = which(group_list == levels(group_list)[2])
dat1 = dat[, group1]
dat2 = dat[, group2]
dat3 = cbind(dat1, dat2)
pvals = apply(exprSet, 1, function(x){
  t.test(as.numeric(x)~group_list)$p.value
})
#p值矯正,如何選擇method?
p.adj = p.adjust(pvals, method = "BH")
avg_1 = rowMeans(dat1)
avg_2 = rowMeans(dat2)
log2FC = avg_2-avg_1
DEG_t.test = cbind(avg_1, avg_2, log2FC, pvals, p.adj)
DEG_t.test=DEG_t.test[order(DEG_t.test[,4]),]
#為什么轉(zhuǎn)換成dataframe?
DEG_t.test=as.data.frame(DEG_t.test)

#19 DEG by limma 
suppressMessages(library(limma)) 
##分組矩陣:design
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(exprSet)
design
##比較矩陣:contrast(制定了誰比誰這個(gè)規(guī)則,這個(gè)矩陣聲明走哺，我們要把progres.組跟stable進(jìn)行差異分析比較)
contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)
contrast.matrix 
##limma DEG steps
fit <- lmFit(exprSet,design)
fit2 <- contrasts.fit(fit, contrast.matrix) ##這一步很重要，大家可以自行看看效果
fit2 <- eBayes(fit2,trend = TRUE)  
## default no trend !!!eBayes() with trend=TRUE
tempOutput = topTable(fit2, coef=1, n=Inf)
nrDEG = na.omit(tempOutput) 
head(nrDEG)
## volcano plot
DEG=nrDEG
logFC_cutoff <- with(DEG,mean(abs( logFC)) + 2*sd(abs( logFC)) )
DEG$change = as.factor(ifelse(DEG$P.Value < 0.05 & abs(DEG$logFC) > logFC_cutoff,
                              ifelse(DEG$logFC > logFC_cutoff ,'UP','DOWN'),'NOT')
)
this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),
                    '\nThe number of up gene is ',nrow(DEG[DEG$change =='UP',]) ,
                    '\nThe number of down gene is ',nrow(DEG[DEG$change =='DOWN',])
)
this_tile
head(DEG)
g = ggplot(data=DEG, aes(x=logFC, y=-log10(P.Value), color=change)) +
  geom_point(alpha=0.4, size=1.75) +
  theme_set(theme_set(theme_bw(base_size=20)))+
  xlab("log2 fold change") + ylab("-log10 p-value") +
  ggtitle( this_tile  ) + theme(plot.title = element_text(size=15,hjust = 0.5))+
  scale_colour_manual(values = c('blue','black','red'))  ## corresponding to the levels(res$change)
print(g)
volcano
#20 對T檢驗(yàn)結(jié)果的P值和limma包差異分析的P值畫散點(diǎn)圖哲虾，看看哪些基因相差很大丙躏。
DEG_t.test=DEG_t.test[rownames(nrDEG),]
#logFC
plot(DEG_t.test[,3],nrDEG[,1]) 
#P.value
plot(DEG_t.test[,4],nrDEG[,4]) 
plot(-log10(DEG_t.test[,4]),-log10(nrDEG[,4]))
圖片.png

圖片.png

圖片.png
library(ggplot2)
library(ggpubr)
my_comparisons <- list(c("stable", "progres."))
dat=data.frame(group=group_list,
               sampleID= names(exprSet['DLEU1',]),
               values= as.numeric(exprSet['DLEU1',]))
ggboxplot(
  dat, x = "group", y = "values",
  color = "group",
  add = "jitter"
)+
  stat_compare_means(comparisons = my_comparisons, method = "t.test")
圖片.png
## 前25個(gè)基因的heatmap 
library(pheatmap)
choose_gene=head(rownames(nrDEG),25)
choose_matrix=exprSet[choose_gene,]
choose_matrix=t(scale(t(choose_matrix)))
pheatmap(choose_matrix)
圖片.png
最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 人面猴
    序言：七十年代末，一起剝皮案震驚了整個(gè)濱河市束凑，隨后出現(xiàn)的幾起案子晒旅，更是在濱河造成了極大的恐慌，老刑警劉巖汪诉，帶你破解...
    沈念sama閱讀 212,816評論 6贊 492
  • 死咒
    序言：濱河連續(xù)發(fā)生了三起死亡事件废恋，死亡現(xiàn)場離奇詭異，居然都是意外死亡扒寄，警方通過查閱死者的電腦和手機(jī)鱼鼓，發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 90,729評論 3贊 385
  • 救了他兩次的神仙讓他今天三更去死
    文/潘曉璐 我一進(jìn)店門，熙熙樓的掌柜王于貴愁眉苦臉地迎上來该编，“玉大人迄本，你說我怎么就攤上這事】慰ⅲ” “怎么了嘉赎？”我有些...
    開封第一講書人閱讀 158,300評論 0贊 348
  • 道士緝兇錄：失蹤的賣姜人
    文/不壞的土叔 我叫張陵，是天一觀的道長于樟。 經(jīng)常有香客問我曹阔，道長，這世上最難降的妖魔是什么隔披？ 我笑而不...
    開封第一講書人閱讀 56,780評論 1贊 285
  • ?港島之戀（遺憾婚禮）
    正文 為了忘掉前任赃份，我火速辦了婚禮，結(jié)果婚禮上，老公的妹妹穿的比我還像新娘抓韩。我一直安慰自己纠永，他們只是感情好，可當(dāng)我...
    茶點(diǎn)故事閱讀 65,890評論 6贊 385
  • 惡毒庶女頂嫁案：這布局不是一般人想出來的
    文/花漫 我一把揭開白布谒拴。 她就那樣靜靜地躺著尝江，像睡著了一般。 火紅的嫁衣襯著肌膚如雪英上。 梳的紋絲不亂的頭發(fā)上炭序，一...
    開封第一講書人閱讀 50,084評論 1贊 291
  • 城市分裂傳說
    那天，我揣著相機(jī)與錄音苍日，去河邊找鬼惭聂。 笑死，一個(gè)胖子當(dāng)著我的面吹牛相恃，可吹牛的內(nèi)容都是我干的辜纲。 我是一名探鬼主播，決...
    沈念sama閱讀 39,151評論 3贊 410
  • 雙鴛鴦連環(huán)套：你想象不到人心有多黑
    文/蒼蘭香墨 我猛地睜開眼拦耐，長吁一口氣：“原來是場噩夢啊……” “哼耕腾！你這毒婦竟也來了？” 一聲冷哼從身側(cè)響起杀糯，我...
    開封第一講書人閱讀 37,912評論 0贊 268
  • 萬榮殺人案實(shí)錄
    序言：老撾萬榮一對情侶失蹤扫俺，失蹤者是張志新（化名）和其女友劉穎，沒想到半個(gè)月后固翰，有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體牵舵，經(jīng)...
    沈念sama閱讀 44,355評論 1贊 303
  • ?護(hù)林員之死
    正文 獨(dú)居荒郊野嶺守林人離奇死亡，尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點(diǎn)故事閱讀 36,666評論 2贊 327
  • ?白月光啟示錄
    正文 我和宋清朗相戀三年倦挂，在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了畸颅。 大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片。...
    茶點(diǎn)故事閱讀 38,809評論 1贊 341
  • 活死人
    序言：一個(gè)原本活蹦亂跳的男人離奇死亡方援，死狀恐怖没炒，靈堂內(nèi)的尸體忽然破棺而出，到底是詐尸還是另有隱情犯戏，我是刑警寧澤送火，帶...
    沈念sama閱讀 34,504評論 4贊 334
  • ?日本核電站爆炸內(nèi)幕
    正文 年R本政府宣布，位于F島的核電站先匪，受9級特大地震影響种吸，放射性物質(zhì)發(fā)生泄漏。R本人自食惡果不足惜呀非，卻給世界環(huán)境...
    茶點(diǎn)故事閱讀 40,150評論 3贊 317
  • 男人毒藥：我在死后第九天來索命
    文/蒙蒙 一坚俗、第九天 我趴在偏房一處隱蔽的房頂上張望镜盯。 院中可真熱鬧，春花似錦猖败、人聲如沸速缆。這莊子的主人今日做“春日...
    開封第一講書人閱讀 30,882評論 0贊 21
  • 一樁弒父案恩闻，背后竟有這般陰謀
    文/蒼蘭香墨 我抬頭看了看天上的太陽艺糜。三九已至，卻和暖如春幢尚，著一層夾襖步出監(jiān)牢的瞬間破停，已是汗流浹背。 一陣腳步聲響...
    開封第一講書人閱讀 32,121評論 1贊 267
  • 情欲美人皮
    我被黑心中介騙來泰國打工尉剩， 沒想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留真慢，地道東北人。 一個(gè)月前我還...
    沈念sama閱讀 46,628評論 2贊 362
  • 代替公主和親
    正文 我出身青樓边涕，卻偏偏與公主長得像，于是被迫代替她去往敵國和親褂微。 傳聞我的和親對象是個(gè)殘疾皇子功蜓，可洞房花燭夜當(dāng)晚...
    茶點(diǎn)故事閱讀 43,724評論 2贊 351