題目：
初級(jí)：http://www.bio-info-trainee.com/3793.html
中級(jí)：http://www.bio-info-trainee.com/3750.html
高級(jí)： http://www.bio-info-trainee.com/3415.html
basic
#1打開 `Rstudio` 告訴我它的工作目錄。
getwd()
#2新建6個(gè)向量，基于不同的`原子類型`。（重點(diǎn)是字符串，數(shù)值，邏輯值）
R中的原子型向量可以存儲(chǔ)以下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對(duì)前面讀取的 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 對(duì)應(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找到下面探針對(duì)應(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á)量分組看其是否影響生存
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 對(duì)應(yīng)的R的bioconductor注釋包，并且安裝它！
通過GSE序列號(hào)下載數(shù)據(jù)→從得到的list中獲取GPL序列號(hào)→找到對(duì)應(yīng)的bioconda注釋包
https://mp.weixin.qq.com/s/sVSsV2fWZOQmNd72Vb3SmQ
9. 下載數(shù)據(jù)集GSE42872的表達(dá)矩陣，并且分別挑選出 所有樣本的(平均表達(dá)量/sd/mad/)最大的探針，并且找到它們對(duì)應(yīng)的基因。
SD:Standard Deviation
MAD：Median absolute deviation，就是先求出給定數(shù)據(jù)的中位數(shù)(注意并非均值)然后原數(shù)列的每個(gè)值與這個(gè)中位數(shù)求出絕對(duì)差,然后新數(shù)列的中位值就是MAD。平均絕對(duì)誤差可以避免誤差相互抵消的問題，因而可以準(zhǔn)確反映實(shí)際預(yù)測(cè)誤差的大小。
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,]
#沒懂,再一次對(duì)expr和ids匹配？
ids=ids[match(rownames(exprSet),ids$probe_id),]
#9. 只保留expr中所有樣品平均表達(dá)量最大的探針
##第一種方法
###by(data, INDICES, FUN)函數(shù)的典型用法： 是將data數(shù)據(jù)框或矩陣按照INDICES因子水平進(jìn)行分組，然后對(duì)每組應(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í)對(duì)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
#對(duì)表達(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：代表的是對(duì)向量中的每個(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 對(duì)樣品聚類并且繪圖，然后添加樣本的臨床表型數(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 對(duì)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
