生信技能樹(shù)聯(lián)盟入門視頻：生信技能樹(shù)視頻課程學(xué)習(xí)路徑，這么好的視頻還免費(fèi)！
專題歷史目錄：3個(gè)學(xué)生的linux視頻學(xué)習(xí)筆記
接下來(lái)介紹R語(yǔ)言：
【生信技能樹(shù)】生信人應(yīng)該這樣學(xué)R語(yǔ)言
R語(yǔ)言
在你開(kāi)始R之旅前，建議你看看下面這兩個(gè)
1. 介紹R語(yǔ)言及Rstudio
- 了解R，Rstudio及R包;安裝的包在packages中檢查
.libPaths() #找安裝路徑 - 幫助文檔，幫忙看路徑
?substring - 定位文件蘸嘶、設(shè)置文件位置
getwd()
setwd()
- plot畫板關(guān)閉
dev.off()
2. R語(yǔ)言基礎(chǔ)變量講解
- 重點(diǎn)就是理解: 五種變量結(jié)構(gòu)（class屬性）
- 我也曾經(jīng)寫過(guò)一點(diǎn)這方面筆記：R語(yǔ)言學(xué)習(xí)筆記
- grep()搜索函數(shù)
index1 = grep('RNA-Seq', a$Assay_Type)
index2 = grepl('RNA-Seq', a$Assay_Type)
b = a[index1,] # 下標(biāo)
b = a[index2,] # 索引
3. 外部數(shù)據(jù)導(dǎo)入導(dǎo)出
- Excel表格到R語(yǔ)言轉(zhuǎn)換
- 去GEO數(shù)據(jù)庫(kù)下載表達(dá)矩陣GSE17215
用excel打開(kāi)后
GEO17215
- 此時(shí)用導(dǎo)入R會(huì)出現(xiàn)這種情況
> b=read.table('GSE17215_series_matrix.txt.gz',sep = '\t',comment.char = 'T',header = T)
Error in scan(file = file, what = what, sep = sep, quote = quote, dec = dec, :
line 29 did not have 2 elements
# 因?yàn)榇嬖谝恍￠_(kāi)頭的文件存在陪汽，此時(shí)經(jīng)過(guò)下面這種處理：
# 帶感嘆號(hào)的不讀取
> b=read.table('GSE17215_series_matrix.txt.gz',sep = '\t',comment.char = '!',header = T)
> View(b)
保存训唱，進(jìn)一步處理GSE17215_series_matrix.csv
# 保存為.CSV格式
write.csv(b,'GSE17215_series_matrix.csv')
# 第一列設(shè)為行名，再去掉第一行
rownames(b)=b[,1]
b=b[,-1]
-
經(jīng)過(guò)上面處理后如下掩缓；
image.png - 畫個(gè)熱圖看看看
b = log2(b)
pheatmap::pheatmap(b[1:10,])
熱圖
- 保存b推薦采用下面這種方式進(jìn)行
save(b,file = 'b_input.Rdata')
load(file = 'b_input.Rdata')
- 建議把excel轉(zhuǎn)成csv來(lái)讀取
4. 中級(jí)變量操作
- 所有函數(shù)有參數(shù)雪情，很多是互通的， eg sort, max, min, fivenum(四分位值)
函數(shù)
# 最大值
> sort(b$GSM431121,decreasing = T)[1]
[1] 15.28882
> max(b$GSM431121)
[1] 15.28882
# 最小值你辣，五分位數(shù)
> min(b$GSM431121)
[1] 3.859587
> fivenum(b$GSM431121)
[1] 3.859587 4.710004 5.952603 8.691293 15.288819
向量化操作
> table(b$GSM431121<5)
FALSE TRUE
14398 7879
> d=b[b$GSM431121<5,]
# 你可以看看你的b，是為有7879行
看b的每行的平均值
> mean(b[1,])
[1] NA
Warning message:
In mean.default(b[1, ]) : 參數(shù)不是數(shù)值也不是邏輯值：回覆NA
> mean(as.numeric(b[1,]))
[1] 10.60797
> as.numeric(b[1,])
[1] 8.911691 9.221081 11.410364 11.325483 11.418782 11.360438
> mean(as.numeric(b[1,]))
[1] 10.60797
采用rowMean函數(shù)尘执，head前6行
> head(rowMeans(b))
1007_s_at 1053_at 117_at 121_at 1255_g_at 1294_at
10.607973 7.925899 5.193894 7.168633 4.275652 5.686036
# 采用for循環(huán)看看
for (i in 1:nrow(b)) {
print(mean(as.numeric(b[i,])))
}
for (i in 1:6) {
print(mean(as.numeric(b[i,])))
}
> for (i in 1:6) {
+ print(mean(as.numeric(b[i,])))
+ }
[1] 10.60797
[1] 7.925899
[1] 5.193894
[1] 7.168633
[1] 4.275652
[1] 5.686036
- 使用apply循環(huán)也可以實(shí)現(xiàn)
x=mean(as.numeric(b[1,]))
apply(b,1,function(x){
mean(x)
})
查找最大值
# 最大值的查找
for (i in 1:nrow(b)) {
print(max(as.numeric(b[i,])))
}
apply(b,1,function(y){
max(y)
})
apply(b,1,max)
# 自寫函數(shù)rowMax查找
rowMax=function(z){
apply(z, 1, max)
}
rowMax(b)
計(jì)算每行的方差舍哄，取最大的50個(gè)，畫熱圖
apply(b, 1, sd)
sort(apply(b, 1, sd),decreasing=T)[1:50]
cg=names(sort(apply(b, 1, sd),decreasing=T)[1:50])
pheatmap::pheatmap(b[cg,])
image.png
- 隨機(jī)選取50ge誊锭，sample函數(shù)
sample(1:nrow(b),50)
pheatmap::pheatmap(b[sample(1:nrow(b),50),])
隨機(jī)50個(gè)
+ abs, sqrt :戒對(duì)治表悬，平方根
+ Log, log10, log2, exp: 對(duì)數(shù)與指數(shù)函數(shù)
+ Sin, cos, tan, acos, atan, atan2: 三角函數(shù)
+ sinh, cosh, tanh, asinha, acosh, aranh：雙曲函數(shù)
+ 集合運(yùn)算，reshape, merge總結(jié)
- 思考一下excel表格里面有變量類型嗎
a = read.table('XXtable.txt', head = T
sep = '\t')
b = read.table('GSE17215_series_matrix.txt.gz',
comment,char = '!', head = T,
sep = '\t')
write.csv(b, 'GSE17215_series_matrix.csv')
write.table(b,'tmp.csv', sep = ',')
##把行名去掉
d = read.csv('GSE17215_series_matrix.csv')
# readline 讀入之后拆分
5. 熱圖
- 隨機(jī)產(chǎn)生a1丧靡，并產(chǎn)生熱圖
a1=rnorm(100) #隨機(jī)產(chǎn)生100個(gè)服從正態(tài)分布的數(shù)
?rnorm
dim(a1)=c(5,20) # 添加維度屬性蟆沫，矩陣matrix
pheatmap(a1)
a1
- 產(chǎn)生a2+熱圖
a2=rnorm(100)+2
dim(a2)=c(5,20)
pheatmap(a2)
a2
- 將a1,a2橫向拼接并畫圖
pheatmap(cbind(a1,a2),cluster_cols = F)
a3=cbind(a1,a2) #橫向拼接
a4=rbind(a1,a2) #縱向拼接
a1,a2
6. 選取差異明顯的基因的表達(dá)量矩陣?yán)L制熱圖
rm(list = ls()) #魔幻操作籽暇，一鍵清空
library(pheatmap)
a1 = rnorm(100)
dim(a1) = c(5,20)
pheatmap(a1)
a2 = rnow(100)+2
dim(a2) = c(5,20)
library(pheatmap)
pheatmap(a1, cluster_rows = F, cluster_cols = F)
pheatmap(cbind(a1,a2))
pheatmap(cbind(a1,a2), show_rownames = F, show_colnames = F)
- 拉平極差值
7. ID轉(zhuǎn)換
b=as.data.frame(a3)
paste('a1',1:20,sep = '_')
paste('a2',1:20,sep = '_')
names(b)=c(paste('a1',1:20,sep = '_'),paste('a2',1:20,sep = '_'))
pheatmap(b,cluster_cols = F)
- ENSG基因ID處理
#strsplit('','[.]') #根據(jù)點(diǎn)號(hào)分割
> strsplit('ENSG0000000003.13','[.]')
[[1]]
[1] "ENSG0000000003" "13"
> strsplit('ENSG0000000003.13','[.]')[[1]]
[1] "ENSG0000000003" "13"
> strsplit('ENSG0000000003.13','[.]')[[1]][1]
[1] "ENSG0000000003"
#ID轉(zhuǎn)換包
library(stringr)
a$ensemble_id=str_split(a$v1,'[1]',simplify = T)[,1]
# duplicated() #去重
- 一些包
org.Hs.eg.db #在包里有基因注釋關(guān)系
8. 任意基因任意癌癥表達(dá)量分組的生存分析
- 在線制作生存曲線的工具，oncolnc【http://www.oncolnc.org/】
- 也可以下載原始數(shù)據(jù)自己分析
#讀取csv文件
rm(list = ls()) #清除所有變量
options(stringsAsFactors = F)
a=read.table('LGG_93663_50_50.csv',header = T,sep =',',fill = T)
# 后續(xù)畫圖
colnames(a)
dat=a
# 圖ggbetweenstats
library(ggstatsplot)
ggbetweenstats(data = dat,x=Group,y=Expression)
# 圖ggsurvplot
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_ARHGAP18_in_LGG.png')
ggsurvplot(sfit,palette = c("#E7B800","#2E9FDF"),
risk.table = TRUE,pval = TRUE,
conf.int = TRUE,xlab="Time in months",
ggtheme = theme_light(),
ncensor.plot=TRUE)
#ggsave('survival_ARHGAP18_in_LGG.png')
ggbetweenstats
ggsurvplot
9. 任意基因任意癌癥表達(dá)量和臨床性狀關(guān)聯(lián)
某個(gè)基因在某個(gè)癌癥的表達(dá)量饭庞，關(guān)聯(lián)臨床信息
Cbioportal【http://www.cbioportal.org/】戒悠，同樣下載ARHGAP18
rm(list = ls()) #清除所有變量
options(stringsAsFactors = F)
a=read.table('plot.txt',header = T,sep = '\t',fill = T)
colnames(a)=c('id','stage','gene','mut')
dat=a
library(ggstatsplot)
ggbetweenstats(data = dat,x=stage,y=gene)
image.png
10. 表達(dá)矩陣樣本的相關(guān)性
看兩個(gè)變量的相關(guān)性
> cor(1:10,1:10)
[1] 1
> a = rnorm(10)
> b = rnorm(10)
> cor(a,b)
[1] -0.1555608
> a = rnorm(10)
> b = 10*a+rnorm(10)
> cor(a,b)
[1] 0.9971822
學(xué)習(xí)airway這個(gè)數(shù)據(jù)包
rm(list = ls()) #清除所有變量
options(stringsAsFactors = F)
library(airway)
data("airway") #加載數(shù)據(jù)
exprSet=assay(airway)
colnames(exprSet)
group_list=colData(airway)[,3]
exprSet=exprSet[apply(exprSet,1,function(x) sum（x>1)>5),]
exprSet=log(edgeR::cpm(exprSet)+1)
exprSet=exprSet[names(sort(apply(exprSet, 1,mad),decreasing = T)[1:500])]
學(xué)習(xí)上面得到的exprSet
image.png
> dim(exprSet)
[1] 64102 8
> cor(exprSet[,1],exprSet[,2]) #查看相關(guān)性
[1] 0.9632268
- 可知第一列和第二列的相關(guān)性較好：
- 相關(guān)性真的較好
- 存在較多的0，那么大部分相關(guān)性由這些0決定
直接查看列樣本之間基因表達(dá)的相關(guān)性
> cor(exprSet)
SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517
SRR1039508 1.0000000 0.9632268 0.9431333 0.8978772 0.9476515 0.9014495
SRR1039509 0.9632268 1.0000000 0.9323701 0.9428196 0.9202060 0.9290967
SRR1039512 0.9431333 0.9323701 1.0000000 0.9592993 0.9291853 0.9203183
SRR1039513 0.8978772 0.9428196 0.9592993 1.0000000 0.8834640 0.9498761
SRR1039516 0.9476515 0.9202060 0.9291853 0.8834640 1.0000000 0.9313570
SRR1039517 0.9014495 0.9290967 0.9203183 0.9498761 0.9313570 1.0000000
SRR1039520 0.9603249 0.9398292 0.9827034 0.9413725 0.9298860 0.9030552
SRR1039521 0.9023008 0.9463991 0.9525262 0.9853059 0.8724714 0.9291902
SRR1039520 SRR1039521
SRR1039508 0.9603249 0.9023008
SRR1039509 0.9398292 0.9463991
SRR1039512 0.9827034 0.9525262
SRR1039513 0.9413725 0.9853059
SRR1039516 0.9298860 0.8724714
SRR1039517 0.9030552 0.9291902
SRR1039520 1.0000000 0.9559571
SRR1039521 0.9559571 1.0000000
- 這樣直接查看是沒(méi)有意義的舟山，通常實(shí)驗(yàn)設(shè)計(jì)是分有control組绸狐，實(shí)驗(yàn)組的
-
通過(guò)熱圖看看
image.png - 添加樣本分組信息
exprSet=assay(airway)
colnames(exprSet)
group_list=colData(airway)[,3]
tmp=data.frame(g=group_list)
rownames(tmp)=colnames(exprSet)
pheatmap::pheatmap(cor(exprSet),annotation_col = tmp)
image.png
篩選沒(méi)有意義的基因：表達(dá)量大于1，樣本數(shù)小于5的
- 先看看第一行的情況
> x=exprSet[1,]
> x
SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517 SRR1039520
679 448 873 408 1138 1047 770
SRR1039521
572
> x.1
Error: object 'x.1' not found
> x>1
SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517 SRR1039520
TRUE TRUE TRUE TRUE TRUE TRUE TRUE
SRR1039521
TRUE
> sum(x>1)
[1] 8
# TRUE的邏輯值為1累盗，F(xiàn)LASE的邏輯值為0
- 按照上面的規(guī)則進(jìn)行篩選
> dim(exprSet)
[1] 64102 8
> exprSet=exprSet[apply(exprSet, 1,function(x)sum(x>1)>5),]
> dim(exprSet)
[1] 19481 8
- 進(jìn)一步對(duì)exprSet處理
exprSet=log(edgeR::cpm(exprSet)+1) #edgeR::cpm(exprSet)去除文庫(kù)大小差異
dim(exprSet)
exprSet=exprSet[names(sort(apply(exprSet, 1,mad),decreasing = T)[1:500]),]
dim(exprSet)
M=cor(log2(exprSet+1)) # +1排除0
tmp=data.frame(g=group_list)
rownames(tmp)=colnames(M)
pheatmap::pheatmap(M,annotation_col = tmp)
image.png
view()
dim() #看維度
11. 芯片表達(dá)矩陣下游分析
12. RNA-seq表達(dá)矩陣差異分析
- 芯片的表達(dá)矩陣 exprSet = exprs(sCLLex)
- RNA-seq表達(dá)矩陣 exprSet = assay(ariway) # assay獲取表達(dá)矩陣
rm(list = ls()) #清除所有變量
options(stringsAsFactors = F)
library(airway)
data("airway") #加載數(shù)據(jù)
exprSet=assay(airway)
colnames(exprSet)
group_list=colData(airway)[,3]
exprSet=exprSet[apply(exprSet, 1,function(x)sum(x>1)>5),]
table(group_list)
> table(group_list)
group_list
trt untrt
4 4
exprSet[1:4,1:4]
> exprSet[1:4,1:4]
SRR1039508 SRR1039509 SRR1039512 SRR1039513
ENSG00000000003 679 448 873 408
ENSG00000000419 467 515 621 365
ENSG00000000457 260 211 263 164
ENSG00000000460 60 55 40 35
dev.off()
boxplot(exprSet)
image.png
13. R語(yǔ)言小作業(yè)
image.png
準(zhǔn)備工作--安裝所需的包
cran_packages <- c('tidyverse',
'ggpubr',
'ggstatsplot'
Biocductor_packages <- c('org.Hs.eg.db',
'hgu133a.db',
'CLL',
'hgu95av2.db',
'survminer',
'survival',
'hugene10sttranscriptcluster',
'limma')
if(length(getOption("CRAN"))==0) options(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")
for (pkg in cran_packages){
if (! require(pkg,character.only=T) ) {
install.packages(pkg,ask = F,update = F)
require(pkg,character.only=T)
}
}
# first prepare BioManager on CRAN
if(length(getOption("CRAN"))==0) options(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")
if(!require("BiocManager")) install.packages("BiocManager",update = F,ask = F)
if(length(getOption("BioC_mirror"))==0) options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
# use BiocManager to install
for (pkg in Biocductor_packages){
if (! require(pkg,character.only=T) ) {
BiocManager::install(pkg,ask = F,update = F)
require(pkg,character.only=T)
}
}
1.找到指定ensembl ID與symbol的對(duì)應(yīng)關(guān)系
> ENSG00000000003.13
> ENSG00000000005.5
> ENSG00000000419.11
> ENSG00000000457.12
> ENSG00000000460.15
> ENSG00000000938.11
- 思路：在注釋包中有g(shù)ene_id與symbol寒矿、gene_id與ensembl_id的對(duì)應(yīng)關(guān)系。
#將以上基因id保存在a.txt若债，存放于工作目錄下符相。
rm(list=ls())
options(stringsAsFactors = F)
read.table('R 10/R1')
a=read.table('R 10/R1')
library(org.Hs.eg.db)
g2s=toTable(org.Hs.egSYMBOL) #gene ID symble
g2e=toTable(org.Hs.egENSEMBL)
library(stringr)
unlist(lapply(a$V1,function(x){
strsplit(as.character(x),'[.]')[[1]][1]
}))
a$ensembl_id=unlist(lapply(a$V1,function(x){
strsplit(as.character(x),'[.]')[[1]][1]
}))
tmp=merge(a,g2e,by='ensembl_id')
tmp=merge(tmp,g2s,by='gene_id')
答案1
2.根據(jù)探針名找對(duì)應(yīng)symbol ID
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
- 思路：找到注釋包中探針與symbol的對(duì)應(yīng)關(guān)系然后取子集
rm(list=ls())
options(stringsAsFactors = F)
a=read.table('e2')
library(hgu133a.db)
ids=toTable(hgu133aSYMBOL)
head(ids)
colnames(a)='probe_id'
tmp=merge(ids,a,by='probe_id')
# 第二種方法
match(a$probe_id,ids$probe_id)
tmp2=ids[match(a$probe_id,ids$probe_id),]
- 附：判斷得到的兩組結(jié)果是否一致
# 法一：
identical(tmp1,tmp2) #返回邏輯值
# 法二：
dplyr::setdiff(tmp1,tmp2) #返回兩組的差別【沒(méi)差就返回空】
3.根據(jù)symbol找基因表達(dá)量并作圖
- 找到R包CLL內(nèi)置的數(shù)據(jù)集的表達(dá)矩陣?yán)锩娴腡P53基因的表達(dá)量，并且繪制在 progres.-stable分組的boxplot圖,并通過(guò) ggpubr 進(jìn)行美化蠢琳。
+探針的三大內(nèi)容：表達(dá)矩陣assay/exprs跃闹、探針信息featureData弛房、樣本信息phenoData
symbol
rm(list=ls())
options(stringsAsFactors = F)
suppressPackageStartupMessages(library(CLL))
data(sCLLex)
sCLLex
exprSet=exprs(sCLLex)
library(hgu95av2.db) #這個(gè)包找到目的基因
ids=toTable(hgu95av2SYMBOL)
exprSet <- exprs(sCLLex) #探針的表達(dá)量
pdata <- pData(sCLLex) #sampleID與disease的對(duì)應(yīng)關(guān)系
p2s <- toTable(hgu95av2SYMBOL) #探針與symbol的對(duì)應(yīng)關(guān)系
p3 <- filter(p2s,symbol=='TP53')
# boxplot [find TP53 has 3 probe IDs]
probe_tp53 <- c("1939_at","1974_s_at","31618_at")
i = 3 #可換1，2
boxplot(exprSet[probe_tp53[i],] ~ pdata$Disease)
image.png
4.BRCA1基因表達(dá)量
- 找到BRCA1基因在TCGA數(shù)據(jù)庫(kù)的乳腺癌數(shù)據(jù)集(Breast Invasive Carcinoma (TCGA, PanCancer Atlas))的表達(dá)情況
- 提示：使用http://www.cbioportal.org/index.do 定位數(shù)據(jù)集：http://www.cbioportal.org/datasets
-
該基因有四個(gè)亞型，用ggbetweenstats作圖比較一下洒闸。
image.png
rm(list=ls())
options(stringsAsFactors = F)
a=read.table('plot.txt',sep = '\t',fill=T,header = T)
colnames(a)=c('id','subtype','expression','mut')
dat=a
library(ggstatsplot)
ggbetweenstats(data = dat,x=subtype,y=expression)
image.png
5 .TP53 生存分析
- 找到TP53基因在TCGA數(shù)據(jù)庫(kù)的乳腺癌數(shù)據(jù)集的表達(dá)量分組看其是否影響生存
- 提示使用：http://www.oncolnc.org/
-
思路：生存分析，TP53表達(dá)量分為高低兩組做圖比較
官網(wǎng)生存分析
rm(list=ls())
options(stringsAsFactors = F)
a=read.table('BRCA_7157_50_50.csv',sep = ',',fill=T,header = T)
dat=a
library(ggstatsplot)
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_taga.png")
TP53
- 和官網(wǎng)給的分析相差不大0.139的0.14
- 也可以進(jìn)一步花復(fù)雜的圖
# complex survplot
ggsurvplot(sfit,palette = c("orange", "navyblue"),
risk.table = T, pval = T,
conf.int = T, xlab = "Time of datys",
ggtheme = theme_light(),
ncensor.plot = T)
image.png
ggbetweenstats(data = dat,
x = Group,
y = Expression)
image.png
b=read.table('BRCA.txt',sep = '\t',fill=T,header = T) #txt手蝎，用'\t'趁餐；excel用','。
colnames(b)=c('Patient','subtype','expression','mut')
head(b)
b$Patient=substring(b$Patient,1,12)
tmp=merge(a,b,by='Patient')
table(tmp$subtype)
dat=tmp[tmp$subtype=='BRCA_Basal',]
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)
image.png
6.從表達(dá)矩陣中提取指定基因畫熱圖
- 下載數(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
- 提示：根據(jù)基因名拿到探針I(yè)D菇绵，縮小表達(dá)矩陣?yán)L制熱圖肄渗，沒(méi)有檢查到的基因直接忽略即可。
rm(list=ls())
options(stringsAsFactors = F)
#下載和表達(dá)矩陣
library(GEOquery)
GSE17125 <- "GSE17215"
f=GSE17125
#這個(gè)包需要注意兩個(gè)配置咬最，一般來(lái)說(shuō)自動(dòng)化的配置是足夠的翎嫡。
#Setting options('download.file.method.GEOquery'='auto')
#Setting options('GEOquery.inmemory.gpl'=FLASE)
if(!file.exists(f)){
gset <- getGEO('GSE17215', destdir = ".", getGPL = F, AnnotGPL = F)
save(gset, file = f)
}
load('GSE17215') #載入數(shù)據(jù)
class(gset)
length(gset)
class(gset[[1]])
# 因?yàn)檫@個(gè)GEO數(shù)據(jù)集只有一個(gè)GPL平臺(tái)，所以下載到的是一個(gè)含有一個(gè)元素的list
a=gset[[1]]
dat=exprs(a)
dim(dat)
#改基因
library(hgu133a.db)
ids=toTable(hgu133aSYMBOL)
head(ids)
dat=dat[ids$probe_id,]
#dat=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)
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]]
# 過(guò)濾不再基因名的
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.相關(guān)性計(jì)算和熱圖
下載數(shù)據(jù)集GSE24673的表達(dá)矩陣計(jì)算樣本的相關(guān)性并且繪制熱圖永乌，需要標(biāo)記上樣本分組信息
相關(guān)性分析：
# 7題
rm(list=ls())
options(stringsAsFactors = F)
f='GSE24673_eSet.Rdata'
#下載和表達(dá)矩陣
library(GEOquery)
#這個(gè)包需要注意兩個(gè)配置惑申，一般來(lái)說(shuō)自動(dòng)化的配置是足夠的。
#Setting options('download.file.method.GEOquery'='auto')
#Setting options('GEOquery.inmemory.gpl'=FLASE)
if(!file.exists(f)){
gset <- getGEO('GSE24673', destdir = ".", getGPL = F, AnnotGPL = F)
save(gset, file = f)
}
load('GSE24673_eSet.Rdata') #載入數(shù)據(jù)
class(gset)
length(gset)
class(gset[[1]])
# 因?yàn)檫@個(gè)GEO數(shù)據(jù)集只有一個(gè)GPL平臺(tái)翅雏，所以下載到的是一個(gè)含有一個(gè)元素的list
a=gset[[1]]
dat=exprs(a)
dim(dat)expr[1:4,1:4]
pdata <- pData(geo[[1]])
# 自己根據(jù)pdata第八列做一個(gè)分組信息矩陣
grp <- c('rbc','rbc','rbc',
'rbn','rbn','rbn',
'rbc','rbc','rbc',
'normal','normal')
grp_df <- data.frame(group=grp)
rownames(grp_df) <- colnames(expr)
new_cor <- cor(expr)
pheatmap::pheatmap(new_cor, annotation_col = grp_df)
8.找到芯片平臺(tái)對(duì)應(yīng)的注釋包
- 找到 GPL6244 platform of Affymetrix Human Gene 1.0 ST Array 對(duì)應(yīng)的R的bioconductor注釋包圈驼，并且安裝它！
- https://mp.weixin.qq.com/s/sVSsV2fWZOQmNd72Vb3SmQ
- http://www.reibang.com/p/40b560755cdf
# 8 題
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) #注釋包都加上.db
options()$repos
options()$BioC_mirror
9.找到指定探針和對(duì)應(yīng)的基因
- 下載數(shù)據(jù)集GSE42872的表達(dá)矩陣望几，并且分別挑選出所有樣本的(平均表達(dá)量/sd/mad/)最大的探針绩脆，并且找到它們對(duì)應(yīng)的基因。
# 9題
rm(list=ls())
options(stringsAsFactors = F)
f='GSE42872_eSet.Rdata'
#下載和表達(dá)矩陣
library(GEOquery)
#這個(gè)包需要注意兩個(gè)配置，一般來(lái)說(shuō)自動(dòng)化的配置是足夠的靴迫。
#Setting options('download.file.method.GEOquery'='auto')
#Setting options('GEOquery.inmemory.gpl'=FLASE)
if(!file.exists(f)){
gset <- getGEO('GSE42872', destdir = ".", getGPL = F, AnnotGPL = F)
save(gset, file = f)
}
load('GSE42872_eSet.Rdata') #載入數(shù)據(jù)
class(gset)
length(gset)
class(gset[[1]])
# 因?yàn)檫@個(gè)GEO數(shù)據(jù)集只有一個(gè)GPL平臺(tái)惕味，所以下載到的是一個(gè)含有一個(gè)元素的list
a=gset[[1]]
dat=exprs(a)
dim(dat)
pd=pData(a)
# 平均表達(dá)量/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]
> # 平均表達(dá)量/sd/mad 最大探針
> sort(apply(dat, 1, mean),decreasing = T)[1]
7978905
14.53288
> sort(apply(dat, 1, sd),decreasing = T)[1]
8133876
3.166429
> sort(apply(dat, 1, mad),decreasing = T)[1]
8133876
4.268561
10.limma 差異分析
- 下載數(shù)據(jù)集GSE42872的表達(dá)矩陣，并且根據(jù)分組使用limma做差異分析玉锌，得到差異結(jié)果矩陣
- 接第九題名挥，得到表型信息，然后用limma做差異分析
# 10題
rm(list=ls())
options(stringsAsFactors = F)
f='GSE42872_eSet.Rdata'
#下載和表達(dá)矩陣
library(GEOquery)
#這個(gè)包需要注意兩個(gè)配置芬沉，一般來(lái)說(shuō)自動(dòng)化的配置是足夠的躺同。
#Setting options('download.file.method.GEOquery'='auto')
#Setting options('GEOquery.inmemory.gpl'=FLASE)
if(!file.exists(f)){
gset <- getGEO('GSE42872', destdir = ".", getGPL = F, AnnotGPL = F)
save(gset, file = f)
}
load('GSE42872_eSet.Rdata') #載入數(shù)據(jù)
class(gset)
length(gset)
class(gset[[1]])
# 因?yàn)檫@個(gè)GEO數(shù)據(jù)集只有一個(gè)GPL平臺(tái)，所以下載到的是一個(gè)含有一個(gè)元素的list
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
exprSet[1:4,1:4] # limma接受log后的表達(dá)矩陣丸逸，這一步需要檢查一下
#DEO 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
## 這個(gè)矩陣聲明蹋艺，我們要把progres.組和stable進(jìn)行差異分析
## step1
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_noted,results.csv",quote=F)
head(nrDEG)
