生信R語言-初探

1、根據(jù)R包org.Hs.eg.db找到ensembol 基因ID對應(yīng)的基因名(symbol)

library(org.Hs.eg.db)
if(!require(stringr))install.packages('stringr')
library(stringr)
g2s=toTable(org.Hs.egSYMBOL)
g2e=toTable(org.Hs.egENSEMBL)
g <- data.frame(id=c("ENSG00000000003.13","ENSG00000000005.5","ENSG00000000419.11","ENSG00000000457.12","ENSG00000000460.15","ENSG00000000938.11"))
a <- unlist(lapply(g$id,function(x){
  x
  str_split(x,'[.]')[[1]][1]}))
geneid <- g2e[match(a,g2e$ensembl_id),]
b <- g2s[match(geneid$gene_id,g2s$gene_id),]
merge(geneid,b,by.x='gene_id',by.y='gene_id')

2、根據(jù)R包hgu133a.db找到探針對應(yīng)的基因名(symbol)

rm(list = ls())
library(hgu133a.db)
ids=toTable(hgu133aSYMBOL)
head(ids)
tz <- c("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")
ids[match(tz,ids$probe_id),]

3、找到R包CLL內(nèi)置的數(shù)據(jù)集的表達(dá)矩陣?yán)锩娴腡P53基因的表達(dá)量，并且繪制在 progress-stable分組的boxplot圖,通過 ggpubr 進行美化

rm(list = ls())
suppressPackageStartupMessages(library(CLL))
data(sCLLex)
pd=pData(sCLLex)
exprSet=exprs(sCLLex) 
library(hgu95av2.db)
ids <- toTable(hgu95av2SYMBOL)
tmp=ids[ids$symbol=='TP53',1]

for(i in 1:length(tmp)){
  boxplot(exprSet[i,]~pd$Disease)
}

library(ggpubr)
lapply(1:length(tmp),function(i){
  a <- cbind(exprSet[i,],pd)
  colnames(a)[1]='type'
  ggboxplot(data = a,x = 'Disease',y = 'type',color = 'Disease')
})
image.png

image.png

image.png

image.png

image.png

image.png

4、找到BRCA1基因在TCGA數(shù)據(jù)庫的乳腺癌數(shù)據(jù)集(Breast Invasive Carcinoma (TCGA, PanCancer Atlas))的表達(dá)情況

rm(list = ls())
options(stringsAsFactors = F)
a <- read.table('plot.txt',sep = '\t',header = T,fill = T)
library(ggplot2)
ggplot(a,aes(x=a$Subtype,y = a$BRCA1..mRNA.Expression.Zscores..RSEM..Batch.normalized.from.Illumina.HiSeq_RNASeqV2.,col=a$Subtype))+geom_boxplot()+geom_jitter()+theme_bw()
image.png

5、找到TP53基因在TCGA數(shù)據(jù)庫的乳腺癌數(shù)據(jù)集的表達(dá)量分組看其是否影響生存

rm(list = ls())
options(stringsAsFactors = F)
a <- read.csv('BRCA_7157_50_50.csv',header = T,fill = T)
library(ggplot2)
library(survival)
library(survminer)
dat=a
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)
ggsurvplot(sfit,palette = c('#E7B800','#2E9FDF'),risk.table = TRUE,pval = TRUE,conf.int = TRUE,xlab='Time in months',ggtheme = theme_light())
image.png

image.png

6、下載數(shù)據(jù)集GSE17215的表達(dá)矩陣并且提取下面的基因畫熱圖。

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)
#install.packages('pheatmap')
library(GEOquery)
f='GSE17215_eSet.Rdata'
if(!file.exists(f)){
  gset <- getGEO('GSE17215',destdir='.',
                 AnnotGPL = F,
                 getGPL=F)
  save(gset,file = f)
}
a <- strsplit('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',' ')[[1]]
library(hgu133a.db)
ids=toTable(hgu133aSYMBOL)
#head(ids)
tz <- ids[match(a,ids$symbol,nomatch = 0),]
load(f)
class(gset)
dat=exprs(gset[[1]])
dim(dat)
dat2 <- dat[tz$probe_id,]
write.csv(dat2,'dat2.csv')
dat2=read.csv('dat2.csv',header = T)
colnames(dat2)[1]='pid'
dat3=merge(dat2,tz,by.x='pid',by.y='probe_id')
dat3=dat3[,-1]
rownames(dat3) <- dat3$symbol
dat3=dat3[,-ncol(dat3)]
dat3=log2(dat3)
library(pheatmap)
pheatmap::pheatmap(dat3)
image.png

7、下載數(shù)據(jù)集GSE24673的表達(dá)矩陣計算樣本的相關(guān)性并且繪制熱圖，需要標(biāo)記上樣本分組信息

rm(list = ls())
options(stringsAsFactors = F)
library(GEOquery)
f='GSE24673_eSet.Rdata'
if(!file.exists(f)){
  gset <- getGEO('GSE24673',destdir='.',
                 AnnotGPL = F,
                 getGPL=F)
  save(gset,file = f)
}
load(f)
class(gset)
dat=exprs(gset[[1]])
pd=pData(gset[[1]])
group_list=c(rep('rbc',3),rep('rbn',3),rep('rbc',3),rep('normal',2))
M=cor(dat)
pheatmap::pheatmap(M)
tmp=data.frame(g=group_list)
rownames(tmp)=colnames(M)
pheatmap::pheatmap(M,annotation_col = tmp)
image.png

image.png

8、找到 GPL6244 platform of Affymetrix Human Gene 1.0 ST Array 對應(yīng)的R的bioconductor注釋包，并且安裝它！

rm(list = ls())
options(stringsAsFactors = F)
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、下載數(shù)據(jù)集GSE42872的表達(dá)矩陣，并且分別挑選出 所有樣本的(平均表達(dá)量/sd/mad/)最大的探針，并且找到它們對應(yīng)的基因。

rm(list = ls())
options(stringsAsFactors = F)
library(GEOquery)
f='GSE42872_eSet.Rdata'
if(!file.exists(f)){
  gset <- getGEO('GSE42872',destdir='.',
                 AnnotGPL = F,
                 getGPL=F)
  save(gset,file = f)
}
load(f)
a=gset[[1]]
dat=exprs(a)
dim(dat)
pd=pData(a)
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]
image.png

10、下載數(shù)據(jù)集GSE42872的表達(dá)矩陣，并且根據(jù)分組使用limma做差異分析，得到差異結(jié)果矩陣

rm(list = ls())
options(stringsAsFactors = F)
library(GEOquery)
f='GSE42872_eSet.Rdata'
if(!file.exists(f)){
  gset <- getGEO('GSE42872',destdir='.',
                 AnnotGPL = F,
                 getGPL=F)
  save(gset,file = f)
}
load(f)
a=gset[[1]]
dat=exprs(a)
dim(dat)
pd=pData(a)
group_list=unlist(lapply(pd$title,function(x){strsplit(x,' ')[[1]][4]}))
exprSet=dat
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 
##step1
fit <- lmFit(exprSet,design)
##step2
fit2 <- contrasts.fit(fit, contrast.matrix) ##這一步很重要，大家可以自行看看效果
fit2 <- eBayes(fit2,trend=TRUE)  ## 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)

不懂的問題歡迎來問小陌哦~

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容