2020.2.6 R语言|Practice4

作者: 哈喽迷人鬼们 | 来源:发表于2020-02-07 20:26 被阅读0次

2020.2.6 R语言|Practice4
学习小组Day4笔记-皇晓燕
Day4 学习小组--张小张
学习小组Day4笔记--扬马延
《学习小组Day4笔记--寒鹤》
R语言基础--数据类型-总结
运用SPECS来灵修：路15
学习小组Day4-沈荣
R 语言脚本操作MySQL 在Windows环境的自动运行
学习小组day4笔记-曾俊辉

2020.2.7 R语言|Practice4

今日学习内容一览

一、Jimmy 生信技能树B站生信小技巧系列课程01——05讲

二、完成所有R语言小作业（中级）

三、关于近期学习R值得巩固基础的几个方面

二、R语言小作业（中级）

找到BRCA1基因在TCGA数据库的[乳腺癌数据集](Breast Invasive Carcinoma (TCGA, PanCancer Atlas))的表达情况——差异分析

rm(list = ls()) 
options(stringsAsFactors = F)
a=read.table('BRC1在乳腺癌中的表达情况.txt', sep = '\t',fill = T,header = T)
colnames(a)=c('id','subtype','expression','mut')
dat=a

#BiocManager::install('ggstasplot')
#BiocManager::install('metaBMA')

library(ggstatsplot)
ggbetweenstats(data=dat,x=subtype,y=expression)
library(ggplot2)
ggsave('plot-BRC1-TCGA.png')

（1）在cbioportal网站上下载BRCA1基因在乳腺癌中的TCGA数据库

（2）根据情况改变表达矩阵的列名（'id','subtype','expression','mut'）

（3）使用ggstatsplot R包中的ggbetweenstats函数，以data为源数据，以subtype为x轴，以expression为y轴画图

（4）使用ggplot2 R包中的ggsave函数可存储源图像

plot-BRC1-TCGA.png

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

rm(list = ls()) 
options(stringsAsFactors = F)
a=read.table('53基因.csv', sep = ',',fill = T,header = T)
dat=a
library(ggplot2)
library(survival)
#BiocManager::install('survminer')
library(survminer)
#BiocManager::install('survfit')
#生存分析代码
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=TRUE)
ggsave('survuval_TP53_in_BRCA_TCGA.png')

（1）在UCSC.xena网站上下载BRCA1基因在乳腺癌中的TCGA数据库并读入R

（2）使用ggplot2、survival、survminer 三个R包对表达矩阵进行分组（Dead or Alive/Low or High）

（3）检查源数据表达矩阵中是否含有Days、Status、Expression信息

（4）使用ggsurvplot函数进行绘图并使用ggsave函数保存

survuval_TP53_in_BRCA_TCGA.png

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

（1）下载GSE17215数据集（下载GSE类数据集均使用此代码，仅改变序号）

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

f='GSE17215_eSet.Rdata'

library(GEOquery)

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]])
a=gset[[1]] 
dat=exprs(a)
dim(dat)

（2）探针ID和基因名的相互转化

ids=toTable(hgu133aSYMBOL)
head(ids)
dat=dat[ids$probe_id,]
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)

（3）筛选出仅含有目的基因的表达矩阵并绘制热图

ng=strsplit(ng,' ')[[1]]
ng %in% rownames(dat)
table(ng %in% rownames(dat))
ng=ng[ng %in% rownames(dat)] ##过滤掉ng不在dat里面的基因名
dat[ng,]
dat=dat[ng,]
dat=log2(dat)
pheatmap::pheatmap(dat,scale='row')

作业6 热图.png

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

（1）下载GSE17215数据集（下载GSE类数据集均使用此代码，仅改变序号）

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]])
a=gset[[1]] 
dat=exprs(a)
dim(dat)

（2）使用pData函数获取临床分组信息，再利用group_list自制分组类型

pd=pData(a)
group_list=c('rbc','rbc','rbc',
             'rbn','rbn','rbn',
             'rbc','rbc','rbc',
             "normal",'normal')

（3）使用cor函数计算列与列之间的相关系数后，向热图中加入分组信息绘图

dat[1:4,1:4]
M=cor(dat)
pheatmap::pheatmap(M)
tmp=data.frame(g=group_list)
row.names(tmp)=colnames(dat)
pheatmap::pheatmap(M,annotation_col = tmp)

作业7热图.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/"))
#所有的注释包都要加上.db后缀
BiocManager::install("你需要下载的bioconductor注释包.db",ask = F,update = F)
options()$repos
options()$BioC_mirror

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

（1）下载GSE42872数据集（下载GSE类数据集均使用此代码，仅改变序号）

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]])
a=gset[[1]] 
dat=exprs(a)
dim(dat)

作业9.png

（2）使用apply函数取每一行的平均值，再使用sort函数进行排序选第一位

sort(apply(dat,1,mean),decreasing = T)[1]
boxplot(dat)
sort(apply(dat,1,sd),decreasing = T)[1]
sort(apply(dat,1,mad),decreasing = T)[1]

作业99.png

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

（1）下载GSE42872数据集（下载GSE类数据集均使用此代码，仅改变序号）

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]])
a=gset[[1]] 
dat=exprs(a)
dim(dat)

（2）为exprSet赋值并检验其是否需要log

exprSet=dat
exprSet[1:4,1:4]

（3）使用limma R包对表达矩阵进行差异分析（较为固定的代码块）

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<-makeContrasts("progres.-stable",levels = design)
#最重要的一步
contrast.matrix

fit <- lmFit(exprSet,design)
fit2 <- contrasts.fit(fit,contrast.matrix)
fit2 <- eBayes(fit2) 
tempOutput = topTable(fit2,coef=1,n=Inf)
nrDEG=na.omit(tempOutput)
head(nrDEG)