R語言小作業(yè)-中級
代碼來自Jimmy的Github 所以更多的精力是理解注釋代碼。
- 根據(jù)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)
#去掉小數(shù)點之后的數(shù)字來提取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')
- 根據(jù)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)
# 根據(jù)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內置的數(shù)據(jù)集的表達矩陣里面的TP53基因的表達量，并且繪制在 progres.-stable分組的boxplot圖
rm(list = ls())
options(stringsAsFactors = F)
suppressPackageStartupMessages(library(CLL)) # 靜默加載功炮，不會出現(xiàn)下面的紅字
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數(shù)據(jù)庫的乳腺癌數(shù)據(jù)集(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且非所有的行中變量數(shù)目相同則用空白填補
# 從一個電子表格中導出的文件通常會把拖尾的空字段忽略掉蹦误。為了讀取這樣的文件芝薇，必須設置參數(shù) fill = TRUE胚嘲。
colnames(a)=c('id','subtype','expression','mut')
dat=a
library(ggstatsplot)
ggbetweenstats(data =dat, x = subtype, y = expression)
image.png
- 找到TP53基因在TCGA數(shù)據(jù)庫的乳腺癌數(shù)據(jù)集的表達量分組看其是否影響生存
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")
## 先定義的兩種顏色 可以根據(jù)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
image.png
image.png
- 下載數(shù)據(jù)集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)
#注意查看下載文件大小，檢查數(shù)據(jù)
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') ##載入數(shù)據(jù)
class(gset)
length(gset)
class(gset[[1]])
#因為這個GEO數(shù)據(jù)集只有一個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只留下了有注釋探針的，觀測數(shù)從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
- 下載數(shù)據(jù)集GSE24673的表達矩陣計算樣本的相關性并且繪制熱圖晾嘶，需要標記上樣本分組信息
rm(list = ls())
options(stringsAsFactors = F)
#注意查看下載文件大小妓雾，檢查數(shù)據(jù)
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') ##載入數(shù)據(jù)
class(gset)
length(gset)
class(gset[[1]])
#因為這個GEO數(shù)據(jù)集只有一個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/"))
##根據(jù)GPL號找到bioc_pachage號碼 然后" .db"下載
BiocManager::install("hugene10sttranscriptcluster.db",ask = F,update = F)
options()$repos
options()$BioC_mirror
- 下載數(shù)據(jù)集GSE42872的表達矩陣机断，并且分別挑選出 所有樣本的(平均表達量/sd/mad/)最大的探針楷拳，并且找到它們對應的基因材部。
rm(list = ls())
options(stringsAsFactors = F)
#注意查看下載文件大小，檢查數(shù)據(jù)
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') ##載入數(shù)據(jù)
class(gset)
length(gset)
class(gset[[1]])
#因為這個GEO數(shù)據(jù)集只有一個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函數(shù) 對每一行求men/sd/mad 然后降序排列 并取出第一個值 即最大值
sort(apply(dat,1,sd),decreasing = T)[1]
sort(apply(dat,1,mad),decreasing = T)[1]
- 下載數(shù)據(jù)集GSE42872的表達矩陣苦丁，并且根據(jù)分組使用limma做差異分析浸颓，得到差異結果矩陣
rm(list = ls())
options(stringsAsFactors = F)
#注意查看下載文件大小，檢查數(shù)據(jù)
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') ##載入數(shù)據(jù)
class(gset)
length(gset)
class(gset[[1]])
#因為這個GEO數(shù)據(jù)集只有一個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] #查看數(shù)據(jù) 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轉換才能用
