首先需要完成R語(yǔ)言練習(xí)題-初級(jí)胜蛉，在http://www.bio-info-trainee.com/3793.html
作業(yè) 1
請(qǐng)根據(jù)R包org.Hs.eg.db找到下面ensembl 基因ID 對(duì)應(yīng)的基因名(symbol)
ENSG00000000003.13
ENSG00000000005.5
ENSG00000000419.11
ENSG00000000457.12
ENSG00000000460.15
ENSG00000000938.11
解：
# 先將待解的ensembl基因ID生成txt，為e1.txt
rm(list = ls())
options(stringsAsFactors = F)
a=read.table('e1.txt')
head(a)
# 載入注釋包（需提前安裝庞呕，見(jiàn)下）
library(org.Hs.eg.db)
# 查看包內(nèi)容
ls("package:org.Hs.eg.db")
g2s=toTable(org.Hs.egSYMBOL);head(g2s)
g2e=toTable(org.Hs.egENSEMBL);head(g2e)
head(g2e)
library(stringr)
# 用strsplit對(duì)a進(jìn)行分割，取"."之前的內(nèi)容肴敛，即ensembl_id
a$ensembl_id=unlist(lapply(a$V1,function(x){
strsplit(as.character(x),'[.]')[[1]][1]
})
)
# 將三者關(guān)聯(lián)
tmp1=merge(a,g2e,by='ensembl_id')
tmp2=merge(tmp1,g2s,by='gene_id')
knitr::opts_chunk$set(echo = TRUE)
image.png
作業(yè) 2
根據(jù)R包hgu133a.db找到下面探針對(duì)應(yīng)的基因名(symbol)
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
解：
rm(list = ls())
options(stringsAsFactors = F)
# 將探針保存為txt归园，命名為e2停团，并讀取
a=read.table('e2.txt')
# 修改列名
colnames(a)='probe_id'
library(hgu133a.db)
# 找到注釋文件內(nèi)探針與基因的對(duì)應(yīng)信息
ids=toTable(hgu133aSYMBOL)
head(ids)
# 合并
tmp1=merge(ids,a,by='probe_id')
# tmp1和tmp2效果相同
tmp2=ids[match(a$probe_id,ids$probe_id),]
# 如果之前不修改a的列名，可用
tmp3=merge(ids,a,by.x='probe_id',by.y="V1")
image.png
作業(yè) 3
找到R包CLL內(nèi)置的數(shù)據(jù)集的表達(dá)矩陣?yán)锩娴腡P53基因的表達(dá)量袭艺，并且繪制在 progres.-stable分組的boxplot圖搀崭。想想如何通過(guò) ggpubr 進(jìn)行美化。
解：
rm(list = ls())
options(stringsAsFactors = F)
# 不提示加載信息
suppressPackageStartupMessages(library(CLL))
# 載入CLL信息
data(sCLLex)
sCLLex
# 獲取表達(dá)矩陣及樣本信息
exprSet=exprs(sCLLex)
pd=pData(sCLLex)
# 加載注釋包
library(hgu95av2.db)
ids=toTable(hgu95av2SYMBOL)
head(ids)
# 在ids直接搜索TP53猾编，得三個(gè)探針瘤睹；或者用代碼獲取（探針數(shù)量多時(shí)推薦）
aaa <- ids[ids$symbol%in%"TP53",][,1]
#畫圖
boxplot(exprSet['1939_at',] ~ pd$Disease)
boxplot(exprSet['1974_s_at',] ~ pd$Disease)
boxplot(exprSet['31618_at',] ~ pd$Disease)
作業(yè) 4
找到BRCA1基因在TCGA數(shù)據(jù)庫(kù)的乳腺癌數(shù)據(jù)集(Breast Invasive Carcinoma (TCGA, PanCancer Atlas))的表達(dá)情況
提示：使用http://www.cbioportal.org/index.do 定位數(shù)據(jù)集：http://www.cbioportal.org/datasets
解：參考
http://www.reibang.com/p/d24a47298a14
作業(yè) 5
找到TP53基因在TCGA數(shù)據(jù)庫(kù)的乳腺癌數(shù)據(jù)集的表達(dá)量分組看其是否影響生存
解：參考
http://www.reibang.com/p/aa727a67948b
作業(yè)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
提示：根據(jù)基因名拿到探針I(yè)D答倡，縮小表達(dá)矩陣?yán)L制熱圖轰传，沒(méi)有檢查到的基因直接忽略即可。
解：
rm(list = ls())
options(stringsAsFactors = F)
# 注意查看下載文件的大小瘪撇，檢查數(shù)據(jù)
f='GSE17215_eSet.Rdata'
# 使用GEOquery包下載GEO數(shù)據(jù)（需提前下載安裝此包）
library(GEOquery)
# 下載可能會(huì)出錯(cuò)获茬，可以試試更改鏡像、翻墻倔既，重啟恕曲，rm（list=ls（），據(jù)說(shuō)早上的網(wǎng)好容易下載...看運(yùn)氣了叉存，下2為更改鏡像
options(repos<- c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/") )
options("BioC_mirror"<- "https://mirrors.ustc.edu.cn/bioc/")
if(!file.exists(f)){
gset <- getGEO('GSE17215', destdir=".",
AnnotGPL = F, ## 注釋文件
getGPL = F) ## 平臺(tái)文件
save(gset,file=f) ## 保存到本地
}
# 先保存后載入码俩，養(yǎng)成好習(xí)慣
load('GSE17215_eSet.Rdata')
# 因?yàn)檫@個(gè)GEO數(shù)據(jù)集只有一個(gè)GPL平臺(tái)，所以下載到的是一個(gè)含有一個(gè)元素的list
class(gset)
length(gset)
class(gset[[1]])
# not list易于操作
a=gset[[1]]
# 獲取表達(dá)矩陣
dat=exprs(a)
dim(dat)
# 加載注釋包
library(hgu133a.db)
# 看一下包里有什么
ls("package:hgu133a.db")
# 獲取探針與基因的對(duì)應(yīng)關(guān)系
ids=toTable(hgu133aSYMBOL)
head(ids)
# 篩選ids中包含的probe_id歼捏，將兩個(gè)表行名設(shè)為同序
dat=dat[ids$probe_id,]
dat[1:4,1:4]
# 取dat表達(dá)量的中位數(shù)稿存，生成ids的median列
ids$median=apply(dat,1,median)
# 將ids按照symbol和median排序（降序），可以看一下前后對(duì)比關(guān)系進(jìn)行理解
ids=ids[order(ids$symbol,ids$median,decreasing = T),]
# 去除重復(fù)的symbol值瞳秽，保留第一位瓣履，即表達(dá)量最大的一位（dim一下看看前后對(duì)比）
ids=ids[!duplicated(ids$symbol),]
# 同上
dat=dat[ids$probe_id,]
# 先同序，再直接改為symbol（dim一下瞧瞧兩者行數(shù)）
rownames(dat)=ids$symbol
dat[1:4,1:4]
dim(dat)
# 好练俐，進(jìn)入正題
ng='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'
# 把基因加上“”袖迎，使之與rownames(dat)格式一致
ng=strsplit(ng,' ')[[1]]
# 看一下，多少基因存在于dat中，可得41個(gè)TRUE&9個(gè)FALSE
table(ng %in% rownames(dat))
# 清洗掉不存在的ng燕锥，注意這一步存在排序（連同下一步理解）
ng=ng[ng %in% rownames(dat)]
# 再改行名
dat=dat[ng,]
# 畫圖
dat=log2(dat)
pheatmap::pheatmap(dat,scale = 'row')
image.png
作業(yè)7
下載數(shù)據(jù)集GSE24673的表達(dá)矩陣計(jì)算樣本的相關(guān)性并且繪制熱圖辜贵，需要標(biāo)記上樣本分組信息
解：
rm(list = ls())
options(stringsAsFactors = F)
GSE_name = 'GSE24673'
options( 'download.file.method.GEOquery' = 'libcurl' ) #windows系統(tǒng)
gset <- getGEO( GSE_name, getGPL = F )
# 保存并載入數(shù)據(jù)
save( gset, file = 'GSE24673_gset.Rdata' )
load(file ="GSE24673_gset.Rdata")
a <- gset[[1]]
dat=exprs(a)
dim(dat)
pd=pData(a)
# 我們通過(guò)觀察樣本信息（pd）來(lái)得知11個(gè)樣本的性質(zhì)（分組信息），建立小分組
group_list=c('rbc','rbc','rbc',
'rbn','rbn','rbn',
'rbc','rbc','rbc',
'normal','normal')
dat[1:4,1:4]
# 相關(guān)性分析并作圖
M=cor(dat)
pheatmap::pheatmap(M)
# 標(biāo)記分組信息归形，相關(guān)性分析并作圖
tmp=data.frame(g=group_list)
rownames(tmp)=colnames(M)
pheatmap::pheatmap(M,annotation_col = tmp)
image.png
image.png
作業(yè)8
找到 GPL6244 platform of Affymetrix Human Gene 1.0 ST Array 對(duì)應(yīng)的R的bioconductor注釋包托慨，并且安裝它！
解：
去菜鳥(niǎo)團(tuán)網(wǎng)站搜索 http://www.bio-info-trainee.com/1399.html
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.db",ask = F,update = F)
options()$repos
options()$BioC_mirror
作業(yè)9
下載數(shù)據(jù)集GSE42872的表達(dá)矩陣暇榴，并且分別挑選出 所有樣本的(平均表達(dá)量/sd/mad/)最大的探針厚棵，并且找到它們對(duì)應(yīng)的基因。
解：
# 經(jīng)典步驟前面已經(jīng)解釋過(guò)很多次了
rm(list = ls())
options(stringsAsFactors = F)
f='GSE42872_eSet.Rdata'
library(GEOquery)
if(!file.exists(f)){
gset <- getGEO('GSE42872', destdir=".",
AnnotGPL = F, ## 注釋文件
getGPL = F) ## 平臺(tái)文件
save(gset,file=f)
}
load('GSE42872_eSet.Rdata')
class(gset)
length(gset)
class(gset[[1]])
a=gset[[1]]
dat=exprs(a)
dim(dat)
pd=pData(a)
# 先畫個(gè)箱圖看看
boxplot(dat)
# 分別挑選出 所有樣本的(平均表達(dá)量/sd/mad/)最大的探針蔼紧，分別為7978905婆硬，8133876和8133876
sort(apply(dat,1,mean),decreasing = T)[1]
sort(apply(dat,1,sd),decreasing = T)[1]
sort(apply(dat,1,mad),decreasing = T)[1]
# 加載bioconduct注釋包（包名有沒(méi)有很熟悉？）
library(hugene10sttranscriptcluster.db)
ls("package:hugene10sttranscriptcluster.db")
ggg <-toTable(hugene10sttranscriptclusterSYMBOL)
# 找一下奸例；尷尬的是7978905找不到...
ggg2 <- ggg[ggg$probe_id%in%7978905,]
ggg2 <- ggg[ggg$probe_id%in%8133876,]
作業(yè)10
下載數(shù)據(jù)集GSE42872的表達(dá)矩陣彬犯，并且根據(jù)分組使用limma做差異分析，得到差異結(jié)果矩陣
解：
rm(list = ls())
options(stringsAsFactors = F)
f='GSE42872_eSet.Rdata'
library(GEOquery)
if(!file.exists(f)){
gset <- getGEO('GSE42872', destdir=".",
AnnotGPL = F, ## 注釋文件
getGPL = F) ## 平臺(tái)文件
save(gset,file=f)
}
load('GSE42872_eSet.Rdata')
class(gset)
length(gset)
class(gset[[1]])
a=gset[[1]]
dat=exprs(a)
dim(dat)
pd=pData(a)
boxplot(dat)
# 從樣本信息提取出分組
group_list=unlist(lapply(pd$title,function(x){
strsplit(x,' ')[[1]][4]
}))
exprSet=dat
exprSet[1:4,1:4]
# 用limma包做差異表達(dá)分析
#差異分析重點(diǎn)在于做好表達(dá)矩陣和分組信息哩至，具體原理可以不用理解
suppressMessages(library(limma))
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(exprSet)
design
contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = "-"),levels = 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)
image.png
參考來(lái)源：生信技能樹(shù)
友情鏈接：
課程分享
生信技能樹(shù)全球公益巡講
（https://mp.weixin.qq.com/s/E9ykuIbc-2Ja9HOY0bn_6g）
B站公益74小時(shí)生信工程師教學(xué)視頻合輯
（https://mp.weixin.qq.com/s/IyFK7l_WBAiUgqQi8O7Hxw）
招學(xué)徒：
（https://mp.weixin.qq.com/s/KgbilzXnFjbKKunuw7NVfw）
歡迎關(guān)注公眾號(hào)：青島生信菜鳥(niǎo)團(tuán)
