[ZHUHAI_Biotrainee]R语言中级作业(全)

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

开始看之前，要做好准备！这篇简书有些长，适合初学者的有充足的时间运行和比较代码看，如果时间有限或者性格比较急的建议将该篇简书先放一放，感谢查阅！

这道题主要用的是“org.Hs.eg.db”,和“stringr"先来了解一下这个包哦！

开始做题啦！
铺垫性质的~第一步：
读入外部文件——为后续的merge做准备，对外部数据进行处理（Eg:发现在g2e中ensembl_id是整数，没有后面的小数部分 ）

rm(list=ls())

a <- read.table("ensembl.txt")
library(stringr)
colnames(a)
str_split(a$V1,"[.]")
[[1]]
[1] "ENSG00000000003" "13"             

[[2]]
[1] "ENSG00000000005" "5"              

[[3]]
[1] "ENSG00000000419" "11"             

[[4]]
[1] "ENSG00000000457" "12"             

[[5]]
[1] "ENSG00000000460" "15"             

[[6]]
[1] "ENSG00000000938" "11"  
unlist(str_split(a$V1,"[.]"))
[1] "ENSG00000000003" "13"              "ENSG00000000005" "5"               "ENSG00000000419"
 [6] "11"              "ENSG00000000457" "12"              "ENSG00000000460" "15"             
[11] "ENSG00000000938" "11"
class(unlist(str_split(a$V1,"[.]")))
unlist(str_split(a$V1,'[.]',simplify = T))    #可以查看？str_split的说明书，simplify=T，返回的是一个矩阵，我们需要的就是一个矩阵的形式 更容易处理
unlist(str_split(a$V1,'[.]',simplify = T))[,1] 
a$ensembl_id <- unlist(str_split(a$V1,'[.]',simplify = T)[,1])
# 这个代码是晓娟老师给的，我没看懂  所以我用的是上面的代码    a<-data.frame(strsplit(gsub('\\.\\d+','',a),'\\s+'))

心得：
①关于str_split函数：在这个代码中，我们可以看到用str_split函数对a$V1进行字符串的分割，我们发现结果是一个list，而list就比较尴尬了，目前我是没有办法对list批量取出一列。所以这时候我们不禁想到如果可以把刚刚的这种list转化成matrix，就十分方便了。敲黑板了！！！这里我认为是很方便的一个参数就是simplify，他是str_split函数的一个参数，可以把character转化成matrix。
②关于data.frame:再次重复data.frame的结构特点。
a.数据框是非常整齐的一个数据结构类型
b.数据框的每一列/每一行的数据结构类型都是相同的。比如每一列/每一行都是char/numeric/logical等。
③如果你看完了②并运行了以上代码，通过思考你就会发现这里的unlist在这个部分的利用率并不高。只是更加方便去理解data.frame
④关于merge以及合并/加一列的相关问题：反复强调的merge是关联有相同的列或者行名时才会用到的一个函数。我在[ZHUHAI_BiotraineeR语言初级题目（下）]中就犯过类似的错误，如果不理解，翻开看刚刚提到的简书！总结如下：
a.merge是关联有相同的列或者行名
b.如果想合并两个既没有相同列/行，又没有相同内容的列/行的data.frame时，不可以用merge。 那么合并data.frame怎么办呢？ EG：data.frame(data1,data2) ######OK，so easy！
c.同b比较类似，如果想在数据框中加一列，怎么办？EG：在a的这个数据框中加ensembl_id这个字符串名字，并且加的这个名字一定要等于一个具备构成数据框的结构

接着c点说，举个例子：

tmp <- str_split(a$V1,"[.]",simplify = T)[,1]       
data.frame(tmp,a) 

上面的这个代码等同于下面的这一个：

a$ensembl_id <- str_split(a$V1,"[.]",simplify = T)[,1]

前面洋洋洒洒的写了那么多，主要是为了深入去理解基础，真正的第一步代码请运行如下代码：

a <- read.table("ensembl.txt")
library(stringr)
colnames(a)
str_split(a$V1,"[.]")
str_split(a$V1,"[.]",simplify = T)
str_split(a$V1,"[.]",simplify = T)[,1]
a$ensembl_id <- str_split(a$V1,"[.]",simplify = T)[,1]

第二步：
关联：merge()

library("org.Hs.eg.db")
keytypes(org.Hs.eg.db)
g2s <- toTable(org.Hs.egSYMBOL)
g2e <- toTable(org.Hs.egENSEMBL)
b <- merge(a,g2e,by="ensembl_id",all=T)
d <- merge(b,g2s,by="gene_id",all=T)

按照题目要求，第一题结束。but,我们关联了两次，观察结果我们发现出现很多次的NA。So,追求完美的我们需要去重，调序

##搜索 d 中 ensembl_id 重复的基因
table(table(d$ensembl_id)>1)
table(d$ensembl_id)[table(d$ensembl_id)>1]
#去重，调序
d <- d[order(d$V1),]
d <- d[!duplicated(d$V1),]
d <- [match(a$V1,d$V1),]

认为有点难度，先给代码！

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

rm(list = ls())
a1 <- read.table("ensembl_1.txt")
a1
library(stringr)
colnames(a1)
str_split(a1$V1,"[_]")
str_split(a1$V1,"[_]",simplify = T)
str_split(a1$V1,"[_]",simplify = T)[,1]
a1$probe_id <- str_split(a1$V1,"[_]",simplify = T)[,1]  

之前我在这一步看了一下a1,感觉多出来了第一列V1，于是我进行了a1 <- a1[,-1],但是删除了之后a1变成了char，因为一列。so,我们不应该删掉第一列，让a1孤零零的成为了char

#探索symbol
library(hgu133a.db)
keytypes(hgu133a.db)
ids=toTable(hgu133aSYMBOL)
head(ids)
b1 <- merge.data.frame(a1,ids,by="probe_id",all=T) 

第二题与第一题的比较：
①两道题非常相似，只是分隔符的差别。
②如果在用str_split()函数的时候，我们只用"" 不加"[]"就只是字符隔。
Eg:第一题：应该用str_split(aV1,".")，达不到要求。

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

rm(list = ls())
library(affy)
library(CLL)
data(sCLLex)
exprSet=exprs(sCLLex)
pd=pData(sCLLex)         #pData什么意思?
library(hgu95av2.db) 
ids=toTable(hgu95av2SYMBOL)
head(ids)
boxplot(exprSet['1939_at',] ~ pd$Disease) ## sig
boxplot(exprSet['1974_s_at',] ~ pd$Disease)
boxplot(exprSet['31618_at',] ~ pd$Disease)
#ggpubr
d<-cbind(as.data.frame(exprSet['1939_at',]),as.data.frame(pd$Disease))
head(d)
library(ggpubr)
colnames(d)<-c("e","f")
p<-ggboxplot(d, x=("f"), y =("e") ,
             color = "f",
             palette = "jco",
             add = "jitter")
opar<-par(no.readonly=T)
#  Add p-value
p + stat_compare_means()
# Change method
p + stat_compare_means(method = "t.test")

4、找到BRCA1基因在TCGA数据库的乳腺癌数据集的表达情况

get到的新技能，请看:
A、学到了一个新的网址: cBioPortal http://www.cbioportal.org/datasets （如果不会用，网页中有一些教程，可以参考）
B、过程如下：

第一步
第二步
第三步
_______________________________Begin_________________________________

方法一：

rm(list = ls())
options(stringsAsFactors = F)
a=read.table('e4-plot.txt',sep = '\t',fill = T,header = T)
View(a)
colnames(a)=c('id','subtype','expression','mut')
dat=a
library(ggstatsplot)
ggbetweenstats(data =dat, x = subtype,  y = expression)
library(ggplot2)
ggsave('plot-again-BRCA1-TCGA-BRCA-cbioportal.png')

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

方法一：

rm(list=ls())
options(stringsAsFactors = F)
a=read.table('BRCA_7157_50_50.csv',header = T,sep = ',',fill = T)
library(ggstatsplot)
ggbetweenstats(data = a,x=Group,y=Expression)
library(ggplot2)
library(survival)
library(survminer)
table(a$Status)
a$Status=ifelse(a$Status=='Dead',1,0)
sfit <- survfit(Surv(Days,Status)~Group,data=a)
sfit
summary(sfit)
ggsurvplot(sfit,conf.int = F,pval = TRUE)

方法二：

rm(list = ls())
options(stringsAsFactors = F)
a=read.table('BRCA_7157_50_50.csv',sep = ',',fill = T,header = T)
#View(a)
dat=a
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=TRUE)
ggsave('survival_TP53_in_BRCA_TCGA.png')

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

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)


library(hgu133a.db)
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)

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的表达矩阵计算样本的相关性并且绘制热图，需要标记上样本分组信息

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]])   #在这一步时发现gest是一个list的格式，如果我们不把[1]取出，后面就没办法对列表进行操作，并且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)
pheatmap::pheatmap(M)
tmp=data.frame(g=group_list)
rownames(tmp)=colnames(M)
pheatmap::pheatmap(M,annotation_col = tmp)

关于第6、7题的心得：
首先，通过做题我们发现第六题和第七题非常相似。
①在第6、7中都要if所执行的代码，那么if所运行的代码是什么意思呢？
这一步的主要作用根据我们了解的GEO的登录信息、平台信息去在RStudio中下载GEO的压缩文件，也就是说我们有了这一步就不需要再自行下载GEO的压缩文件了。
②load('GSE24673_eSet.Rdata') 的命令主要是在压缩文件中解压出为“Rdata”格式的文件
③class(gset[[1]]) 在这一步时发现gest是一个list的格式，如果我们不把[1]取出，后面就没办法对列表进行操作，并且gset[[1]]里面有很多内置的数据信息
③exprs表示获取a的表达矩阵,为后续的cor做准备
④cor(x)函数为做相关性的函数，如果x不是矩阵或者向量，我们就没办法做出相关性。
⑤pheatmap与heatmap 两个包是不一样的，利用这两个包画出来的图也是完全不同的。so,一定要明确你的包呀~ 不得不说我之前用heatmap的包画的图巨丑。。

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、下载数据集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]
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)
# 这个包需要注意两个配置，一般来说自动化的配置是足够的。
#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)
group_list=unlist(lapply(pd$title,function(x){
  strsplit(x,' ')[[1]][4]
}))


exprSet=dat
exprSet[1:4,1:4]
# 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<-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)