注册网站域名,网站备案拍照好麻烦,做网站技巧,做seo网站优化价格上节#xff0c;我们使用结核病基因数据#xff0c;做了一个数据预处理的实操案例。例子中结核类型#xff0c;包括结核#xff0c;潜隐进展#xff0c;对照和潜隐#xff0c;四个类别。本节延续上个数据#xff0c;进行了差异分析。 差异分析 计算差异指标step12
加载… 上节我们使用结核病基因数据做了一个数据预处理的实操案例。例子中结核类型包括结核潜隐进展对照和潜隐四个类别。本节延续上个数据进行了差异分析。 差异分析 计算差异指标step12
加载数据
load(dataset_TB_LTBI_step8.Rdata)
构建差异比较矩阵
#样本列表
group_listgroup_data_TB_LTBI$group_more #构建分组
designmodel.matrix(~0factor(group_list))
colnames(design)levels(factor(group_list))#head(dataset_TB_LTBI)row.names(design)colnames(dataset_TB_LTBI)
design #得到分组矩阵0代表不是1代表是#str(design)library(limma)
##差异比较矩阵
contrast_matrixmakeContrasts(paste0(c(LTBI,TB),collapse -),levels design)
计算差异基因指标
#step:lmFit
fitlmFit(dataset_TB_LTBI,design)
fit2contrasts.fit(fit,contrast_matrix)
#step:eBayes
fit3eBayes(fit2)#step3:topTable
tempoutputtopTable(fit3,coef 1,nInf)
DEG_Mna.omit(tempoutput) #得到差异分析矩阵重点看logFC和P值
head(DEG_M) #查看数据logFC AveExpr t P.Value adj.P.Val B
ASPHD2 -1.452777 8.415563 -12.38370 5.885193e-22 5.868863e-18 39.30255
C1QC -3.978887 5.971935 -12.34993 6.954041e-22 5.868863e-18 39.14037
GBP1P1 -4.075057 5.607978 -12.24397 1.174622e-21 6.608814e-18 38.63087
GBP6 -3.225604 4.393248 -11.93968 5.320543e-21 1.692866e-17 37.16200
SDC3 -2.374911 7.388880 -11.92896 5.612049e-21 1.692866e-17 37.11012
LHFPL2 -1.705514 8.411180 -11.91494 6.017652e-21 1.692866e-17 37.04225#绘制前40个基因在不同样本之间的热图
library(pheatmap)
#绘制前40个基因在不同样本之间的热图
f40_genehead(rownames(DEG_M),40)
f40_subset_matrixdataset_TB_LTBI[f40_gene,]
head(f40_subset_matrix)
f40_subset_matrixxt(scale(t(f40_subset_matrix))) #数据标准化。。。数据标准化和归一化的区别平移和压缩
pheatmap(f40_subset_matrixx) #出图
差异分析 结果过滤筛选step13 res DEG[,c(logFC,P.Value,adj.P.Val)]colnames(res)-c(logFC,PValue,padj)colnames(res)
library(dplyr)
FC_filter 0.585
P_filter0.05
all_diff res %% filter(abs(logFC)FC_filter) %% filter(padjP_filter)
res$id rownames(res)
resselect(res,id,everything())
#write.table(res,all_diff.txt,sep\t,quoteF)up_diffres %% filter(logFCFC_filter) %% filter(padjP_filter)
up_diff$id rownames(up_diff)
up_diffselect(up_diff,id,everything())
#write.table(up_diff,up_diff.txt,sep\t,quoteF)down_diffres %% filter(logFC -FC_filter ) %% filter(padjP_filter)
down_diff$id rownames(down_diff)
down_diffselect(down_diff,id,everything())
#write.table(down_diff,down_diff.txt,sep\t,quoteF)group_data_clean -function(data){# colnames(data)[c(9,10,11)] c(logFC,PValue,padj)data[which(data$padj %in% NA),sig] - no diffdata[which(data$logFC FC_filter data$padj 0.05),sig] - updata[which(data$logFC -FC_filter data$padj 0.05),sig] - downdata[which(abs(data$logFC) FC_filter | data$padj 0.05),sig] - no diffcat( 上调,nrow(data[data$sig %in% up, ]))cat( 下调,nrow(data[data$sig %in% down, ]))cat( no fiff,nrow(data[data$sig %in% no diff, ]))# filter_data subset(data, data$sig up | data$sig down)# filter_data$Geneid - rownames(filter_data)return(data)
}
limma_clean_res group_data_clean(res)#上调 1381 下调 1432 no fiff 14066rownames(all_diff)dataset_TB_LTBI_DEG dataset_TB_LTBI[rownames(all_diff),]
dim(dataset_TB_LTBI_DEG) #[1] 2813 102#
#数据保存
#
save(DEG,res,all_diff,limma_clean_res,dataset_TB_LTBI_DEG,file DEG_TB_LTBI_step13.Rdata)
#
#数据保存
#差异分析 绘制火山图step14
library(ggplot2)data - limma_clean_res#################
# ggplot2绘制火山图
data$label - c(rownames(data)[1:10],rep(NA,nrow(data) - 10))
#sizeGrWindow(12, 9)
pdf(file差异基因火山图step14.pdf, width 9, height 6)
ggplot(data,aes(logFC,-log10(PValue),color sig)) xlab(log2FC) geom_point(size 0.6) scale_color_manual(valuesc(#00AFBB,#999999,#FC4E07)) geom_vline(xintercept c(-1,1), linetype dashed) geom_hline(yintercept -log10(0.05), linetype dashed) theme(title element_text(size 15), text element_text(size 15)) theme_classic() geom_text(aes(label label),size 3, vjust 1,hjust -0.1)dev.off()差异基因分析完毕下面我们可以观察一下这些基因富集在哪些通路之上。
