【生信技能树】R语言中级练习题

2020-03-29 本文已影响0人函数IsLazy

R语言小作业-中级
代码来自Jimmy的Github 所以更多的精力是理解注释代码。

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

ENSG00000000003.13
ENSG00000000005.5
ENSG00000000419.11
ENSG00000000457.12
ENSG00000000460.15
ENSG00000000938.11

rm(list = ls())
options(stringsAsFactors = F)
a=read.table("e1.txt")  #把上面的ensemble ID放到e1.txt
library(org.Hs.eg.db)
ls("package:org.Hs.eg.db")
g2s = toTable(org.Hs.egSYMBOL);head(g2s)
g2e = toTable(org.Hs.egENSEMBL)
#去掉小数点之后的数字来提取ensemble_ID
a$ensembl_id=unlist(lapply(a$V1,function(x){
  strsplit(x,"[.]")[[1]][1]
})
)
#合并ID来提取基因名
tmp=merge(a,g2e,by='ensembl_id')
tmp=merge(tmp,g2s,by='gene_id')

根据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)
b = read.table("e2.txt")   # 上面放入e2.txt
colnames(b)='probe_id'  #改一下名字 方便合并
library(hgu133a.db)
ids = toTable(hgu133aSYMBOL)
head(ids)
# 根据ID名来合并文件
tmp = merge(ids,b,by='probe_id')
#也可以用match
match(b$probe_id,ids$probe_id)
tmp2=ids[match(b$probe_id,ids$probe_id),]

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

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

suppressPackageStartupMessages(library(CLL)) # 静默加载，不会出现下面的红字
data(sCLLex)
exprSet = exprs(sCLLex) # 提取表达矩阵
pd = pData(sCLLex)  # 检查分组信息
library(hgu95av2.db))
ids = toTable(hgu95av2SYMBOL) # 把探针信息提取出来，并且搜索TP53对应的探针名，有三个
head(ids)
boxplot(exprSet['1939_at',]~pd$Disease)  #s这个看起来差异最大 为了避免这种多个探针对应同一个基因的情况 后面会有一个算法 分别求均值 然后保留均值最大的探针
boxplot(exprSet['1974_s_at',]~pd$Disease)
boxplot(exprSet['31618_at',]~pd$Disease)

image.png

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

rm(list = ls())
options(stringsAsFactors = F)
a=read.table('e4-plot.txt',sep = '\t',fill = T,header = T)
# 在read.table（）中 fill= 如果为T且非所有的行中变量数目相同则用空白填补
# 从一个电子表格中导出的文件通常会把拖尾的空字段忽略掉。为了读取这样的文件，必须设置参数 fill = TRUE。
colnames(a)=c('id','subtype','expression','mut')
dat=a
library(ggstatsplot)
ggbetweenstats(data =dat, x = subtype,  y = expression)

image.png

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

rm(list = ls())
options(stringsAsFactors = F)
a=read.table("e5.csv",header = T, sep = ",",fill = T)
colnames(a)
dat=a
library(ggstatsplot)
ggbetweenstats(data=dat,x=Group,y=Expression) #分组查看表达情况

library(ggplot2)
library(survival)
library(survminer)
table(dat$Status)
dat$Status =ifelse(dat$Status=="Dead",1,0)
sfit <-survfit(Surv(Days,Status)~Group,data=dat)
sfit
summary(sfit)
ggsurvplot(sfit,conf.int = F,pval = T)
ggsave("survival_TP53_in_BRCA_TCGA.png")

## 先定义的两种颜色 可以根据RGB的配色表自行调整
ggsurvplot(sfit,palette = c("#0000FF","#FF0000"),
           risk.table=TRUE,pval=TRUE,
           conf.int=TRUE,xlab="Time in monthes",
           ggtheme=theme_light(),
           ncensor.plot=TRUE)
ggsave("2.png")

image.png

下载数据集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

rm(list = ls())
options(stringsAsFactors = F)
#注意查看下载文件大小，检查数据
f='GSE17215_eSet.Rdata'

library(GEOquery)
#这个包注意两个配制，一般来说自动化的配置是足够的
#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)
dim(dat) #[1] 22277     6


##注释基因ID
library(hgu133a.db)
ids=toTable(hgu133aSYMBOL)
head(ids)
dat=dat[ids$probe_id,]
dim(dat)  # [1] 19811     6  dat只留下了有注释探针的，观测数从222277--19811
dat[1:4,1:4]
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)  #  [1] 12402     6  因为存在多个探针对应一个基因的情况因此去除重复的 保留均值最大的一个

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% row.names(dat))
ng=ng[ng %in% row.names(dat)]
dat=dat[ng,]
dim(dat)  # [1] 41  6  dat只包括上述的基因
dat=log2(dat)
pheatmap::pheatmap(dat,,scale = 'row')

image.png

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

rm(list = ls())
options(stringsAsFactors = F)
#注意查看下载文件大小，检查数据
f='GSE24673_eSet.Rdata'

library(GEOquery)
#这个包注意两个配制，一般来说自动化的配置是足够的
#Setting options('download.file.method.GEOquery'='auto')
#Setting options('GEOquery.inmemory.gpl'=FALSE)
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) # 提取表达矩阵
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)  ##相关性
tmp=data.frame(g=group_list)
rownames(tmp)=colnames(M)
pheatmap::pheatmap(M,annotation = tmp)

image.png

找到 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/"))
##根据GPL号找到bioc_pachage号码 然后" .db"下载
BiocManager::install("hugene10sttranscriptcluster.db",ask = F,update = F) 
options()$repos
options()$BioC_mirror

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

rm(list = ls())
options(stringsAsFactors = F)
#注意查看下载文件大小，检查数据
f='GSE42872_eSet.Rdata'

library(GEOquery)
#这个包注意两个配制，一般来说自动化的配置是足够的
#Setting options('download.file.method.GEOquery'='auto')
#Setting options('GEOquery.inmemory.gpl'=FALSE)
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) ##分组信息
##(平均表达量/sd/mad/)最大的探针
boxplot(dat)
sort(apply(dat,1,mean),decreasing = T)[1]  #apply函数 对每一行求men/sd/mad 然后降序排列 并取出第一个值 即最大值
sort(apply(dat,1,sd),decreasing = T)[1]
sort(apply(dat,1,mad),decreasing = T)[1]

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

rm(list = ls())
options(stringsAsFactors = F)
#注意查看下载文件大小，检查数据
f='GSE42872_eSet.Rdata'
library(GEOquery)
#这个包注意两个配制，一般来说自动化的配置是足够的
#Setting options('download.file.method.GEOquery'='auto')
#Setting options('GEOquery.inmemory.gpl'=FALSE)
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) ##分组信息
boxplot(dat)
##提取分组名 title以空格分的第四个元素
group_list=unlist(lapply(pd$title,function(x){
strsplit(x,' ')[[1]][4]
}))

exprSet=dat
exprSet[1:4,1:4] #查看数据 limma接收到额一定是log之后的矩阵

##DEG by  limma
suppressMessages(library(limma))
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(exprSet)
design
contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = "-"),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)

还需要加上ID转换才能用

【生信技能树】R语言中级练习题

猜你喜欢

热点阅读