bilibili视频-《生信人应该这样学R语言》-P13-R语言
2019-06-06 本文已影响0人
小梦游仙境
R语言小作业-中级
1.请根据R包org.Hs.eg.db找到下面ensembl 基因ID 对应的基因名(symbol)
ENSG00000000003.13
ENSG00000000005.5
ENSG00000000419.11
ENSG00000000457.12
ENSG00000000460.15
ENSG00000000938.11
新建一个e1.txt的文档
rm(list=ls())
options(stringsAsFactors = F)#如果不加后面会在strsplit(x,'[.]')时提示：non-character argument.
a=read.table('e1.txt')
head(a)
library(org.Hs.eg.db)
ls("package:org.Hs.eg.db") #图1
g2s <- toTable(org.Hs.egSYMBOL) #图2
g2e <- toTable(org.Hs.egENSEMBL) #图2
library(stringr)#只要"."之前的ensembl ID
#或者用lapply函数
lapply(a$V1,function(x){
strsplit(as.character(x),'[.]')[[1]][1]#以"."为分割
})#在没加options(stringsAsFactors = F的情况下加as.character一样
unlist(lapply(a$V1,function(x){
+ strsplit((x),'[.]')[[1]][1]
+ })) #图3
a$ensembl_id=unlist(lapply(a$V1,function(x){
strsplit((x),'[.]')[[1]][1]
}))
tmp <- merge(a,g2e, by="ensembl_id")#图4
tmp <- merge(tmp,g2s, by="gene_id")#图5
图1
图2
图3
图4
图5
2. 根据R包hgu133a.db找到下面探针对应的基因名(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
rm(list=ls())
options(stringsAsFactors = F)
a=read.table('e2.txt')#此时a是新建的e2.txt
colnames(a)='probe_id' #图7 a的列名为V1，而与ids的列名不同，因此把a的列名改为与ids相同的'probe_id'，用于后面merge。
library(hgu133a.db)
ls("package:hgu133a.db")#图6
ids=toTable(hgu133aSYMBOL)#图7
head(ids)
tmp1=merge(ids,a,by='probe_id')#图7
tmp2=ids[match(a$'probe_id',ids$'probe_id'),]#图7 用match函数得到前一个向量在后一个向量的坐标,再在ids中取出来，','左边即按照这些坐标把这坐标对应的这行取出来
## 附：判断得到的两组结果是否一致
# 法一：
identical(tmp1,tmp2) #返回逻辑值
# 法二：
dplyr::setdiff(tmp1,tmp2) #返回两组的差别【没差就返回空】
图6
图7
3.找到R包CLL内置的数据集的表达矩阵里面的TP53基因的表达量，并且绘制在 progres.-stable分组的boxplot图
rm(list=ls())
options(stringsAsFactors = F)
suppressPackageStartupMessages(library(CLL))
data(sCLLex)
sCLLex
exprSet=exprs(sCLLex) #可通过exprs来获取这个对象的表达矩阵exprSet
# sCLLex
exprSet <- exprs(sCLLex) #探针的表达量
pd<- pData(sCLLex) #sampleID与disease的对应关系-分组信息
p2s <- toTable(hgu95av2SYMBOL) #探针与symbol的对应关系
p3 <- filter(p2s,symbol=='TP53')
library(hgu95av2.db)
ids=toTable(hgu95av2SYMBOL)#图8
head(ids)
boxplot(exprSet['1939_at',]~pd$Disease)#图9
boxplot(exprSet['1974_s_at',]~pd$Disease)#图9
boxplot(exprSet['31618_at',]~pd$Disease)#图9
##接下来这段摘自小洁（https://www.jianshu.com/p/1766eae6f97a），留下来
probe_tp53 <- c("1939_at","1974_s_at","31618_at")
i = 3 #可换1，2
boxplot(exprSet[probe_tp53[i],] ~ pdata$Disease)
#用ggpubr作图
#http://www.sthda.com/english/articles/24-ggpubr-publication-ready-plots/
exp_tab <- rownames_to_column(as.data.frame(exprSet))
exp_tab2 <- gather(exp_tab,
key = 'sample',
value = 'exp',-1)
pdata <- rownames_to_column(pdata)
exp_tab3 <- merge(exp_tab2,pdata,by.x='sample',by.y='rowname')
i=1 ###可换1，2
dev.off()
p <- ggboxplot(filter(exp_tab3,rowname==probe_tp53[i]),
x = 'Disease',
y = 'exp',
color = "Disease", palette =c("#00AFBB", "#E7B800", "#FC4E07"),
add = "jitter", shape = "Disease")
p
图8-发现TP53有3个probe_id
图9
4.找到BRCA1基因在TCGA数据库的乳腺癌数据集(Breast Invasive Carcinoma (TCGA, PanCancer Atlas))的表达情况
提示：使用http://www.cbioportal.org/index.do 定位数据集：http://www.cbioportal.org/datasets见图10
图10
#下载数据后自己做一遍
rm(list=ls())
options(stringsAsFactors = F)
#ID，四个亚型，表达量
a= read.table("e4-plot.txt", sep = "\t",fill=T,header=T)
## boxplot
colnames(a) <- c("id", "subtype", "expression", "mut")
dat=a
library(ggstatsplot)
ggbetweenstats(data = dat,
x = subtype,
y = expression)
library(ggplot2)
ggsave("e4-BRCA1-TCGA.png")#图11ggsave保存需要下载ggplot2这个包
图11
5.找到TP53基因在TCGA数据库的乳腺癌数据集的表达量分组看其是否影响生存
提示使用：http://www.oncolnc.org/见图12、图13
图12
图13
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)
dat$Status=ifelse(dat$Status=='Dead',1,0)#图14 将status改名，Dead改为1，Alive就改为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')#图15
图14
图15
head(a)
b=read.table('e4-plot.txt',sep = '\t',fill = T,header = T)
colnames(b)=c('Patient','subtype','expression','mut')#图17
head(b)#此处的b其实是上一题中的a，为了与本题中的a的Patient ID相一致，此时把b的第一列的名字命名为Patient.
b$Patient=substring(b$Patient,1,12)#图16记住substring取字符用法
tmp=merge(a,b,by='Patient')#图17
图16
图17
可以检查下merge后的tmp肯定会有减少，这时需要注意下看看还有多少个观测，如图显示还有800多，还不错。见图18
图18
#如果不写循环
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) #图19
图19
6.下载数据集GSE17215的表达矩阵并且提取下面的基因画热图
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
提示：根据基因名拿到探针ID，缩小表达矩阵绘制热图，没有检查到的基因直接忽略即可。
##下载表达矩阵
rm(list = ls()) ## 魔幻操作，一键清空~
options(stringsAsFactors = F)
# 注意查看下载文件的大小，检查数据
f='GSE17215_eSet.Rdata'
library(GEOquery)#加载GEO这个包
# 这个包需要注意两个配置，一般来说自动化的配置是足够的。
#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') ## 载入数据
class(gset)
length(gset)
class(gset[[1]])
# 因为这个GEO数据集只有一个GPL平台，所以下载到的是一个含有一个元素的list
a=gset[[1]]
dat=exprs(a)#图21-1-1（图21的第1个图）
dim(dat)
#再进一步过滤基因hgu133a.db这个包
library(hgu133a.db)#图20
ids=toTable(hgu133aSYMBOL)
head(ids)
dat=dat[ids$probe_id,]#图21-2（图21的第2个图）
dat[1:4,1:4]
ids$median=apply(dat,1,median)#图22-1.对dat每一行求median值
ids=ids[order(ids$symbol,ids$median,decreasing = T),]#图22-2.order是按从小到达排序，加了decreasing = T是按从大到小排序。ids$symbol是将symbol基因名按照字母顺序从最后即"z"排序的，然后每个z开头的再按ids$median排序
ids=ids[!duplicated(ids$symbol),]#图22-3.将ids$symbol去重后作为新的行，ids取这些行
dat=dat[ids$probe_id,]#图23-3和图23-2.dat这个表达矩阵取出新的ids的probi_id这整行为一个新dat
rownames(dat)=ids$symbol#图23-2
dat[1:4,1:4]
dim(dat)
图20
图21
图22
图23
# 想要的基因
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]]#图24
ng
table(ng %in% rownames(dat))#看我们想要的基因有哪些在dat（GSE17215）这个基因集里
ng=ng[ng %in% rownames(dat)]
dat=dat[ng,]
dat=log2(dat)
pheatmap::pheatmap(dat,scale = 'row')#图25.把row用scale归一化
图24
图25
7.下载数据集GSE24673的表达矩阵计算样本的相关性并且绘制热图，需要标记上样本分组信息
rm(list = ls()) ## 魔幻操作，一键清空~
options(stringsAsFactors = F)
# 注意查看下载文件的大小，检查数据
f='GSE24673_eSet.Rdata'
library(GEOquery)
if(!file.exists(f)){
gset <- getGEO('GSE24673', destdir=".",
AnnotGPL = F,
getGPL = F)
save(gset,file=f)
}
load('GSE24673_eSet.Rdata') ## 载入数据
class(gset)
length(gset)
class(gset[[1]])
# 因为这个GEO数据集只有一个GPL平台，所以下载到的是一个含有一个元素的list
a=gset[[1]]
dat=exprs(a)#图26.用exprs取出表达矩阵
dim(dat)
pd=pData(a)#图27.用pData取出临床信息，即可以得到分组信息
group_list=c('rbc','rbc','rbc', #看了pd后创建分组信息
'rbn','rbn','rbn',
'rbc','rbc','rbc',
'normal','normal')
dat[1:4,1:4]
M=cor(dat)#图28.看样本和样本之间的相关性
pheatmap::pheatmap(M)
tmp=data.frame(g=group_list)#图29
rownames(tmp)=colnames(M)
pheatmap::pheatmap(M,annotation_col = tmp)#图28 分组信息用annotation_col.
图26
图27
图28
图29
图30
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("请输入自己找到的R包",ask = F,update = F)
options()$repos
options()$BioC_mirror
9.下载数据集GSE42872的表达矩阵，并且分别挑选出所有样本的(平均表达量/sd/mad/)最大的探针，并且找到它们对应的基因。
####依然同上一题一样 下载矩阵
rm(list = ls())
options(stringsAsFactors = F)
f='GSE42872_eSet.Rdata'
library(GEOquery)
if(!file.exists(f)){
gset <- getGEO('GSE42872', destdir=".",
AnnotGPL = F,
getGPL = F)
save(gset,file=f)
}
load('GSE42872_eSet.Rdata') ## 载入数据
class(gset)
length(gset)
class(gset[[1]])
# 因为这个GEO数据集只有一个GPL平台，所以下载到的是一个含有一个元素的list
a=gset[[1]]
dat=exprs(a)#用exprs取出表达矩阵
dim(dat)
pd=pData(a)#用pData取出临床信息，即可以得到分组信息
# (平均表达量/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.下载数据集GSE42872的表达矩阵，并且根据分组使用limma做差异分析，得到差异结果矩阵
rm(list = ls())
options(stringsAsFactors = F)
f='GSE42872_eSet.Rdata'
library(GEOquery)
if(!file.exists(f)){
gset <- getGEO('GSE42872', destdir=".",
AnnotGPL = F,
getGPL = F)
save(gset,file=f)
}
load('GSE42872_eSet.Rdata') ## 载入数据
class(gset)
length(gset)
class(gset[[1]])
# 因为这个GEO数据集只有一个GPL平台，所以下载到的是一个含有一个元素的list
a=gset[[1]]
dat=exprs(a)
dim(dat)
pd=pData(a)#图31
# (平均表达量/sd/mad/)最大的探针
boxplot(dat)
group_list=unlist(lapply(pd$title,function(x)
strsplit(x,' ')[[1]][4]
))
exprSet=dat#目前的矩阵是dat，转换成limma包里代码用的exprSet一致
exprSet[1:4,1:4]
# DEG by limma
suppressMessages(library(limma))
design <- model.matrix(~0+factor(group_list))#图31
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进行差异分析比较
##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)
图31
