作業(yè) 1
請(qǐng)根據(jù)R包org.Hs.eg.db找到下面ensembl 基因ID 對(duì)應(yīng)的基因名(symbol)
ENSG00000000003
ENSG00000000005
ENSG00000000419
ENSG00000000457
ENSG00000000460
ENSG00000000938
提示：
library(org.Hs.eg.db)
g2s=toTable(org.Hs.egSYMBOL)
g2e=toTable(org.Hs.egENSEMBL)
g2se <- merge(g2s,g2e,by.x="gene_id",by.y="gene_id")
index <- c("ENSG00000000003","ENSG00000000005","ENSG00000000419",
"ENSG00000000457","ENSG00000000460","ENSG00000000938")
g2se[1:5,1:3]
# gene_id symbol ensembl_id
#1 1 A1BG ENSG00000121410
#2 10 NAT2 ENSG00000156006
#3 100 ADA ENSG00000196839
#4 1000 CDH2 ENSG00000170558
#5 10000 AKT3 ENSG00000117020
g2se[g2se$ensembl_id %in% index,]
gene_id symbol ensembl_id
# gene_id symbol ensembl_id
# 11530 2268 FGR ENSG00000000938
# 21666 55732 C1orf112 ENSG00000000460
# 22360 57147 SCYL3 ENSG00000000457
# 23819 64102 TNMD ENSG00000000005
# 25603 7105 TSPAN6 ENSG00000000003
# 29268 8813 DPM1 ENSG00000000419
作業(yè) 2
作業(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
提示：
library(hgu133a.db)
ids=toTable(hgu133aSYMBOL)
head(ids)
# probe_id symbol
# 1 1053_at RFC2
# 2 117_at HSPA6
# 3 121_at PAX8
# 4 1255_g_at GUCA1A
# 5 1316_at THRA
# 6 1320_at PTPN21
index2 <- 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")
ids[ids$probe_id %in% index2,]
# probe_id symbol
# 1 1053_at RFC2
# 2 117_at HSPA6
# 3 121_at PAX8
# 4 1255_g_at GUCA1A
# 5 1316_at THRA
# 6 1320_at PTPN21
# 7 1405_i_at CCL5
# 8 1431_at CYP2E1
# 9 1438_at EPHB3
# 10 1487_at ESRRA
# 11 1494_f_at CYP2A6
# 12 1598_g_at GAS6
# 13 160020_at MMP14
# 14 1729_at TRADD
# 15 177_at PLD1
作業(yè) 3
找到R包CLL內(nèi)置的數(shù)據(jù)集的表達(dá)矩陣?yán)锩娴腡P53基因的表達(dá)量阻桅，并且繪制在 progres.-stable分組的boxplot圖
提示：
suppressPackageStartupMessages(library(CLL))
data(sCLLex)
sCLLex
exprSet=exprs(sCLLex)
group_list <- pData(sCLLex)
library(hgu95av2.db)
想想如何通過 ggpubr 進(jìn)行美化嫂沉。
ids3 <- toTable(hgu95av2SYMBOL)
ids3[ids3$symbol %in% "TP53",]
# probe_id symbol
# 966 1939_at TP53
# 997 1974_s_at TP53
# 1420 31618_at TP53
TP53 <- c(ids3[ids3$symbol %in% "TP53",][,1])
exprSet_TP53 <- exprSet[rownames(exprSet) %in% TP53,]
library(reshape2)
exprSet_TP53_L <- melt(exprSet_TP53)
colnames(exprSet_TP53_L) <- c("probe","sample","value")
exprSet_TP53_L$group <- rep(group_list$Disease,each=nrow(exprSet_TP53))
library(ggplot2)
ggplot(data = exprSet_TP53_L)+geom_boxplot(mapping = aes(x=group,
y=value))
作業(yè) 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
自己上網(wǎng)操作即可
作業(yè) 5
找到TP53基因在TCGA數(shù)據(jù)庫的乳腺癌數(shù)據(jù)集的表達(dá)量分組看其是否影響生存 提示使用：http://www.oncolnc.org/
自己上網(wǎng)操作即可
作業(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制熱圖，沒有檢查到的基因直接忽略即可搔啊。
library(GEOquery)
gset <- getGEO("GSE17215",destdir = ".",AnnotGPL = F,
getGPL = F)
exprSet6 <- exprs(gset[[1]])
Pdata6 <- pData(gset[[1]])
library(hthgu133a.db)
ids6 <- toTable(hthgu133aSYMBOL)
index6 <- c("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")
probe_idex <- ids6[ids6$symbol %in% index6,][,1]
exprSet_filter <- exprSet6[rownames(exprSet6) %in% probe_idex,]
exprSet_filter= t(scale(t(exprSet_filter)))
pheatmap::pheatmap(exprSet_filter)
image.png
作業(yè)7
下載數(shù)據(jù)集GSE24673的表達(dá)矩陣計(jì)算樣本的相關(guān)性并且繪制熱圖漫蛔，需要標(biāo)記上樣本分組信息
gset7 <- getGEO("GSE24673",destdir = ".",AnnotGPL = F,
getGPL = F)
##沒有分組情況下表達(dá)相關(guān)性初探
exprSet7 <- exprs(gset7[[1]])
Pdata7 <- pData(gset7[[1]])
cor_sample <- cor(exprSet7)
pheatmap::pheatmap(cor_sample)
#分組情況下表達(dá)相關(guān)性
group_list7 <- as.character(Pdata7[,8])
library(stringr)
group_list7=str_split(group_list7," -",simplify = T)[,2]
group_list7[10:11]=c("healthy","healthy")
cor_sample <- cor(exprSet7)
annotation_col = data.frame(
group = group_list7
)
rownames(annotation_col) =colnames(exprSet7)
pheatmap::pheatmap(cor_sample,annotation_col = annotation_col)
image.png
作業(yè)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("hugene10sttranscriptcluster.db",ask = F,update = F)
options()$repos
options()$BioC_mirror
#我本身有這個(gè)包剃毒，所有就直接加載了
library(hugene10sttranscriptcluster.db)
作業(yè)9
下載數(shù)據(jù)集GSE42872的表達(dá)矩陣赘阀，并且分別挑選出 所有樣本的(平均表達(dá)量/sd/mad/)最大的探針脑奠，并且找到它們對(duì)應(yīng)的基因。
library(GEOquery)
gset <- getGEO("GSE42872",destdir = ".",AnnotGPL = F,
getGPL = F)
###注意變量名字不要和前面幾道題的重復(fù)了
exprSet9 <- exprs(gset[[1]])
Pdata9 <- pData(gset[[1]])
library(hugene10sttranscriptcluster.db)
ids9 <- toTable(hugene10sttranscriptclusterSYMBOL)
tail(sort(table(ids9$symbol)))
#可以看到有些symbol對(duì)應(yīng)了多個(gè)探針
# RPL41 UBTFL1 CDK11B UBE2D3 IGKC LRRFIP1
# 6 6 8 8 10 10
table(sort(table(ids9$symbol)))
#18072個(gè)symol和探針是一一對(duì)應(yīng)的
# 1 2 3 4 5 6 8 10
# 18072 599 132 16 5 6 2 2
ids9=ids9[ids$symbol != "",]
#順序一致了
# 1 7896759 LINC01128
# 2 7896761 SAMD11
# 3 7896779 KLHL17
# 4 7896798 PLEKHN1
# 5 7896817 ISG15
ids9=ids9[ids9$probe_id %in% rownames(exprSet9),]
dat9=exprSet9
dim(dat9)
#過濾前有33279個(gè)探針
# 33297 6
dat9=dat9[ids9$probe_id,]
dim(dat9)
#過濾后19827個(gè)探針
# 19827 6
dat9[1:5,1:5]
# GSM1052615 GSM1052616 GSM1052617 GSM1052618 GSM1052619
# 7896759 8.75126 8.61650 8.81149 8.32067 8.41445
# 7896761 8.39069 8.52617 8.43338 9.17284 9.10216
# 7896779 8.20228 8.30886 8.18518 8.13322 8.06453
# 7896798 8.41004 8.37679 8.27521 8.34524 8.35557
# 7896817 7.72204 7.74572 7.78022 7.72308 7.53797
ids9[1:5,1:2]
####按照平均值
ids9$mean <- apply(dat9,1,mean)
ids9$max <- apply(dat9,1,max)
ids9$sd <- apply(dat9,1,sd)
ids9_mean <- ids9
ids9_max <- ids9
ids_sd <- ids9
ids9_mean=ids9_mean[order(ids9$symbol,ids9$mean,decreasing = T),]
dim(ids9_mean)
#去重復(fù)之前
# 19827 5
ids9_mean=ids9_mean[!duplicated(ids9_mean$symbol),]
dat9=dat9[ids9_mean$probe_id,]
dim(dat9)
#去重復(fù)之后
# 18834 6
rownames(dat9) <- ids9_mean$symbol
###后面的方法和之前一樣，注意重新導(dǎo)入dat9酸休。按照最大值方法
dat9=exprSet9
dim(dat9)
dat9=dat9[ids9$probe_id,]
dim(dat9)
ids9_max=ids9_max[order(ids9_max$symbol,ids9_max$max,decreasing = T),]
dim(ids9_max)
ids9_max=ids9_max[!duplicated(ids9_max$symbol),]
dat9=dat9[ids9_max$probe_id,]
dim(dat9)
rownames(dat9) <- ids9_max$symbol
###按照方差
dat9=exprSet9
dim(dat9)
dat9=dat9[ids9$probe_id,]
dim(dat9)
ids_sd=ids_sd[order(ids_sd$symbol,ids_sd$sd,decreasing = T),]
dim(ids_sd)
ids_sd=ids_sd[!duplicated(ids_sd$symbol),]
dat9=dat9[ids_sd$probe_id,]
dim(dat9)
rownames(dat9) <- ids_sd$symbol
作業(yè)10
下載數(shù)據(jù)集GSE42872的表達(dá)矩陣雨席，并且根據(jù)分組使用limma做差異分析吠式，得到差異結(jié)果矩陣
This entry was posted in 未分類 by ulwvfje. Bookmark the permalink.
#準(zhǔn)備好三個(gè)文件：過濾去重后的表達(dá)矩陣（上一步的dat9）、design分組文件糙置、contrast.matrix比較文件
exprSet9=dat9
Pdata9 <- pData(gset[[1]])
group_list9 <- str_split(Pdata9$title," ",simplify = T)[,4]
library(limma)
design <- model.matrix(~0+factor(group_list9))
colnames(design)=levels(factor(group_list9))
head(design)
# Control Vemurafenib
# 1 1 0
# 2 1 0
# 3 1 0
# 4 0 1
# 5 0 1
# 6 0 1
exprSet=dat9
rownames(design)=colnames(exprSet)
design
# Control Vemurafenib
# GSM1052615 1 0
# GSM1052616 1 0
# GSM1052617 1 0
# GSM1052618 0 1
# GSM1052619 0 1
# GSM1052620 0 1
# attr(,"assign")
# [1] 1 1
# attr(,"contrasts")
# attr(,"contrasts")$`factor(group_list9)`
# [1] "contr.treatment"
contrast.matrix<-makeContrasts("Vemurafenib-Control",
levels = design)
contrast.matrix ##這個(gè)矩陣聲明谤饭，我們要把 Tumor 組跟 Normal 進(jìn)行差異分析比較
# Contrasts
# Levels Vemurafenib-Control
# Control -1
# Vemurafenib 1
##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)
# logFC AveExpr t P.Value adj.P.Val B
# CD36 5.780170 7.370282 79.72556 1.209610e-16 2.278179e-12 26.75898
# DUSP6 -4.212683 9.106625 -62.43810 1.804535e-15 1.323536e-11 25.01000
# DCT 5.633027 8.763220 61.56547 2.108212e-15 1.323536e-11 24.89904
# SPRY2 -3.801663 9.726468 -53.95479 9.056119e-15 4.264073e-11 23.80849
# MOXD1 3.263063 10.171635 47.08154 4.074111e-14 1.305678e-10 22.59432
# ETV4 -3.843247 9.667077 -46.99304 4.159535e-14 1.305678e-10 22.57698
exprSet["CD36",]
#驗(yàn)證嗤疯，說明結(jié)果是正確的
# GSM1052615 GSM1052616 GSM1052617 GSM1052618 GSM1052619 GSM1052620
# 4.54610 4.40210 4.49239 10.25060 10.21480 10.31570
