R语言小作业-中级

学习了一段时间的R，开始做些题，来加深所学的知识。

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

1.新建TXT，将题中基因复制到TXT，并命名。然后导入到R。

a<-data.table::fread("1.txt",header=F)

2.导入R后发现得进行小数点分割，由于lappy和str_split比较复杂，因此采用separate函数，

a2 <- separate(a,V1,into = "ensembl_id",sep = "[.]")

3.分割后剩下的就是和参考基因组进行对接

ibrary(org.Hs.eg.db)
g2s=toTable(org.Hs.egSYMBOL)
g2e=toTable(org.Hs.egENSEMBL)
c<-merge(g2e,a2,by="ensembl_id")
d<-merge(c,g2s,by="gene_id")

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

思路和上一题一样

a<-data.table::fread("1.txt",header=F)
colnames(a)<-c("probe_id")
d<-merge(a,ids,by="probe_id")

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

此题是应了解什么是expressionSet,

rm(list=ls())
suppressPackageStartupMessages(library(CLL))
data(sCLLex)
sCLLex
exprSet=exprs(sCLLex) 
library(hgu95av2.db)
pd=pData(sCLLex)#查看对像里的分组信息
ids=toTable(hgu95av2SYMBOL)
head(ids)
boxplot(exprSet["1939_at",]~pd$Disease)

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

cbioportal地址：http://www.cbioportal.org/
1.按照网址下载所需要的数据plot.txt
2.将下载的数据导入到R，进行整理

a <- data.table::fread("plot.txt")
colnames(a) = c('id','subtype','expression','mutant')
suppressPackageStartupMessages(library(ggstatsplot))
ggbetweenstats(data = BRCA1, x="subtype", y="expression")

image.png

5.找到TP53基因在TCGA数据库的乳腺癌数据集的表达量分组看其是否影响生存提示使用：http://www.oncolnc.org/

按照B站jimmy大神的操作找到源数据文件，导入进R后作图

a<-data.table::fread("BRCA_7157_50_50.csv")
class(a)
View(a)
ggbetweenstats(data = a, x='Status', y='Expression')

image.png

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

suppressPackageStartupMessages(library(GEOquery))
gse17215 = getGEO("GSE17215", destdir=".", AnnotGPL = F, getGPL = F) 
class(gse17215)
length(gse17215)
class(gse17215[[1]])
a=gse17215[[1]]
dat<-exprs(a)
class(dat)
dim(dat)
head(dat)
library(hgu133a.db)
ids=toTable(hgu133aSYMBOL)
colnames(data)="probe_id"
ids$median=apply(dat,1,median)
dat=dat[ids$probe_id,]
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
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')

image.png

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

gse24673 = getGEO("GSE24673", destdir=".", AnnotGPL = F, getGPL = F)
class(gse24673)
a=gse24673[[1]]
dat=exprs(a)
class(a)
dim(dat)
pd=pData(a)
group_list=c('rbc','rbc','rbc',
             'rbn','rbn','rbn',
             'rbc','rbc','rbc',
             'normal','normal')
# 相关性分析并作图
M=cor(dat)
View(M)
pheatmap::pheatmap(M)

image.png

tmp=data.frame(g=group_list)
class(tmp)
View(tmp)
rownames(tmp)=colnames(M)
pheatmap::pheatmap(M,annotation_col = tmp)

image.png

8.找到 GPL6244 platform of Affymetrix Human Gene 1.0 ST Array 对应的R的bioconductor注释包，并且安装它！

hugene10sttranscriptcluster
BiocManager::install("hugene10sttranscriptcluster.db")

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

表达矩阵的下载已经看过很多次了，但是还是忘
rm(list = ls())  
options(stringsAsFactors = F)
library(GEOquery)
 gset <- getGEO('GSE42872', destdir=".",
                 AnnotGPL = F,     
                 getGPL = F)      
a=gset[[1]]
dat=exprs(a)
pd=pData(a)
# 分别挑选出 所有样本的(平均表达量/sd/mad/)最大的探针
sort(apply(dat,1,mean),decreasing = T)[1]
sort(apply(dat,1,sd),decreasing = T)[1]
sort(apply(dat,1,mad),decreasing = T)[1]
# 加载bioconduct注释包
library(hugene10sttranscriptcluster.db)
ls("package:hugene10sttranscriptcluster.db")
ggg <-toTable(hugene10sttranscriptclusterSYMBOL)
# 找一下；尴尬的是7978905找不到...
ggg2 <- ggg[ggg$probe_id%in%7978905,]
ggg2 <- ggg[ggg$probe_id%in%8133876,]
后来发现基因注释包和表达矩阵就没有相同的探针

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

rm(list = ls())  
options(stringsAsFactors = F)
library(GEOquery)
gset <- getGEO('GSE42872', destdir=".",
               AnnotGPL = F,     
               getGPL = F)      
class(gset)
a<-gset[[1]]
b<-exprs(a)
pd=pData(a)
library(stringr)
group_list=unlist(lapply(pd$title,function(x){
  str_split(x,' ')[[1]][4]
}))
exprSet=b

suppressMessages(library(limma)) 
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(exprSet)

contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = "-"),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)

image.png