Day4-小龙

R语言
emmm，会用，今天就贴一点代码吧
这是一篇生信技能树的学徒作业，我试着整理了一下，并跑出来啦~这样就学会这个套路，并且学会使用代码。
具体步骤如下：
1.准备相关包和环境
2.下载数据
3.分组、挑选样本、注释探针
4.检查数据质量
5.差异分析与可视化
6.进行富集分析

1.准备相关包和环境+下载数据

# 惯例情况环境，设置镜像 
rm(list = ls()) 
options()$repos options("repos"=c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")) 
options()$BioC_mirror options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/") 
options("repos"=c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
 
# 调用需要用到的包 (GEOmirror是jimmy大神的包，巨好用)
library(GEOquery) 
library(GEOmirror)
 
# 一句代码直接搞定GSE数据集下载 
setwd('F:/R/geo')
eSet<- geoChina('GSE76275') 
str(eSet)
 
#刚刚查看了eSet是一个list,我们直接取第一个，并取出表达矩阵exprSet，画个boxplot看一下数据怎么样 eSet<- eSet[[1]]   
exprSet <- exprs(eSet) 
dim(exprSet) 
head(exprSet[,1:4]) 
boxplot(exprSet,las=2)

boxplot

样本间差异不是很大，可以不用做log

3.分组，挑选探针，注释

# 接下来获取临床信息，对样本进行分组 
pdata<-pData(eSet) 
characteristics <- "triple-negative status:ch1" 
group_list <- pdata[, characteristics] 
table(group_list)
 
not_TN_expr = exprSet[ , grep( "not TN", group_list )] 
TN_expr = exprSet[ , !(colnames(exprSet) %in% colnames(not_TN_expr))]
 
newAssayData = cbind( not_TN_expr, TN_expr)
group_list = c(rep( 'control', ncol( not_TN_expr ) ), rep( 'TNBC',  ncol( TN_expr ) ) ) 
dim( newAssayData ) 
newAssayData[ 1:5, 1:5 ] 
table( group_list )
library(AnnoProbe) 
gpl='GPL570' 
probe2gene=idmap(gpl) 
head(probe2gene)
 
# 去除没有注释的探针 
exprSet =newAssayData[rownames(newAssayData) %in% probe2gene[,1],] 
probe2gene=probe2gene[match(rownames(exprSet), probe2gene[,1]),] 
dim(exprSet) 
dim(probe2gene) 
tail(sort(table(probe2gene[,2])), n = 12L)
 
# 相同的基因只保留最大的那个
 {  MAX = by(exprSet, probe2gene[,2],            
function(x) rownames(x)[ which.max(rowMeans(x))])  
MAX = as.character(MAX)  
exprSet = exprSet[rownames(exprSet) %in% MAX,]  
rownames( exprSet ) = probe2gene[ match( rownames( exprSet ), probe2gene[ , 1 ] ), 2 ] } dim(exprSet) 
exprSet[1:5,1:5] 
AssayData = exprSet 
save(AssayData, group_list, file='final_AssayData.Rdata')

今天就先分享到这，有需要的再私聊分享后续代码哦