[R語言練習(xí)題]R語言中級練習(xí)題

這是生信技能樹論壇R語言的中級測試題，其中大部分是GEO數(shù)據(jù)挖掘和TCGA數(shù)據(jù)庫的一些操作，雖然我目前用不到，但是卻可以通過查看大佬的代碼來提高自己的認(rèn)識水平，也可以對其他領(lǐng)域有個粗淺的認(rèn)識。

1.根據(jù)R包org.Hs.eg.db找到下面ensembl 基因ID 對應(yīng)的基因名(symbol)
ENSG00000000003.13
ENSG00000000005.5
ENSG00000000419.11
ENSG00000000457.12
ENSG00000000460.15
ENSG00000000938.11

BiocManager::install("org.Hs.eg.db")
library(stringr)
library(org.Hs.eg.db)
options(stringsAsFactors = F)
a <- read.csv("s2-1.txt",header = F)
g2s=toTable(org.Hs.egSYMBOL)
g2e=toTable(org.Hs.egENSEMBL)
class(a)
b <- unlist(lapply(a$V1,function(x){
  str_split(x,"[.]")[[1]][1]}))
b <- as.data.frame(b)
c <- merge(b,g2e,by.x="b",by.y="ensembl_id")
d <- merge(c,g2s,by.x="gene_id",by.y="gene_id")
head(d)

2.根據(jù)R包hgu133a.db找到下面探針對應(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

BiocManager::install("hgu133a.db")
library(hgu133a.db)
library(stringr)
ids=toTable(hgu133aSYMBOL)
head(ids)
options(stringsAsFactors = F)
a <- read.csv("s2-2.txt",sep = "\t",header = F)
#方法1，使用merge()
b <- merge(a,ids,by.x="V1",by.y="probe_id")
head(b)
#方法二，使用match(),返回的是前面的向量在后面的坐標(biāo)
c <- ids[match(a$V1,ids$probe_id),]
b==c #判斷bc是否相等

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

options(stringsAsFactors = F)
rm(list=ls())
BiocManager::install("CLL")
suppressPackageStartupMessages(library(CLL))
data(sCLLex)
sCLLex
exprSet=exprs(sCLLex) 
BiocManager::install("hgu95av2.db")
library(hgu95av2.db)
ids=toTable(hgu95av2SYMBOL)

pd <- pData(sCLLex) #查看分組信息
#五個以下不用寫循環(huán)
boxplot(exprSet["1939_at",]~pd$Disease)
boxplot(exprSet["1974_s_at",]~pd$Disease)
boxplot(exprSet["31618_at",]~pd$Disease)
Figure3.png

4.找到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

rm(list = ls())
options(stringsAsFactors = F)
#fill=T，不管如何都會讀取文件
a=read.table('e4-plot.txt',sep = '\t',fill = T,header = T)

colnames(a)=c('id','subtype','expression','mut')
BiocManager::install("ggstatsplot")
#畫boxplot，帶統(tǒng)計值的
library(ggstatsplot)
ggbetweenstats(data =a, x = subtype,  y = expression)
#b與a外觀一致，不知a其他兩列的作用是什么
b <- data.frame(subtype=a$subtype,expression=a$expression)
ggbetweenstats(data =b, x = subtype,  y = expression)
#保存
library(ggplot2)
ggsave('plot-again-BRCA1-TCGA-BRCA-cbioportal.png')
Figure4.png

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

提示使用：http://www.oncolnc.org/.整體不顯著的時候要看其亞型是否顯著。

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)
#將status列dead值變?yōu)?.alive值為0
dat$Status=ifelse(dat$Status=='Dead',1,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')
#查看亞型之間是否顯著
b=read.table('e4-plot.txt',sep = '\t',fill = T,header = T)
colnames(b)=c('Patient','subtype','expression','mut')
head(b)
b$Patient=substring(b$Patient,1,12)
tmp=merge(a,b,by='Patient')

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)  

#ifelse()的使用
ifelse(test, yes, no)
Figure5.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
提示：根據(jù)基因名拿到探針I(yè)D，縮小表達(dá)矩陣?yán)L制熱圖，沒有檢查到的基因直接忽略即可。

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

library(GEOquery)
# 這個包需要注意兩個配置，一般來說自動化的配置是足夠的。
#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)       ## 平臺文件
  save(gset,file=f)   ## 保存到本地
}
load('GSE17215_eSet.Rdata')  ## 載入數(shù)據(jù)
class(gset)
length(gset)
class(gset[[1]])
# 因為這個GEO數(shù)據(jù)集只有一個GPL平臺，所以下載到的是一個含有一個元素的list
a=gset[[1]]
dat=exprs(a)
dim(dat)


library(hgu133a.db)
ids=toTable(hgu133aSYMBOL)
head(ids)
dat=dat[ids$probe_id,]
dat[1:4,1:4] 
ids$median=apply(dat,1,median)
ids=ids[order(ids$symbol,ids$median,decreasing = T),]
ids=ids[!duplicated(ids$symbol),]
dat=dat[ids$probe_id,]
rownames(dat)=ids$symbol
dat[1:4,1:4]  
dim(dat)

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'
#拆分成合適的character后進(jìn)行提取
ng=strsplit(ng,' ')[[1]]
#%in%判斷在不在，match也可以
table(ng %in%  rownames(dat))
#取返回值為T的
ng=ng[ng %in%  rownames(dat)]
dat=dat[ng,]
#因為熱圖的值比較大，所以需要log一下
dat=log2(dat)
#scale=‘row’對row進(jìn)行歸一化
pheatmap::pheatmap(dat,scale = 'row')

7.下載數(shù)據(jù)集GSE24673的表達(dá)矩陣計算樣本的相關(guān)性并且繪制熱圖，需要標(biāo)記上樣本分組信息.【這里展示了一個操作，ctrl+F，對選定植進(jìn)行一鍵替換】

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

library(GEOquery)
# 這個包需要注意兩個配置，一般來說自動化的配置是足夠的。
#Setting options('download.file.method.GEOquery'='auto')
#Setting options('GEOquery.inmemory.gpl'=FALSE)
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]])
# 因為這個GEO數(shù)據(jù)集只有一個GPL平臺，所以下載到的是一個含有一個元素的list
a=gset[[1]]
dat=exprs(a)
dim(dat)
#查看分組信息
pd=pData(a)
#因為這里的分組信息亂，干脆就自己創(chuàng)建list
group_list=c('rbc','rbc','rbc',
             'rbn','rbn','rbn',
             'rbc','rbc','rbc',
             'normal','normal')
#做相關(guān)性之前要查看部分?jǐn)?shù)據(jù)的相關(guān)性
dat[1:4,1:4]
M=cor(dat)
pheatmap::pheatmap(M)
tmp=data.frame(g=group_list)
rownames(tmp)=colnames(M)
pheatmap::pheatmap(M,annotation_col = tmp)

8.找到 GPL6244 platform of Affymetrix Human Gene 1.0 ST Array 對應(yīng)的R的bioconductor注釋包，并且安裝它！
注意：所有的注釋R包后面都要加.db

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

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

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

library(GEOquery)
# 這個包需要注意兩個配置，一般來說自動化的配置是足夠的。
#Setting options('download.file.method.GEOquery'='auto')
#Setting options('GEOquery.inmemory.gpl'=FALSE)
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]])
# 因為這個GEO數(shù)據(jù)集只有一個GPL平臺，所以下載到的是一個含有一個元素的list
a=gset[[1]]
dat=exprs(a)
dim(dat)
pd=pData(a)
# (平均表達(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)
# 注意查看下載文件的大小，檢查數(shù)據(jù) 
f='GSE42872_eSet.Rdata'

library(GEOquery)
# 這個包需要注意兩個配置，一般來說自動化的配置是足夠的。
#Setting options('download.file.method.GEOquery'='auto')
#Setting options('GEOquery.inmemory.gpl'=FALSE)
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]])
# 因為這個GEO數(shù)據(jù)集只有一個GPL平臺，所以下載到的是一個含有一個元素的list
a=gset[[1]]
dat=exprs(a)
dim(dat)
pd=pData(a)
# (平均表達(dá)量/sd/mad/)最大的探針
boxplot(dat)
group_list=unlist(lapply(pd$title,function(x){
  strsplit(x,' ')[[1]][4]
}))


exprSet=dat
exprSet[1:4,1:4]
# DEG by limma limma接受的是log后的矩陣
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<-makeContrasts("progres.-stable",levels = design)
contrast.matrix 
##這個矩陣聲明，我們要把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)

啟示
1.我們不生產(chǎn)代碼，我們只是代碼的板運工。
2.理解R包原理后會修改一些較復(fù)雜圖的參數(shù)。
3.apply,%in%,match(),字符串操作(base和stringr，拆分用的較多)，sort，ifelse。
4.普通boxplot，帶統(tǒng)計值的boxplot，生存分析。
5.與網(wǎng)頁小工具的交替使用，可提高效率。

原文鏈接：http://www.bio-info-trainee.com/3750.html
B站視頻解答：https://www.bilibili.com/video/av25643438?p=13

Github原始代碼：https://github.com/jmzeng1314/R_bilibili

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容