R語(yǔ)言小作業(yè)-中級(jí)

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

新建一個(gè)e1.txt的文檔

rm(list=ls())
options(stringsAsFactors = F)#如果不加后面會(huì)在strsplit(x,'[.]')時(shí)提示：non-character argument.
a=read.table('e1.txt')
head(a)
library(org.Hs.eg.db)
ls("package:org.Hs.eg.db") #圖1
g2s <- toTable(org.Hs.egSYMBOL) #圖2
g2e <- toTable(org.Hs.egENSEMBL) #圖2
library(stringr)#只要"."之前的ensembl ID
#或者用lapply函數(shù)
lapply(a$V1,function(x){
  strsplit(as.character(x),'[.]')[[1]][1]#以"."為分割
  })#在沒(méi)加options(stringsAsFactors = F的情況下加as.character一樣
unlist(lapply(a$V1,function(x){
+   strsplit((x),'[.]')[[1]][1]
+ })) #圖3
a$ensembl_id=unlist(lapply(a$V1,function(x){
  strsplit((x),'[.]')[[1]][1]
}))
tmp <- merge(a,g2e, by="ensembl_id")#圖4
tmp <- merge(tmp,g2s, by="gene_id")#圖5

圖1

圖2

圖3

圖4

圖5

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)
a=read.table('e2.txt')#此時(shí)a是新建的e2.txt
colnames(a)='probe_id' #圖7 a的列名為V1，而與ids的列名不同，因此把a(bǔ)的列名改為與ids相同的'probe_id'，用于后面merge。
library(hgu133a.db)
ls("package:hgu133a.db")#圖6
ids=toTable(hgu133aSYMBOL)#圖7
head(ids)
tmp1=merge(ids,a,by='probe_id')#圖7
tmp2=ids[match(a$'probe_id',ids$'probe_id'),]#圖7 用match函數(shù)得到前一個(gè)向量在后一個(gè)向量的坐標(biāo),再在ids中取出來(lái)，','左邊即按照這些坐標(biāo)把這坐標(biāo)對(duì)應(yīng)的這行取出來(lái)
## 附：判斷得到的兩組結(jié)果是否一致
# 法一：
identical(tmp1,tmp2) #返回邏輯值
# 法二：
dplyr::setdiff(tmp1,tmp2) #返回兩組的差別【沒(méi)差就返回空】

圖6

圖7

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

rm(list=ls())
options(stringsAsFactors = F)
suppressPackageStartupMessages(library(CLL))
data(sCLLex)
sCLLex
exprSet=exprs(sCLLex) #可通過(guò)exprs來(lái)獲取這個(gè)對(duì)象的表達(dá)矩陣exprSet
# sCLLex
exprSet <- exprs(sCLLex) #探針的表達(dá)量
pd<- pData(sCLLex) #sampleID與disease的對(duì)應(yīng)關(guān)系-分組信息
p2s <- toTable(hgu95av2SYMBOL) #探針與symbol的對(duì)應(yīng)關(guān)系
p3 <- filter(p2s,symbol=='TP53')
library(hgu95av2.db)
ids=toTable(hgu95av2SYMBOL)#圖8
head(ids)
boxplot(exprSet['1939_at',]~pd$Disease)#圖9
boxplot(exprSet['1974_s_at',]~pd$Disease)#圖9
boxplot(exprSet['31618_at',]~pd$Disease)#圖9

##接下來(lái)這段摘自小潔（http://www.itdecent.cn/p/1766eae6f97a），留下來(lái)
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

圖8-發(fā)現(xiàn)TP53有3個(gè)probe_id

圖9

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見(jiàn)圖10

圖10

#下載數(shù)據(jù)后自己做一遍
rm(list=ls())
options(stringsAsFactors = F)
#ID，四個(gè)亞型，表達(dá)量
a= read.table("e4-plot.txt", sep = "\t",fill=T,header=T)
## boxplot
colnames(a) <- c("id", "subtype", "expression", "mut")
dat=a
library(ggstatsplot)
ggbetweenstats(data = dat, 
               x = subtype,
               y = expression)
library(ggplot2)
ggsave("e4-BRCA1-TCGA.png")#圖11ggsave保存需要下載ggplot2這個(gè)包

圖11

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

提示使用：http://www.oncolnc.org/見(jiàn)圖12、圖13

圖12

圖13

rm(list = ls())
options(stringsAsFactors = F)
a=read.table('BRCA_7157_50_50.csv',sep = ',',fill = T,header = T)
dat=a
library(ggplot2)
library(survival)
library(survminer) 
table(dat$Status)
dat$Status=ifelse(dat$Status=='Dead',1,0)#圖14 將status改名，Dead改為1，Alive就改為0
sfit <- survfit(Surv(Days, Status)~Group, data=dat)
sfit
summary(sfit)
ggsurvplot(sfit, conf.int=F, pval=TRUE)
ggsave('survival_TP53_in_BRCA_TCGA.png')#圖15

圖14

圖15

head(a)
b=read.table('e4-plot.txt',sep = '\t',fill = T,header = T)
colnames(b)=c('Patient','subtype','expression','mut')#圖17
head(b)#此處的b其實(shí)是上一題中的a，為了與本題中的a的Patient ID相一致，此時(shí)把b的第一列的名字命名為Patient.
b$Patient=substring(b$Patient,1,12)#圖16記住substring取字符用法
tmp=merge(a,b,by='Patient')#圖17

圖16

圖17

可以檢查下merge后的tmp肯定會(huì)有減少，這時(shí)需要注意下看看還有多少個(gè)觀測(cè)，如圖顯示還有800多，還不錯(cuò)。見(jiàn)圖18

圖18

#如果不寫循環(huán)
table(tmp$subtype)

type='BRCA_LumB'
x=tmp[tmp$subtype==type,] 
library(ggplot2)
library(survival)
library(survminer) 
#table(x$Status)
x$Status=ifelse(x$Status=='Dead',1,0)
sfit <- survfit(Surv(Days, Status)~Group, data=x)
sfit
summary(sfit)
ggsurvplot(sfit, conf.int=F, pval=TRUE)  #圖19

圖19

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)有檢查到的基因直接忽略即可。

##下載表達(dá)矩陣
rm(list = ls())  ## 魔幻操作，一鍵清空~
options(stringsAsFactors = F)
# 注意查看下載文件的大小，檢查數(shù)據(jù) 
f='GSE17215_eSet.Rdata'

library(GEOquery)#加載GEO這個(gè)包
# 這個(gè)包需要注意兩個(gè)配置，一般來(lái)說(shuō)自動(dòng)化的配置是足夠的。
#Setting options('download.file.method.GEOquery'='auto')
#Setting options('GEOquery.inmemory.gpl'=FALSE)
if(!file.exists(f)){
  gset <- getGEO('GSE17215', destdir=".",
                 AnnotGPL = F,     ## 注釋文件
                 getGPL = F)       ## 平臺(tái)文件
  save(gset,file=f)   ## 保存到本地
}
load('GSE17215_eSet.Rdata')  ## 載入數(shù)據(jù)
class(gset)
length(gset)
class(gset[[1]])
# 因?yàn)檫@個(gè)GEO數(shù)據(jù)集只有一個(gè)GPL平臺(tái)，所以下載到的是一個(gè)含有一個(gè)元素的list
a=gset[[1]]
dat=exprs(a)#圖21-1-1（圖21的第1個(gè)圖）
dim(dat)


#再進(jìn)一步過(guò)濾基因hgu133a.db這個(gè)包
library(hgu133a.db)#圖20
ids=toTable(hgu133aSYMBOL)
head(ids)
dat=dat[ids$probe_id,]#圖21-2（圖21的第2個(gè)圖）
dat[1:4,1:4] 
ids$median=apply(dat,1,median)#圖22-1.對(duì)dat每一行求median值
ids=ids[order(ids$symbol,ids$median,decreasing = T),]#圖22-2.order是按從小到達(dá)排序，加了decreasing = T是按從大到小排序。ids$symbol是將symbol基因名按照字母順序從最后即"z"排序的，然后每個(gè)z開(kāi)頭的再按ids$median排序
ids=ids[!duplicated(ids$symbol),]#圖22-3.將ids$symbol去重后作為新的行，ids取這些行
dat=dat[ids$probe_id,]#圖23-3和圖23-2.dat這個(gè)表達(dá)矩陣取出新的ids的probi_id這整行為一個(gè)新dat
rownames(dat)=ids$symbol#圖23-2
dat[1:4,1:4]  
dim(dat)

圖20

圖21

圖22

圖23

# 想要的基因
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'
ng=strsplit(ng,' ')[[1]]#圖24
ng
table(ng %in%  rownames(dat))#看我們想要的基因有哪些在dat（GSE17215）這個(gè)基因集里
ng=ng[ng %in%  rownames(dat)]
dat=dat[ng,]
dat=log2(dat)
pheatmap::pheatmap(dat,scale = 'row')#圖25.把row用scale歸一化

圖24

圖25

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

rm(list = ls())  ## 魔幻操作，一鍵清空~
options(stringsAsFactors = F)
# 注意查看下載文件的大小，檢查數(shù)據(jù) 
f='GSE24673_eSet.Rdata'
library(GEOquery)
if(!file.exists(f)){
  gset <- getGEO('GSE24673', destdir=".",
                 AnnotGPL = F,     
                 getGPL = F)       
  save(gset,file=f)  
}
load('GSE24673_eSet.Rdata')  ## 載入數(shù)據(jù)
class(gset)
length(gset)
class(gset[[1]])
# 因?yàn)檫@個(gè)GEO數(shù)據(jù)集只有一個(gè)GPL平臺(tái)，所以下載到的是一個(gè)含有一個(gè)元素的list
a=gset[[1]]
dat=exprs(a)#圖26.用exprs取出表達(dá)矩陣
dim(dat)
pd=pData(a)#圖27.用pData取出臨床信息，即可以得到分組信息
group_list=c('rbc','rbc','rbc', #看了pd后創(chuàng)建分組信息
             'rbn','rbn','rbn',
             'rbc','rbc','rbc',
             'normal','normal')
dat[1:4,1:4]
M=cor(dat)#圖28.看樣本和樣本之間的相關(guān)性
pheatmap::pheatmap(M)
tmp=data.frame(g=group_list)#圖29
rownames(tmp)=colnames(M)
pheatmap::pheatmap(M,annotation_col = tmp)#圖28 分組信息用annotation_col.

圖26

圖27

圖28

圖29

圖30

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

#提示
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("請(qǐng)輸入自己找到的R包",ask = F,update = F)
options()$repos
options()$BioC_mirror

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

####依然同上一題一樣 下載矩陣
rm(list = ls()) 
options(stringsAsFactors = F)

f='GSE42872_eSet.Rdata' 
library(GEOquery)
if(!file.exists(f)){
  gset <- getGEO('GSE42872', destdir=".",
                 AnnotGPL = F,    
                 getGPL = F)       
  save(gset,file=f)  
}
load('GSE42872_eSet.Rdata')  ## 載入數(shù)據(jù)
class(gset)
length(gset)
class(gset[[1]])
# 因?yàn)檫@個(gè)GEO數(shù)據(jù)集只有一個(gè)GPL平臺(tái)，所以下載到的是一個(gè)含有一個(gè)元素的list
a=gset[[1]]
dat=exprs(a)#用exprs取出表達(dá)矩陣
dim(dat)
pd=pData(a)#用pData取出臨床信息，即可以得到分組信息
# (平均表達(dá)量/sd/mad/)最大的探針
boxplot(dat)
sort(apply(dat,1,mean),decreasing = T)[1]
sort(apply(dat,1,sd),decreasing = T)[1]
sort(apply(dat,1,mad),decreasing = T)[1]

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)     
  save(gset,file=f)   
}
load('GSE42872_eSet.Rdata')  ## 載入數(shù)據(jù)
class(gset)
length(gset)
class(gset[[1]])
# 因?yàn)檫@個(gè)GEO數(shù)據(jù)集只有一個(gè)GPL平臺(tái)，所以下載到的是一個(gè)含有一個(gè)元素的list
a=gset[[1]]
dat=exprs(a)
dim(dat)
pd=pData(a)#圖31
# (平均表達(dá)量/sd/mad/)最大的探針
boxplot(dat)
group_list=unlist(lapply(pd$title,function(x)
  strsplit(x,' ')[[1]][4]
))


exprSet=dat#目前的矩陣是dat，轉(zhuǎn)換成limma包里代碼用的exprSet一致
exprSet[1:4,1:4]
# DEG by limma 
suppressMessages(library(limma)) 
design <- model.matrix(~0+factor(group_list))#圖31
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(exprSet)
design
contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)
contrast.matrix<-makeContrasts("progres.-stable",levels = design)
contrast.matrix 
##這個(gè)矩陣聲明，我們要把progres.組跟stable進(jìn)行差異分析比較
##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)

圖31