中级题目 - 简书

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

rm(list = ls())
options(stringsAsFactors = F)
ENSG <- read.table('ensembl.txt')
library(org.Hs.eg.db)
ls('package:org.Hs.eg.db')
SYMBOL <- toTable(org.Hs.egSYMBOL )
ENSEMBL <- toTable(org.Hs.egENSEMBL)
library(tidyr)
library(dplyr)
ENSG <- separate(ENSG,V1,into = 'ensembl_id',sep = '[.]')
a <- merge(ENSG,ENSEMBL,by='ensembl_id')
b <- merge(a,SYMBOL,by='gene_id')##用merge而不用match?因为是两个data之间?

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

find_ids <- read.table('ids.txt')
library(hgu133a.db)
ls('package:hgu133a.db')
ids <- toTable(hgu133aSYMBOL)
colnames(find_ids)='probe_id'
c <- merge(find_ids,ids,by='probe_id')

3.找到R包CLL内置的数据集的表达矩阵里面的TP53基因的表达量，并且绘制在 progres.-stable分组的boxplot图..如何修改图,后面再系统学吧

library(CLL)
library(affy)
data(package='CLL')
data(sCLLex)
exprset <- exprs(sCLLex)
library(hgu95av2.db)   ??两个包有什么不一样?和前面那个?
ids=toTable(hgu95av2SYMBOL)
d <- ids['ids$symbol=='TP53'',]
d
e1 <- exprset['1939_at',]
e2 <-exprset['1974_s_at',]
e3 <- exprset['31618_at',]
e <- t(data.frame(e1,e2,e3))
pd <- pData(sCLLex)
boxplot(e1~pd$Disease)
boxplot(e2~pd$Disease)
boxplot(e3~pd$Disease)

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

提示：使用http://www.cbioportal.org/index.do 定位数据集：http://www.cbioportal.org/datasets###

options(stringsAsFactors = F)
BRCA1 <- read.table('BRCA1.txt',sep = '\t',fill = T,header = T)
head(BRCA1)
colnames(BRCA1) <- c('Sample.Id','Subtype','Expression','Mutations')
install.packages('ggstatsplot')
library(ggstatsplot)
ggbetweenstats(data =BRCA1, x = Subtype,  y = Expression)

image.png

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

提示使用：http://www.oncolnc.org/

image.png

网站可以直接做出来了?

BRCA_7157<-read.csv('BRCA_7157_50_50.csv')
library(ggplot2)
library(survival)
install.packages("survminer")
library(survminer) 
table(BRCA_7157$Status)
BRCA_7157$Status=ifelse(BRCA_7157$Status=='Dead',1,0)
sfit <- survfit(Surv(Days, Status)~Group, data=BRCA_7157)
ggsurvplot(sfit, conf.int=F, pval=TRUE)

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

options(stringsAsFactors = F)
GSE = "GSE17215_eSet.Rdata"
head(GSE,10)##???
library(GEOquery)
if(!file.exists(GSE)){
  gset <- getGEO('GSE17215', destdir=".",
                 AnnotGPL = F,     ## 注释文件
                 getGPL = F)       ## 平台文件
  save(gset,file=f)   ## 保存到本地
}
##得到对应GEO号的表达矩阵，注释信息，样本信息等

load('GSE17215_eSet.Rdata')  ## 载入数据
class(gset)
length(gset)
class(gset[[1]])
# 因为这个GEO数据集只有一个GPL平台，所以下载到的是一个含有一个元素的list
a=gset[[1]]
dat=exprs(a)
dim(dat)

head(dat)
library(hgu133a.db)
ids=toTable(hgu133aSYMBOL)
head(ids)
dat=dat[ids$probe_id,]
head(dat)
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'
ng=strsplit(ng,' ')[[1]]
table(ng %in%  rownames(dat))
ng=ng[ng %in%  rownames(dat)]
dat=dat[ng,]
dat=log2(dat)
pheatmap::pheatmap(dat,scale = 'row')

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

a=gset[[1]]
dat=exprs(a)
dim(dat)
pd=pData(a)
group_list=c('rbc','rbc','rbc',
             'rbn','rbn','rbn',
             'rbc','rbc','rbc',
             'normal','normal')
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 对应的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

9

rm(list=ls())
options(stringsAsFactors = F)

library(AnnoProbe)
library(Biobase)
gset <- geoChina('GSE42872')
eSet <- gset[[1]]
exprSet <-  exprs(eSet)

library(hugene10sttranscriptcluster.db)?????

IDs <- toTable(hugene10sttranscriptclusterSYMBOL)

IDs_means_max <- sort(apply(exprSet,1,mean),decreasing = T)[1]
names(IDs_means_max)
#[1] "7978905"


IDs_sd_max <- sort(apply(exprSet,1,sd),decreasing = T)[1]
names(IDs_sd_max)
#[1] "8133876"

IDs_mad_max <- sort(apply(exprSet,1,mad),decreasing = T)[1]
names(IDs_mad_max)
#[1] "8133876"

table(IDs$probe_id == names(IDs_means_max))
#FALSE 
#19814 
table(IDs$probe_id == names(IDs_sd_max))
#FALSE  TRUE 
#19813     1 
IDs[IDs$probe_id == names(IDs_sd_max),]
#      probe_id symbol
#16463  8133876   CD36
table(IDs$probe_id == IDs_mad_max)
#FALSE  TRUE 
#19813     1 
#IDs[IDs$probe_id == names(IDs_mad_max),]
#      probe_id symbol
#16463  8133876   CD36
IDs_means <- sort(apply(exprSet,1,mean),decreasing = T)
i <- 1
while (!(names(IDs_means)[i] %in% IDs$probe_id)) { i <- i+1}
i
#[1] 6
names(IDs_means)[i]
#[1] "7953385"
IDs[IDs$probe_id==names(IDs_means)[i],]
#      probe_id symbol
#4004  7953385  GAPDH

10

###10


#####

rm(list=ls())
options(stringsAsFactors = F)

library(AnnoProbe)
library(Biobase)
gset <- geoChina('GSE42872')
eset <- gset[[1]]
exprset <- exprs(eset)
pd <- pData(eset)

#goujianfenzuxinxi
group <- factor(rep(c('control','treat'),each = 3),levels = c('control','treat'))
group
design <- model.matrix(~0+group)
design


colnames(design)=c('control','treat')
rownames(design) <-  colnames(group)

##biaozhunhua
norm <- voom(exprset,design,plot = T)

##nihe 
fit <- lmFit(norm,design,method = 'ls')
contrast <- makeContrasts('treat-control',levels = design)

fit2 <- contrasts.fit(fit,contrast)
fit2 <- eBayes(fit2)
qqt(fit2$t, df = fit2$df.prior+fit2$df.residual, pch = 16, cex = 0.2)
qqt
abline(0,1)

diffgene <- topTable(fit2,number = Inf,adjust.method = 'fdr')
write.table(diffgene, 'gene_diff.txt', col.names = NA, sep = '\t', quote = FALSE)

diffgene[which(diffgene$logFC >= 1 & diffgene$adj.P.Val < 0.01),'sig'] <- 'red'
diffgene[which(diffgene$logFC <= -1 & diffgene$adj.P.Val < 0.01),'sig'] <- 'blue'
diffgene[which(abs(diffgene$logFC) < 1 | diffgene$adj.P.Val >= 0.01),'sig'] <- 'gray'
log2FoldChange <- diffgene$logFC
FDR <- diffgene$adj.P.Val


plot(log2FoldChange, -log10(FDR), pch = 20, col = diffgene$sig)
abline(v = 1, lty = 2)
abline(v = -1, lty = 2)
abline(h = -log(0.01, 10), lty = 2)