文章所有内容均来自生信技能树“生信马拉松-数据挖掘班”授课内容个人整理,如需转载请注明出处。
基因表达芯片、转录组、单细胞、表观遗传、突变……
输入数据是数值型matrix
/data.frame
颜色的变化表示数值的大小
输入数据是一个连续型vector
和一个有重复值的离散型vector
拓展内容:箱线图的介绍
拓展内容:FC与logFC Foldchange(FC):处理组平均值/对照组平均值 logFoldchange(logFC): Foldchange取log2
注意取logFC的分组,切不可做反
1.芯片数据差异分析的起点是一个取过log的matrix,如果拿到的是未log得矩阵,需要自行log
2.P.Value值越小/-log10(P.Value)越大,越有信心认为差异显著
根据这些主成分对样本进行聚类,代表样本的点在坐标轴上距离越远,说明样本差异越大
关注点: 1.同一分组是否分成一簇(组内重复性好) 2.中心点之间是否有距离(组间差别大)
基因表达芯片的原理:探针的表达量代表基因的表达量
基因表达芯片的原理图
options("repos"="https://mirrors.ustc.edu.cn/CRAN/") #设置CRAN镜像
if(!require("BiocManager")) install.packages("BiocManager",update = F,ask = F)
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/") #设置Bioconductor镜像
#CRAN中需要安装的包
cran_packages <- c('tidyr',
'tibble',
'dplyr',
'stringr',
'ggplot2',
'ggpubr',
'factoextra',
'FactoMineR',
'devtools',
'cowplot',
'patchwork',
'basetheme',
'paletteer',
'AnnoProbe',
'ggthemes',
'VennDiagram',
'tinyarray')
#Bioconductor中需要安装的包
Biocductor_packages <- c('GEOquery',
'hgu133plus2.db',
'ggnewscale',
"limma",
"impute",
"GSEABase",
"GSVA",
"clusterProfiler",
"org.Hs.eg.db",
"preprocessCore",
"enrichplot")
for (pkg in cran_packages){
if (! require(pkg,character.only=T,quietly = T) ) {
install.packages(pkg,ask = F,update = F)
require(pkg,character.only=T)
}
} #character.only=T是为了精准的识别包名的字符串
for (pkg in Biocductor_packages){
if (! require(pkg,character.only=T,quietly = T) ) {
BiocManager::install(pkg,ask = F,update = F) #不询问不更新,减少麻烦,确保代码顺利运行
require(pkg,character.only=T)
}
}
#前面的所有提示和报错都先不要管。主要看这里
for (pkg in c(Biocductor_packages,cran_packages)){
require(pkg,character.only=T)
}
#重复运行,如果没有问题控制台则不再出warning信息
#没有任何提示就是成功了,如果有warning xx包不存在,用library检查一下。
#library报错,就单独安装。
###前期准备工作###
rm(list = ls()) #清空环境变量
#打破下载时间的限制,改前60秒,改后10w秒
options(timeout = 100000)
options(scipen = 20) #不要以科学计数法表示 scipen的数值代表阈值,默认1从10万开始,每加1阈值扩大10倍
#传统下载方式
library(GEOquery)
eSet = getGEO("GSE7305", destdir = '.', getGPL = F)
#网速太慢,下不下来怎么办
#1.从网页上下载/发链接让别人帮忙下,放在工作目录里
#2.试试geoChina,只能下载2019年前的表达芯片数据
library(AnnoProbe)
eSet = geoChina("GSE7305") #代替第“eSet = getGEO("GSE7305", destdir = '.', getGPL = F)”
#研究一下这个eSet
class(eSet) #从GEO下载的series为1个List
length(eSet)
eSet = eSet[[1]] #去除list外壳
class(eSet) #由Biobase包里面的“ExpressionSet”对象
简单的对象可以直接用@/$提取子集,复杂对象需要看帮助文档利用函数提取 ★何时用@/$:直接在环境data.frame中点最前面的三角符号查看
exp <- exprs(eSet) #Biobase中特定提取子集的函数
dim(exp) #看行、列数量 若出现异常“0”,GEO确认原始信息
range(exp) #看数据范围决定是否需要log(有特别大的,超过50的),是否有负值,异常值
# exp = log2(exp+1) #需要log才log,不需要就在这句代码前加#号 “+1”防止负数和“0”值出现
boxplot(exp,las = 2) #看是否有异常样本
拓展内容 1.根据箱线图区分正常和异常样本
★★★★★★处理异常样本的方法★★★★★ ①直接删除异常样本 ②
exp =limma :: normalizeBetweenArrays(exp)
2.表达矩阵里的负值情况 ①取过log,有负值——正常 ②没取过log,有负值——错误数据 ③有一半负值——做了标准化 后两种情况一般弃用,非要用的话需要处理原始数据(不推荐新手操作) 附:不同格式原始数据的处理方法链接
pd <- pData(eSet)
p = identical(rownames(pd),colnames(exp));p #鉴定是否一致
if(!p) {
s = intersect(rownames(pd),colnames(exp)) #取交集找出相同部分
exp = exp[,s] #提取exp的列名
pd = pd[s,] #提取pd的行名
}
###如果只有两个分组不需要此段###
k = pd$source_name_ch1 %in% c("Ctrl in adherent culture",
"DPN treatment in adherent culture")
#%in%中“”内的内容必须与pd中对应列完全一致
pd = pd[k,]
exp = exp[,k]
identical(rownames(pd),colnames(exp)) #再次确认临床数据的行基因数据的列名是否一致
dim(exp)
dim(pd)
table(pd$source_name_ch1) #查看提取后的分组名称
gpl_number <- eSet@annotation;gpl_number
save(pd,exp,gpl_number,file = "step1output.Rdata")
以上内容均引用自生信技能树
原创声明:本文系作者授权腾讯云开发者社区发表,未经许可,不得转载。
如有侵权,请联系 cloudcommunity@tencent.com 删除。
原创声明:本文系作者授权腾讯云开发者社区发表,未经许可,不得转载。
如有侵权,请联系 cloudcommunity@tencent.com 删除。