临沂品牌网站建设公司,温州企业自助建站系统,东方城乡与住房建设部网站,怎么制作链接视频教程前面五节#xff0c;我们使用阿尔兹海默症数据做了一个数据预处理案例#xff0c;包括如下内容#xff1a;
GEO生信数据挖掘#xff08;一#xff09;数据集下载和初步观察
GEO生信数据挖掘#xff08;二#xff09;下载基因芯片平台文件及注释
GEO生信数据挖掘…
前面五节我们使用阿尔兹海默症数据做了一个数据预处理案例包括如下内容
GEO生信数据挖掘一数据集下载和初步观察
GEO生信数据挖掘二下载基因芯片平台文件及注释
GEO生信数据挖掘三芯片探针ID与基因名映射处理
GEO生信数据挖掘四数据清洗离群值处理、低表达基因、归一化、log2处理
GEO生信数据挖掘五提取临床信息构建分组分组数据可视化绘制层次聚类图绘制PCA图 本节目录
结核病基因表达数据GSE107994观察
临床形状数据预处理
基因表达数据预处理
绘图观察数据 结核病基因表达数据GSE107994观察
由于在数据分析过程你拿的数据样式可能会有不同本节我们以结核病基因表达数据GSE107994为例做一个实践案例。该数据集的临床形状数据和基因表达数据是单独分开的读取和处理都需自己改动代码。 先来看看基因表达数据这个探针注释工作已经完成了不需要处理。 再看看临床形状数据需要手工删除前面的注释把后半部分规整的数据保留下来。 临床形状数据预处理
# 手工删除前面的注释读文件转置 pdata - t(read.delim(GSE107994_series_matrix_clean.txt, header TRUE, sep \t))
# 手工删除前面的注释读文件转置
pdata - t(read.delim(GSE107994_series_matrix_clean.txt, header TRUE, sep \t))
pdata -pdata[-1,]
pdata_info pdata[,c(1,7)]
colnames(pdata_info) c(geo_accession,type)#观察样本类型的取值都有哪些结核潜隐进展对照和潜隐
unique(pdata_info[,2])
#Leicester_Active_TB Longitudnal_Leicester_LTBI_Progressor Leicester_Control #Leicester_LTBI group_data as.data.frame(pdata_info)
处理前
处理后 增加不同的类型标签根据需要选取实验组和对照组 # 使用grepl函数判断字符串是否包含Control并进行相应的修改
group_data$group_easy - ifelse(grepl(Control, group_data$type ), Control, TB)# 使用grepl函数判断字符串是否包含特定内容然后进行相应的修改
group_data$group_easy - ifelse(grepl(Control, group_data$type), Control,ifelse(grepl(LTBI, group_data$type), LTBI,TB))
# 使用grepl函数判断字符串是否包含特定内容然后进行相应的修改
group_data$group_more - ifelse(grepl(Control, group_data$type), Control,ifelse(grepl(LTBI_Progressor, group_data$type), LTBI_Progressor,ifelse(grepl(LTBI, group_data$type), LTBI,TB)))#尝试把进展组排除出去save(group_data,file group_data.Rdata)例如 我们可以进行 TB结核 和LTBI(潜隐结核)实验对照分析。 基因表达数据预处理 读取数据集
install.packages(openxlsx)
library(openxlsx)# 读基因表达矩阵第一列为基因名ID
gse_info- as.data.frame(read.xlsx(GSE107994_Raw.xlsx, sheet 1))
colnames(gse_info)后续运行代码过程中发现基因名称中有全数字的情况这里做删除操作。
library(dplyr)
dim(gse_info)
#基因里面有数字
gse_info - gse_info[!grepl(^\\d$, gse_info$ID), ] #有效#基因名全为空
gse_info gse_info[gse_info$ID ! ,] #无剔除
dim(gse_info) #[1] 58023 176#负值处理
gse_info[gse_info 0] - 0.0001#重复值检查
table(duplicated(gse_info$ID))
分组数据条件筛选TB结核 和LTBI(潜隐结核)
#
#
#type分组数据条件筛选step3
##预处理之前先筛选出TB组和LTBI组 的数据
unique(group_data[,group_more]) #TB LTBI_Progressor Control LTBI #TB LTBI 对照则剔除 LTBI_Progressor Control
geo_accession_TB_LTBI - group_data[group_data$group_more LTBI_Progressor | group_data$group_more Control,geo_accession]
gse_TB_FTBI gse_info[,!(names(gse_info) %in% geo_accession_TB_LTBI)]gse_TB_FTBI低表达过滤(平均值小于1) #
#
#删除 低表达(平均值小于1)基因 step4
#
##新增一列计算平均
gene_avg_expression - rowMeans(gse_TB_FTBI[, -1]) # 计算每个基因的平均表达量排除第一列基因名
#仅去除在所有样本里表达量都为零的基因平均值小于1
gse_TB_FTBI_filtered_genes_1 - gse_TB_FTBI[gene_avg_expression 1, ]
低表达过滤方案二保留样本表达的排名前50%的基因
#
#
#删除 低表达(排名前50%)基因 step5
#
##仅保留在一半以上样本里表达的基因# 计算基因表达矩阵每个基因的平均值
gene_means - rowMeans(gse_TB_FTBI_filtered_genes_1[,-1])# 计算基因平均值的排序百分位数
gene_percentiles - rank(gene_means) / length(gene_means)# 获取阈值
threshold - 0.25 # 删除后25%的阈值
#threshold - 0.5 # 删除后50%的阈值
# 根据阈值筛选低表达基因
gse_TB_FTBI_filtered_genes_2 - gse_TB_FTBI_filtered_genes_1[gene_percentiles threshold, ]# 打印筛选后的基因表达矩阵
dim(gse_TB_FTBI_filtered_genes_2) #[1] 17049 176删除重复基因取平均
#
#
#重复基因取平均值 step6
#
#dim(filtered_genes_2)
table(duplicated(filtered_genes_2$ID))#把重复的Symbol取平均值
averaged_data - aggregate(.~ID , filtered_genes_2, mean, na.action na.pass) ##把重复的Symbol取平均值#把行名命名为SYMBOL
row.names(averaged_data) - averaged_data$ID
dim(averaged_data)#去掉缺失值
matrix_na na.omit(averaged_data) #删除Symbol列一般是第一列
matrix_final - subset(matrix_na, select -1)
dim(matrix_final) #[1] 22687 175
离群值处理 #
#
#离群值处理 step7
#
##数据离群处理
#处理极端值
#定义向量极端值处理函数
#用于处理异常值将超出一定范围的值替换为中位数以减少异常值对后续分析的影响。
dljdzfunction(x) {DOWNBquantile(x,0.25)-1.5*(quantile(x,0.75)-quantile(x,0.25))UPBquantile(x,0.75)1.5*(quantile(x,0.75)-quantile(x,0.25))x[which(xDOWNB)]quantile(x,0.5)x[which(xUPB)]quantile(x,0.5)return(x)
}#第一列设置为行名
matrix_leavematrix_final_TB_LTBIboxplot(matrix_leave,outlineFALSE, notchT, las2) ##出箱线图
dim(matrix_leave)#处理离群值
matrix_leave_resapply(matrix_leave,2,dljdz)boxplot(matrix_leave_res,outlineFALSE, notchT, las2) ##出箱线图
dim(matrix_leave_res)
log2 处理
#
#
#log2 处理 step8
#
## limma的函数归一化矫正差异 表达矩阵自动log2化#1.归一化不是绝对必要的但是推荐进行归一化。
#有重复的样本中应该不具备生物学意义的外部因素会影响单个样品的表达
#例如中第一批制备的样品会总体上表达高于第二批制备的样品假设所有样品表达值的范围和分布都应当相似
#需要进行归一化来确保整个实验中每个样本的表达分布都相似。
#2.归一化要在log2标准化之前做library(limma) exprSetnormalizeBetweenArrays(matrix_leave_res)boxplot(exprSet,outlineFALSE, notchT, las2) ##出箱线图## 这步把矩阵转换为数据框很重要
class(exprSet) ##注释此时数据的格式是矩阵Matrix
exprSet - as.data.frame(exprSet)#标准化 表达矩阵自动log2化
qx - as.numeric(quantile(exprSet, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rmT))
LogC - (qx[5] 100) ||(qx[6]-qx[1] 50 qx[2] 0) ||(qx[2] 0 qx[2] 1 qx[4] 1 qx[4] 2)#负值全部置为空
#exprSet[exprSet 0] - 0.0001
#去掉缺失值
#exprSet na.omit(exprSet) #15654
#save (exprSet,file waitlog_data_TB_LTBI.Rdata)## 开始判断
if (LogC) { exprSet [which(exprSet 0)] - NaN## 取log2exprSet_clean - log2(exprSet1) #是否加一 加一的话不产生负值###%%print(log2 transform finished)
}else{print(log2 transform not needed)
}boxplot(exprSet_clean,outlineFALSE, notchT, las2) ##出箱线图dataset_TB_LTBI exprSet_clean
绘图观察数据 #
#
#对照组不同颜色画箱线图 step9
#
## 使用grepl函数判断字符串是否包含LTBI并进行颜色标记为了画图
group_data_TB_LTBI$group_color - ifelse(grepl(LTBI, group_data_TB_LTBI$group_more), yellow, blue)#画箱线图查看数据分布
group_list_color group_data_TB_LTBI$group_color
boxplot( data.frame(dataset_TB_LTBI),outlineFALSE,notchT,colgroup_list_color,las2)dev.off()#
#
#绘制层次聚类图 step10
#
#
##检查表达矩阵的样本名称和分租信息的样本名称顺序,是否一致对应
colnames(dataset_TB_LTBI)
group_data_TB_LTBI$geo_accessionexprSet dataset_TB_LTBI
#修改GSM的名字改为分组信息
#colnames(exprSet)paste(group_list,1:ncol(exprSet),sep )#定义nodePar
nodeParlist(lab.cex0.6,pchc(NA,19),cex0.7,colblue)
#聚类
hchclust(dist(t(exprSet))) #t()的意思是转置#绘图
plot(as.dendrogram(hc),nodePar nodePar,horiz TRUE)dev.off()#
#
#绘制PCA散点样本可视化图 step11
#
###PCA图
#install.packages(ggfortify)
library(ggfortify)
dfas.data.frame(t(exprSet)) #转置后就变成了矩阵
dim(df) #查看数据维度
dim(exprSet)df$groupgroup_data_TB_LTBI$group_more #加入样本分组信息
autoplot(prcomp(df[,1:ncol(df)-1]),datadf,colourgroup) #PCA散点图dev.off() 至此我们对两个数据集进行了预处理工作下面我们可以对处理完毕的数据进行差异分析了。
