http://www.bio-info-trainee.com/3750.html

1、根据R包org.Hs.eg.db找到下面ensembl 基因ID 对应的基因名(symbol)

library(org.Hs.eg.db)
library(stringr)
#根据R包org.Hs.eg.db找到下面ensembl 基因ID 对应的基因名(symbol)
g2e = toTable(org.Hs.egENSEMBL)
g2s = toTable(org.Hs.egSYMBOL)
t = as.data.frame(c("ENSG00000000003.13",
      "ENSG00000000005.5",
      "ENSG00000000419.11",
     "ENSG00000000457.12",
      "ENSG00000000460.15",
      "ENSG00000000938.11"))
colnames(t) <- "ensemble"
t =str_sub(t$ensemble,1,15)
##这里有坑，要把后面的.xx这个版本号去掉
##调整两个数据框数据然后合并
identical(g2e$gene_id,g2s$gene_id)#看看顺序是不是一样，不一样就调整
table(g2s$gene_id %in% g2e$gene_id)
t_g2s <- g2s[g2s$gene_id %in% g2e$gene_id,]
t_g2e <- g2e[g2e$gene_id %in% t_g2s$gene_id,]
library(dplyr)
t_all = left_join(t_g2e,t_g2s,by = "gene_id")
result <- t_all$symbol[t_all$ensembl_id %in% t]

2、根据R包hgu133a.db找到下面探针对应的基因名(symbol)

#根据R包hgu133a.db找到下面探针对应的基因名(symbol)
library(hgu133a.db)
t = 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")
ls("package:hgu133a.db")#看看包里面有什么数据
symbol = toTable(hgu133aSYMBOL)
result = symbol$symbol[symbol$probe_id %in% t]

image.png

3、找到R包CLL内置的数据集的表达矩阵里面的TP53基因的表达量，并且绘制在 progres.-stable分组的boxplot图

#找到R包CLL内置的数据集的表达矩阵里面的TP53基因的表达量，并且绘制在 progres.-stable分组的boxplot图
BiocManager::install("CLL")
suppressPackageStartupMessages(library(CLL))#这个主要是是可以避免加载过多无用包

print(data(package = "CLL"))#看R包里面的数据集
sCLLex
pd <- pData(sCLLex)#不知道这个函数的可以？pData一下
exprSet=exprs(sCLLex) 
head(exprSet)
BiocManager::install("hgu95av2.db")
library("hgu95av2.db")#看看包里面的信息，注释使用这个包，我们要加载找出TP53对应的探针
ENSEMBL = toTable(hgu95av2ENSEMBL)
SYMbol = toTable(hgu95av2SYMBOL)#搜索可以发现有3个探针
?boxplot
boxplot(exprSet["31618_at",]~pd$Disease,data = exprSet,col = c("yellow","red"))
#用ggplot2美化  
library(ggplot2)
library(ggpubr)
exp <- exprSet["31618_at",]
group<- pd$Disease
result <- data.frame(exp,group)
p <-ggplot(result,aes(group,exp))
p + geom_boxplot(stat = "boxplot",aes(colour = group)) + xlab("Group") + ylab("Expression") +stat_compare_means(method = "t.test")

4、找到BRCA1基因在TCGA数据库的乳腺癌数据集(Breast Invasive Carcinoma (TCGA, PanCancer Atlas))的表达情况

http://www.cbioportal.org/datasets

5、找到TP53基因在TCGA数据库的乳腺癌数据集的表达量分组看其是否影响生存

http://www.oncolnc.org/
http://www.oncolnc.org/kaplan/?lower=25&upper=25&cancer=BRCA&gene_id=672&raw=BRCA1&species=mRNA
下载数据，放进R，作图

BRCA <- read.csv("BRCA_672_25_25.csv")
library(stringr)
library(survminer)
library(survival)
BRCA$Status <- str_replace_all(BRCA$Status,"Alive","0")
BRCA$Status <- str_replace_all(BRCA$Status,"Dead","1")
BRCA$Status <- as.numeric(BRCA$Status)
res.cox <- coxph(Surv(Days,Status) ~ Group, data = BRCA)
res.cox
fit <- survfit(Surv(Days,Status) ~ Group, data = BRCA)
ggsurvplot(fit,pval = TRUE)

image.png

6、下载数据集GSE17215的表达矩阵并且提取下面的基因画热图

library(GEOquery)
gse17215 = getGEO("GSE17215", destdir=".", AnnotGPL = F, getGPL = F) 
class(gse17215)
class(gse17215[[1]])
exp <- exprs(gse17215[[1]])
head(exp)
exp = log2(exp+1)#因为数据没有标准化
exp <- exprs(gse17215[[1]])
head(exp)
#http://www.bio-info-trainee.com/1399.html找到平台的注释包
gpl <- eSet[[1]]@annotation#看平台信息
BiocManager::install("hgu133a.db")
library(hgu133a.db)
ls("package:hgu133a.db")
symbol <- toTable(hgu133aSYMBOL)
genename <- "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"
gn <- strsplit(genename," ")[[1]]
#把探针和基因名字对应
symbol1 <- symbol[symbol$symbol %in% gn,]
exp1 <- data.frame(exp[rownames(exp) %in% symbol1$probe_id,])
symbol2 <- symbol1[match(rownames(exp1),symbol1$probe_id),]
exp1$symbol <- symbol2$symbol
exp1 <-exp1[!duplicated(exp1$symbol),]#行名不能重复
rownames(exp1) <- exp1$symbol  
exp1 <- exp1[,-ncol(exp1)]
pheatmap::pheatmap(exp1)

image.png

7、下载数据集GSE24673的表达矩阵计算样本的相关性并且绘制热图，需要标记上样本分组信息

library(GEOquery)
GSE24673 <- GEOquery::getGEO("GSE24673", destdir=".", AnnotGPL = F, getGPL = F) 
data <- data.frame(exprs(GSE24673[[1]]))
pd <- pData(GSE24673[[1]])
data2 <- rbind(pd$`cell type:ch1`,data)
library(pheatmap)
annotation_col <- data.frame(factor(pd$source_name_ch1))
colnames(annotation_col) <- "type" 
rownames(annotation_col) <- colnames(data)
pheatmap(data,annotation_col = annotation_col,cluster_cols = F,annotation_names_col = F,show_rownames = F)

data.cor <- cor(data)#计算相关性
pheatmap(data.cor,annotation_col = annotation_col,cluster_cols = T,annotation_names_col = F,show_rownames = F)

image.png

8、找到 GPL6244 platform of Affymetrix Human Gene 1.0 ST Array 对应的R的bioconductor注释包，并且安装它！

http://www.bio-info-trainee.com/1399.html

hugene10sttranscriptcluster
BiocManager::install( "hugene10sttranscriptcluster")

9、下载数据集GSE42872的表达矩阵，并且分别挑选出 所有样本的(平均表达量/sd/mad/)最大的探针，并且找到它们对应的基因。

library(GEOquery)
GSE42872 <- getGEO("GSE42872", AnnotGPL = F,getGPL = F)
data <- as.data.frame(exprs(GSE42872[[1]]))
pd <- pData(GSE42872[[1]])
head(data)

#排序
me <- names(sort(apply(data,1,mean),decreasing=T)[1])
sd <- names(sort(apply(data,1,sd),decreasing=T)[1])
mad <- names(sort(apply(data,1,mad),decreasing=T)[1])

#找探针对应的基因
BiocManager::install( "hugene10sttranscriptcluster.db")
library("hugene10sttranscriptcluster.db")
ls("package:hugene10sttranscriptcluster.db")
SYMBOY <- toTable(hugene10sttranscriptclusterSYMBOL)
select(hugene10sttranscriptcluster.db, keys=c(me,sd,mad), columns = c("SYMBOL"))
SYMBOY[SYMBOY$probe_id %in% c(me,sd,mad),]
#发现"7978905"没有对应的基因

image.png

10、下载数据集GSE42872的表达矩阵，并且根据分组使用limma做差异分析，得到差异结果矩阵

#下载数据集GSE42872的表达矩阵，并且根据分组使用limma做差异分析，得到差异结果矩阵
library(limma)
library(GEOquery)
library(stringr)
GSE42872 <- getGEO("GSE42872", AnnotGPL = F,getGPL = F)
data <- as.data.frame(exprs(GSE42872[[1]]))
pd <- pData(GSE42872[[1]])
head(data)
pdgroup <- pd$title
grouplist <- unlist(lapply(pdgroup, function(x){strsplit(x," ")[[1]][4]}))
grouplist <- factor(grouplist,levels = c("Control","Vemurafenib"))

#差异
library(limma)
design=model.matrix(~grouplist)
fit=lmFit(data,design)
fit=eBayes(fit)
deg=topTable(fit,coef=2,number = Inf)
#整理数据
library(dplyr)
library("hugene10sttranscriptcluster.db")
SYMBOY <- toTable(hugene10sttranscriptclusterSYMBOL)
deg <- mutate(deg,probe_id = rownames(deg))
deg <- inner_join(deg,SYMBOY,by = "probe_id")
deg <- deg[!duplicated(deg$symbol),]
#标记上下调基因
logFC_t=1
P.Value_t = 0.01
k1 = (deg$P.Value < P.Value_t)&(deg$logFC < -logFC_t)
k2 = (deg$P.Value < P.Value_t)&(deg$logFC > logFC_t)
change = ifelse(k1,"down",ifelse(k2,"up","stable"))

deg <- mutate(deg,change)