https://www.bilibili.com/video/av25643438
嗶哩嗶哩，劈里啪啦，1.5倍速飛快聽完，然后主要整理了R語言最新十道題，壓力就是第一生產(chǎn)力啊。http://www.bio-info-trainee.com/3750.html

準(zhǔn)備工作--安裝所需的包

cran_packages <- c('tidyverse',
                   'ggpubr',
                   'ggstatsplot')
Biocductor_packages <- c('org.Hs.eg.db',
                         'hgu133a.db',
                         'CLL',
                         'hgu95av2.db',
                         'survminer',
                         'survival',
                         'hugene10sttranscriptcluster',
                         'limma')
if(length(getOption("CRAN"))==0) options(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")
for (pkg in cran_packages){
  if (! require(pkg,character.only=T) ) {
    install.packages(pkg,ask = F,update = F)
    require(pkg,character.only=T) 
  }
}
# first prepare BioManager on CRAN
if(length(getOption("CRAN"))==0) options(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")
if(!require("BiocManager")) install.packages("BiocManager",update = F,ask = F)
if(length(getOption("BioC_mirror"))==0) options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
# use BiocManager to install
for (pkg in Biocductor_packages){
  if (! require(pkg,character.only=T) ) {
    BiocManager::install(pkg,ask = F,update = F)
    require(pkg,character.only=T) 
  }
}

1.找到指定ensembl ID與symbol的對應(yīng)關(guān)系

ENSG00000000003.13
ENSG00000000005.5
ENSG00000000419.11
ENSG00000000457.12
ENSG00000000460.15
ENSG00000000938.11

思路：在注釋包中有g(shù)ene_id與symbol、gene_id與ensembl_id的對應(yīng)關(guān)系。

#將以上基因id保存在a.txt，存放于工作目錄下。
rm(list=ls())
options(stringsAsFactors = F)

a=read.table('e1.txt')

g2s <- toTable(org.Hs.egSYMBOL)
g2e <- toTable(org.Hs.egENSEMBL)
a$V1 = apply(a[1], 1,function(x){
  str_split(x,'[.]')[[1]][1]
}) %>% unlist
colnames(a) <- 'ensembl_id'
tmp <- merge(a,g2e, by="ensembl_id")
result <- merge(tmp,g2s, by="gene_id")[-1]

2.根據(jù)探針名找對應(yīng)symbol ID

1053_at
117_at
121_at
1255_g_at
1316_at
1320_at
1405_i_at
1431_at
1438_at
1487_at
1494_f_at
1598_g_at
160020_at
1729_at
177_at

思路：找到注釋包中探針與symbol的對應(yīng)關(guān)系然后取子集

rm(list=ls())
options(stringsAsFactors = F)
library(hgu133a.db)
p2s=toTable(hgu133aSYMBOL)
a=read.table('e2.txt')
colnames(a) <- colnames(p2s)[1]
# 方法一：利用merge
tmp1 <- merge(a,p2s, by="probe_id")
# 方法二：利用match得到第一組向量在第二組中的坐標(biāo)
tmp2 <- p2s[match(a$probe_id,p2s$probe_id),]
## 附：判斷得到的兩組結(jié)果是否一致
# 法一：
identical(tmp1,tmp2) #返回邏輯值
# 法二：
dplyr::setdiff(tmp1,tmp2) #返回兩組的差別【沒差就返回空】

3.根據(jù)symbol找基因表達(dá)量并作圖

找到R包CLL內(nèi)置的數(shù)據(jù)集的表達(dá)矩陣?yán)锩娴腡P53基因的表達(dá)量，并且繪制在 progres.-stable分組的boxplot圖,并通過 ggpubr 進(jìn)行美化。
探針的三大內(nèi)容：表達(dá)矩陣assay/exprs、探針信息featureData、樣本信息phenoData

# 從內(nèi)置數(shù)據(jù)集的表達(dá)矩陣中找TP53基因的表達(dá)量
rm(list=ls())
options(stringsAsFactors = F)
suppressMessages(library(CLL))
data(sCLLex)
# sCLLex
exprSet <- exprs(sCLLex) #探針的表達(dá)量
pdata <- pData(sCLLex) #sampleID與disease的對應(yīng)關(guān)系
p2s <- toTable(hgu95av2SYMBOL) #探針與symbol的對應(yīng)關(guān)系
p3 <- filter(p2s,symbol=='TP53')
# boxplot [find TP53 has 3 probe IDs]
probe_tp53 <- c("1939_at","1974_s_at","31618_at")
i = 3 #可換1，2
boxplot(exprSet[probe_tp53[i],] ~ pdata$Disease)

#用ggpubr作圖
#http://www.sthda.com/english/articles/24-ggpubr-publication-ready-plots/
exp_tab <- rownames_to_column(as.data.frame(exprSet))
exp_tab2 <- gather(exp_tab,
                   key = 'sample',
                   value = 'exp',-1)
pdata <- rownames_to_column(pdata)
exp_tab3 <- merge(exp_tab2,pdata,by.x='sample',by.y='rowname')
i=1 ###可換1，2
dev.off()
p <- ggboxplot(filter(exp_tab3,rowname==probe_tp53[i]), 
               x = 'Disease',
               y = 'exp',
               color = "Disease", palette =c("#00AFBB", "#E7B800", "#FC4E07"),
               add = "jitter", shape = "Disease")
p

4.BRCA1基因表達(dá)量

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

提示：使用http://www.cbioportal.org/index.do 定位數(shù)據(jù)集：http://www.cbioportal.org/datasets
該基因有四個亞型，用ggbetweenstats作圖比較一下。

rm(list=ls())
options(stringsAsFactors = F)
#ID，四個亞型，表達(dá)量
f <- read.csv("e4-plot.txt", sep = "\t")
## boxplot
colnames(f) <- c("id", "subtype", "expression", "mut")
da <- f
library(ggstatsplot)
ggbetweenstats(data = da, 
               x = subtype,
               y = expression)
library(ggplot2)
ggsave("e4-BRCA1-TCGA.png")

5 .TP53 生存分析

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

提示使用：http://www.oncolnc.org/
思路：生存分析，TP53表達(dá)量分為高低兩組做圖比較

# Use http://www.oncolnc.org/ to get raw csv da

rm(list=ls())
options(stringsAsFactors = F)

f <- read.csv('e5-BRCA_7157_50_50.csv')
library(ggstatsplot)
ggbetweenstats(data = da,
               x = Group,
               y = Expression)
da <- f

library(ggplot2)
library(survival)
library(survminer)
table(da$Status)

da$Status <- ifelse(da$Status == "Dead", 1, 0)
survf <- survfit(Surv(Days,Status)~Group, data=da)
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")

6.從表達(dá)矩陣中提取指定基因畫熱圖

下載數(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

提示：根據(jù)基因名拿到探針I(yè)D，縮小表達(dá)矩陣?yán)L制熱圖，沒有檢查到的基因直接忽略即可。

## Exercise 5: Retrive genes from GEO to plot heatmap

rm(list=ls())
options(stringsAsFactors = F)
#下載和表達(dá)矩陣
library(GEOquery)
GSE <- "GSE17215"
if(!file.exists(GSE)){
    geo <- getGEO(GSE, destdir = '.', getGPL = F, AnnotGPL = F)
    save(geo, file = paste0(GSE,'.eSet.Rdata'))
}
load(paste0(GSE,'.eSet.Rdata'))
expr <- exprs(geo[[1]])
p2s=toTable(hgu133aSYMBOL);head(p2s)
expr <- expr[p2s$probe_id,] #有的id找不到注釋直接刪掉
gp <- "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"
gp <- strsplit(gp, ' ')[[1]]
tmp <- dplyr::filter(p2s, p2s$symbol %in% gp)
tmp2 <- tibble::rownames_to_column(data.frame(expr),"probe_id")
tmp3 <- merge(tmp,tmp2,by="probe_id")

tmp3$median <- apply(tmp3[,3:ncol(tmp3)], 1, median)
tmp3 <- tmp3[order(tmp3$symbol, tmp3$median, decreasing = T),]
tmp3 <- tmp3[!duplicated(tmp3$symbol),]
rownames(tmp3) <- tmp3$symbol
tmp3 <- tmp3[,-c(1,2,ncol(tmp3))]

gp_expr <- log2(tmp3)
pheatmap::pheatmap(gp_expr, scale = "row") 

7.相關(guān)性計(jì)算和熱圖

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

rm(list=ls())
options(stringsAsFactors = F)

library(GEOquery)
GSE <- "GSE24673"
if(!file.exists(GSE)){
  geo <- getGEO(GSE, destdir = '.', getGPL = F, AnnotGPL = F)
  save(geo, file = paste0(GSE,'.eSet.Rdata'))
}
load(paste0(GSE,'.eSet.Rdata'))
expr <- exprs(geo[[1]])
dim(expr)
expr[1:4,1:4]
pdata <- pData(geo[[1]]) 
# 自己根據(jù)pdata第八列做一個分組信息矩陣
grp <- c('rbc','rbc','rbc',
         'rbn','rbn','rbn',
         'rbc','rbc','rbc',
         'normal','normal')

grp_df <- data.frame(group=grp)
rownames(grp_df) <- colnames(expr)
new_cor <- cor(expr)
pheatmap::pheatmap(new_cor, annotation_col = grp_df)

8.找到芯片平臺對應(yīng)的注釋包

找到 GPL6244 platform of Affymetrix Human Gene 1.0 ST Array 對應(yīng)的R的bioconductor注釋包，并且安裝它！
https://mp.weixin.qq.com/s/sVSsV2fWZOQmNd72Vb3SmQ
http://m.itdecent.cn/p/40b560755cdf

options()$repos
options()$BioC_mirror 
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
BiocManager::install("hugene10sttranscriptcluster",ask = F,update = F)
options()$repos
options()$BioC_mirror

9.找到指定探針和對應(yīng)的基因

rm(list=ls())
options(stringsAsFactors = F)

library(GEOquery)
GSE <- "GSE42872"
if(!file.exists(GSE)){
  geo <- getGEO(GSE, destdir = '.', getGPL = F, AnnotGPL = F)
  save(geo, file = paste0(GSE,'.eSet.Rdata'))
}
load(paste0(GSE,'.eSet.Rdata'))
expr <- exprs(geo[[1]])

dim(expr)
expr[1:4,1:4]

# 選出所有樣本的(mean/sd/mad/)最大的探針

sort(apply(expr,1,mean),decreasing = T)[1]
sort(apply(expr,1,sd),decreasing = T)[1]
sort(apply(expr,1,mad),decreasing = T)[1]

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

10.limma 差異分析

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

# 接第九題，得到表型信息，然后用limma做差異分析
pdata <- pData(geo[[1]]) 
grp <- unlist(lapply(pdata$title, function(x){
  strsplit(x, ' ')[[1]][4]
}))

suppressMessages(library(limma))
#limma needs：表達(dá)矩陣（expr：需要用log后的矩陣）、分組矩陣（design）、比較矩陣（contrast）
#先做一個分組矩陣～design，說明progres是哪幾個樣本，stable又是哪幾個
design <- model.matrix(~0+factor(grp))
colnames(design) <- levels(factor(grp))
rownames(design) <- colnames(expr)
design
#再做一個比較矩陣【一般是case比control】
contrast<-makeContrasts(paste0(unique(grp),collapse = "-"),levels = design)
contrast

##step1
fit <- lmFit(expr,design)
##step2
fit2 <- contrasts.fit(fit, contrast) 
fit2 <- eBayes(fit2)  
##step3
mtx = topTable(fit2, coef=1, n=Inf)
DEG_mtx = na.omit(mtx) 
View(DEG_mtx)

# 火山圖

DEG=DEG_mtx

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)
vol1 = 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=15,hjust = 0.5))+
  scale_colour_manual(values = c('blue','black','red')) # according to the levels(DEG$change)
print(vol1)

微信公眾號生信星球同步更新我的文章，歡迎大家掃碼關(guān)注！