當(dāng)看完幾遍R語(yǔ)言書(shū)籍時(shí)躯嫉，可以適當(dāng)做些題目，而jimmy在生信技能樹(shù)的這幾道題目是很適合練習(xí)雳灵，每一行代碼自己都要著去理解，在Rstudio中可以不斷修改和試錯(cuò)闸盔，直到弄明白悯辙。http://www.bio-info-trainee.com/3750.html 其實(shí)這些題目有很多好的答案了。
但佛系筆記嘛迎吵，可能在做筆記前花了很多時(shí)間才搞懂為什么是這樣寫(xiě)那樣寫(xiě)躲撰，身邊有很多小伙伴可以一起學(xué)習(xí)，終于做出來(lái)的感覺(jué)挺爽的击费。
第一題
請(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
提示：
library(org.Hs.eg.db)
g2s=toTable(org.Hs.egSYMBOL)
g2e=toTable(org.Hs.egENSEMBL)
先甩個(gè)生信技能樹(shù)博客鏈接—了解常用數(shù)據(jù)庫(kù)ID表示方式 http://www.biotrainee.com/thread-411-1-1.html
答1
> rm(list=ls()) #刪除已存在變量拢蛋，可刪可不刪，看個(gè)人習(xí)慣
> if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("org.Hs.eg.db") # Bioconductor 命令行下載安裝org.Hs.eg.db 包
> suppressMessages(library(org.Hs.eg.db)) #去除加載包時(shí)的警告蔫巩，可用可不用
> ls("package:org.Hs.eg.db") #查看包的內(nèi)容
[1] "org.Hs.eg" "org.Hs.eg.db"
[3] "org.Hs.eg_dbconn" "org.Hs.eg_dbfile"
[5] "org.Hs.eg_dbInfo" "org.Hs.eg_dbschema"
[7] "org.Hs.egACCNUM" "org.Hs.egACCNUM2EG"
[9] "org.Hs.egALIAS2EG" "org.Hs.egCHR"
[11] "org.Hs.egCHRLENGTHS" "org.Hs.egCHRLOC"
[13] "org.Hs.egCHRLOCEND" "org.Hs.egENSEMBL"
[15] "org.Hs.egENSEMBL2EG" "org.Hs.egENSEMBLPROT"
[17] "org.Hs.egENSEMBLPROT2EG" "org.Hs.egENSEMBLTRANS"
[19] "org.Hs.egENSEMBLTRANS2EG" "org.Hs.egENZYME"
[21] "org.Hs.egENZYME2EG" "org.Hs.egGENENAME"
[23] "org.Hs.egGO" "org.Hs.egGO2ALLEGS"
[25] "org.Hs.egGO2EG" "org.Hs.egMAP"
[27] "org.Hs.egMAP2EG" "org.Hs.egMAPCOUNTS"
[29] "org.Hs.egOMIM" "org.Hs.egOMIM2EG"
[31] "org.Hs.egORGANISM" "org.Hs.egPATH"
[33] "org.Hs.egPATH2EG" "org.Hs.egPFAM"
[35] "org.Hs.egPMID" "org.Hs.egPMID2EG"
[37] "org.Hs.egPROSITE" "org.Hs.egREFSEQ"
[39] "org.Hs.egREFSEQ2EG" "org.Hs.egSYMBOL"
[41] "org.Hs.egSYMBOL2EG" "org.Hs.egUCSCKG"
[43] "org.Hs.egUNIGENE" "org.Hs.egUNIGENE2EG"
[45] "org.Hs.egUNIPROT"
> library("org.Hs.eg.db") #加載包
> keytypes(org.Hs.eg.db) # org.Hs.eg.db 包里面的數(shù)據(jù)有哪些類(lèi)型谆棱，可以看到這題中我們需要的是ENSEMBL，和SYMBOL
[1] "ACCNUM" "ALIAS" "ENSEMBL" "ENSEMBLPROT" "ENSEMBLTRANS"
[6] "ENTREZID" "ENZYME" "EVIDENCE" "EVIDENCEALL" "GENENAME"
[11] "GO" "GOALL" "IPI" "MAP" "OMIM"
[16] "ONTOLOGY" "ONTOLOGYALL" "PATH" "PFAM" "PMID"
[21] "PROSITE" "REFSEQ" "SYMBOL" "UCSCKG" "UNIGENE"
[26] "UNIPROT"
> x<-c("ENSG00000000003","ENSG00000000005","ENSG00000000419","ENSG00000000457","ENSG00000000460","ENSG00000000938")
# 將這些ENSEMBL ID 賦值給x
> cols <- c("SYMBOL", "GENENAME")
> z <- select(org.Hs.eg.db,keys=x,columns=cols,keytype="ENSEMBL")
# select 函數(shù)圆仔，將ensembl ID 對(duì)應(yīng)到 gene name,賦值給z
> View(z) #查看z
第二題
根據(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提示：
library(hgu133a.db)
ids=toTable(hgu133aSYMBOL)
head(ids)
答2
> if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("hgu95av2.db") #Bioconductor 下載安裝hgu95av2.db 包垃瞧，
> suppressMessages(library(hgu133a.db)) #去除加載包的警告
> library(hgu133a.db) #加載包
> keytypes(hgu133a.db) #查看包的的內(nèi)容
[1] "ACCNUM" "ALIAS" "ENSEMBL" "ENSEMBLPROT" "ENSEMBLTRANS"
[6] "ENTREZID" "ENZYME" "EVIDENCE" "EVIDENCEALL" "GENENAME"
[11] "GO" "GOALL" "IPI" "MAP" "OMIM"
[16] "ONTOLOGY" "ONTOLOGYALL" "PATH" "PFAM" "PMID"
[21] "PROBEID" "PROSITE" "REFSEQ" "SYMBOL" "UCSCKG"
[26] "UNIGENE" "UNIPROT"
> x <-c("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")
> cols <- c("PROBEID","SYMBOL")
> tmp<- select(hgu133a.db,keys=keys(hgu133a.db),columns = cols) #select 函數(shù)得到
> tmp2<-tmp[match(x,tmp$PROBEID),] #match 函數(shù)
> View(tmp2)
第三題:
找到R包CLL內(nèi)置的數(shù)據(jù)集的表達(dá)矩陣?yán)锩娴腡P53基因的表達(dá)量，并且繪制在 progres.-stable分組的boxplot圖提示：
suppressPackageStartupMessages(library(CLL))
data(sCLLex)
sCLLex
exprSet=exprs(sCLLex)
library(hgu95av2.db)
想想如何通過(guò) ggpubr 進(jìn)行美化坪郭。
答3
> if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("CLL") #下載安裝CLL包
> suppressMessages(library(CLL))
> library(CLL) #加載CLL包
> data(sCLLex) #調(diào)用sCLLex
> library(hgu95av2.db) #加載hgu95av2.db包
> exprSet <- exprs(sCLLex)
> pdata <- pData(sCLLex)
> p2s <- toTable(hgu95av2SYMBOL)
> p3 <- filter(p2s,symbol=='TP53')
第四題
找到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
可以在網(wǎng)站中直接查看个从，也可以用將數(shù)據(jù)下載后，R代碼得到
第五題
找到TP53基因在TCGA數(shù)據(jù)庫(kù)的乳腺癌數(shù)據(jù)集的表達(dá)量分組看其是否影響生存
找到官網(wǎng)歪沃，[http://www.oncolnc.org/](http://www.oncolnc.org/)
第六題
下載數(shù)據(jù)集GSE17215的表達(dá)矩陣并且提取下面的基因畫(huà)熱圖提示：根據(jù)基因名拿到探針I(yè)D嗦锐，縮小表達(dá)矩陣?yán)L制熱圖，沒(méi)有檢查到的基因直接忽略即可沪曙。
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
第七題
下載數(shù)據(jù)集GSE24673的表達(dá)矩陣計(jì)算樣本的相關(guān)性并且繪制熱圖奕污，需要標(biāo)記上樣本分組信息
rm(list=ls()) #刪除已有變量
options(stringsAsFactors = F)
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install('GEOquery')
library(GEOquery) #下載安裝GEOquery
GSE <- "GSE24673" 下載GSE24673 數(shù)據(jù)集
if(!file.exists(GSE)){
geo <- getGEO(GSE, destdir = '.', getGPL = F, AnnotGPL = F)
save(geo, file = paste0(GSE,'.eSet.Rdata'))
} #if 函數(shù) 得到并保存GSE,'.eSet.Rdata' 命名為pasteo文件，賦值給geo
load(paste0(GSE,'.eSet.Rdata')) #加載pasteo
expr <- exprs(geo[[1]]) #獲得geo的exprimentData數(shù)據(jù)珊蟀，并賦值給expr
dim(expr) 查看expr有幾行幾列
expr[1:4,1:4] #查看expr中一至四行菊值，一至四列的內(nèi)容
pdata <- pData(geo[[1]]) #提取phenoData的數(shù)據(jù)，并賦值給pdata
# 自己根據(jù)pdata第八列做一個(gè)分組信息矩陣
grp <- c('rbc','rbc','rbc',
'rbn','rbn','rbn',
'rbc','rbc','rbc',
'normal','normal')
grp_df <- data.frame(group=grp) #將grp 得到一個(gè)數(shù)據(jù)框
rownames(grp_df) <- colnames(expr) #列名等于expr中的名字
new_cor <- cor(expr)
pheatmap::pheatmap(new_cor, annotation_col = grp_df) #熱圖
第八題
找到 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/")#Bioconductor 鏡像
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")) #清華鏡像
BiocManager::install("hugene10sttranscriptcluster",ask = F,update =T) #R包是"hugene10sttranscriptcluster"，更新到最新版本
options()$repos
options()$BioC_mirror
#參考芯片與注釋包的關(guān)系 [http://www.reibang.com/p/f6906ba703a0](http://www.reibang.com/p/f6906ba703a0)
第九題
下載數(shù)據(jù)集GSE42872的表達(dá)矩陣磅崭，并且分別挑選出 所有樣本的(平均表達(dá)量/sd/mad/)最大的探針儿子，并且找到它們對(duì)應(yīng)的基因。
library(GEOquery) #加載GEOquery包
GSE<-GSE42872
if(!file.exists(GSE)){
geo <- getGEO(GSE, destdir = '.', getGPL = F, AnnotGPL = F)
save(geo, file = paste0(GSE,'.eSet.Rdata'))
}
exprSet = exprs(gset[[1]]) #獲得表達(dá)矩陣
# exprSet <- exprSet[apply(exprSet, 1, function(x) sum(x>1) > 5),]
# exprSet <- log(edgeR::cpm(exprSet) + 1)
# 處理中平均值最大的探針
max_mean <- names(sort(apply(exprSet,1, mean), decreasing = T))[1]
# 處理中sd最大的探針
max_sd <- names(sort(apply(exprSet,1, sd), decreasing = T))[1]
# 處理中mad最大的探針
max_mad <- names(sort(apply(exprSet,1, mad), decreasing = T))[1]
# 查詢(xún)探針I(yè)D對(duì)應(yīng)基因的symbol
library(hugene10sttranscriptcluster.db)
select(hugene10sttranscriptcluster.db,
keys = c(max_mean, max_sd, max_mad),
columns = c("SYMBOL"),
keytype = "PROBEID")
第十題
下載數(shù)據(jù)集GSE42872的表達(dá)矩陣砸喻，并且根據(jù)分組使用limma做差異分析柔逼，得到差異結(jié)果矩陣
gset = getGEO('GSE42872', GSEMatrix = TRUE, AnnotGPL = FALSE, getGPL= FALSE)
s <- toTable(hugene10sttranscriptclusterSYMBOL)
exprSet = exprs(gset[[1]])
group_list <- as.character(c(rep("Control",3),rep("Vemurafenib",3)))
dim(exprSet)
# 過(guò)濾【只保留和注釋文件探針id相同的探針】
efilt <- exprSet[rownames(exprSet)%in%s$probe_id,]
dim(efilt)
## 整合1【目的：保證一個(gè)基因?qū)?yīng)一個(gè)探針；如果基因和探針一一對(duì)應(yīng)很好說(shuō)割岛，但如果一個(gè)基因?qū)?yīng)多個(gè)探針：每個(gè)探針取一行的均值-》對(duì)應(yīng)同一基因的探針取表達(dá)量最大的探針-》按照基因名給他們建索引愉适，因?yàn)槭前凑栈騺?lái)過(guò)濾探針（不用s$probe_id構(gòu)建索引的原因是，看清楚我們的目的是讓注釋包的一個(gè)基因?qū)?yīng)我們自己表達(dá)矩陣的一個(gè)探針癣漆。如果用s$probe_id那么結(jié)果就成了讓注釋包的一個(gè)探針對(duì)應(yīng)我們自己表達(dá)矩陣的一個(gè)探針维咸，當(dāng)然這樣也運(yùn)行不成功，因?yàn)樽约罕磉_(dá)矩陣的探針過(guò)濾后的數(shù)量和注釋包的探針數(shù)量不相等惠爽，這樣沒(méi)法一一對(duì)應(yīng)癌蓖。但基因名數(shù)量是不變的，什么是索引婚肆？以不變應(yīng)萬(wàn)變的就是索引）】
maxp = by(efilt,s$symbol,function(x) rownames(x)[which.max(rowMeans(x))])
uniprobes = as.character(maxp)
efilt=efilt[rownames(efilt)%in%uniprobes,]
## 整合2【目的：將我們表達(dá)矩陣的行名換成剛才一對(duì)一的基因名租副，并且match這個(gè)函數(shù)保證了表達(dá)矩陣和注釋包的順序是一致的】
rownames(efilt)=s[match(rownames(efilt),s$probe_id),2]
# 差異分析
suppressMessages(library(limma))
#limma需要三個(gè)矩陣：表達(dá)矩陣（efilt）、分組矩陣（design）较性、比較矩陣（contrast）
#先做一個(gè)分組矩陣～design用僧，說(shuō)明MAO是哪幾個(gè)樣本，MNO又是哪幾個(gè)赞咙，其中1代表“是”
design <- model.matrix(~0+factor(group_list))
colnames(design) <- levels(factor(group_list))
rownames(design) <- colnames(efilt)
design
#再做一個(gè)比較矩陣【一般是case比control】
contrast<-makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)
contrast
DEG <- function(efilt,design,contrast){
##step1
fit <- lmFit(efilt,design)
##step2
fit2 <- contrasts.fit(fit, contrast)
fit2 <- eBayes(fit2)
##step3
mtx = topTable(fit2, coef=1, n=Inf)
deg_mtx = na.omit(mtx)
return(deg_mtx)
}
DEG_mtx <- DEG(efilt,design,contrast) #得到全部的差異基因矩陣
