数据分析与数据挖掘研究之一 (上)

简介: 之前做过一些数据分析与数据挖掘相关的工作,最近抽空将之前做的内容简单整理一下,方便查看,主要使用R语言和PERL脚本语言,使用TCGA和ICGC数据库中的临床数据,做类似的分析可以参考一下,如果想查看详细内容与数据可以通过本人的Gitee及Github仓库下载,链接于篇尾附上。

一、标题:Effect of HSP90AB1 on the local immune response of hepatocellular carcinoma and it realtionship to prognosis(HSP90β对肝癌局部免疫的影响及对肝癌患者预后的影响)



二、部分代码及结果展示:


1、整理TCGA数据库肝细胞癌临床数据的部分PERL脚本


use strict;
#use warnings;
use XML::Simple;
opendir(RD, ".") or die $!;
my @dirs=readdir(RD);
closedir(RD);
open(WF,">clinical.xls") or die $!;
print WF "Id\tfutime\tfustat\tAge\tGender\tGrade\tStage\tT\tM\tN\n";
foreach my $dir(@dirs){
  #print $dir . "\n";
  next if($dir eq '.');
  next if($dir eq '..');
  #print $dir . "\n";
  if(-d $dir){
    opendir(RD,"$dir") or die $!;
    while(my $xmlfile=readdir(RD)){
      if($xmlfile=~/\.xml$/){
        #print "$dir\\$xmlfile\n";
        my $userxs = XML::Simple->new(KeyAttr => "name");
        my $userxml="";
        if(-f "$dir/$xmlfile"){
          $userxml = $userxs->XMLin("$dir/$xmlfile");
        }else{
          $userxml = $userxs->XMLin("$dir\$xmlfile");
        }
        # print output
        #open(WF,">dumper.txt") or die $!;
        #print WF Dumper($userxml);
        #close(WF);
        my $disease_code=$userxml->{'admin:admin'}{'admin:disease_code'}{'content'};   #get disease code
        my $disease_code_lc=lc($disease_code);
        my $patient_key=$disease_code_lc . ':patient';                                #ucec:patient
        my $follow_key=$disease_code_lc . ':follow_ups';
        my $patient_barcode=$userxml->{$patient_key}{'shared:bcr_patient_barcode'}{'content'};  #TCGA-AX-A1CJ
        my $gender=$userxml->{$patient_key}{'shared:gender'}{'content'};      #male/female
        my $age=$userxml->{$patient_key}{'clin_shared:age_at_initial_pathologic_diagnosis'}{'content'};
        my $race=$userxml->{$patient_key}{'clin_shared:race_list'}{'clin_shared:race'}{'content'};  #white/black
        my $grade=$userxml->{$patient_key}{'shared:neoplasm_histologic_grade'}{'content'};  #G1/G2/G3
        my $clinical_stage=$userxml->{$patient_key}{'shared_stage:stage_event'}{'shared_stage:clinical_stage'}{'content'};  #stage I
        my $clinical_T=$userxml->{$patient_key}{'shared_stage:stage_event'}{'shared_stage:tnm_categories'}{'shared_stage:clinical_categories'}{'shared_stage:clinical_T'}{'content'};
        my $clinical_M=$userxml->{$patient_key}{'shared_stage:stage_event'}{'shared_stage:tnm_categories'}{'shared_stage:clinical_categories'}{'shared_stage:clinical_M'}{'content'};
        my $clinical_N=$userxml->{$patient_key}{'shared_stage:stage_event'}{'shared_stage:tnm_categories'}{'shared_stage:clinical_categories'}{'shared_stage:clinical_N'}{'content'};
        my $pathologic_stage=$userxml->{$patient_key}{'shared_stage:stage_event'}{'shared_stage:pathologic_stage'}{'content'};  #stage I
        my $pathologic_T=$userxml->{$patient_key}{'shared_stage:stage_event'}{'shared_stage:tnm_categories'}{'shared_stage:pathologic_categories'}{'shared_stage:pathologic_T'}{'content'};
        my $pathologic_M=$userxml->{$patient_key}{'shared_stage:stage_event'}{'shared_stage:tnm_categories'}{'shared_stage:pathologic_categories'}{'shared_stage:pathologic_M'}{'content'};
        my $pathologic_N=$userxml->{$patient_key}{'shared_stage:stage_event'}{'shared_stage:tnm_categories'}{'shared_stage:pathologic_categories'}{'shared_stage:pathologic_N'}{'content'};
        $gender=(defined $gender)?$gender:"unknow";
        $age=(defined $age)?$age:"unknow";
        $race=(defined $race)?$race:"unknow";
        $grade=(defined $grade)?$grade:"unknow";
        $clinical_stage=(defined $clinical_stage)?$clinical_stage:"unknow";
        $clinical_T=(defined $clinical_T)?$clinical_T:"unknow";
        $clinical_M=(defined $clinical_M)?$clinical_M:"unknow";
        $clinical_N=(defined $clinical_N)?$clinical_N:"unknow";
        $pathologic_stage=(defined $pathologic_stage)?$pathologic_stage:"unknow";
        $pathologic_T=(defined $pathologic_T)?$pathologic_T:"unknow";
        $pathologic_M=(defined $pathologic_M)?$pathologic_M:"unknow";
        $pathologic_N=(defined $pathologic_N)?$pathologic_N:"unknow";
        my $survivalTime="";
        my $vital_status=$userxml->{$patient_key}{'clin_shared:vital_status'}{'content'};
        my $followup=$userxml->{$patient_key}{'clin_shared:days_to_last_followup'}{'content'};
        my $death=$userxml->{$patient_key}{'clin_shared:days_to_death'}{'content'};
        if($vital_status eq 'Alive'){
          $survivalTime="$followup\t0";
        }
        else{
          $survivalTime="$death\t1";
        }
        for my $i(keys %{$userxml->{$patient_key}{$follow_key}}){
          eval{
              $followup=$userxml->{$patient_key}{$follow_key}{$i}{'clin_shared:days_to_last_followup'}{'content'};
              $vital_status=$userxml->{$patient_key}{$follow_key}{$i}{'clin_shared:vital_status'}{'content'};
              $death=$userxml->{$patient_key}{$follow_key}{$i}{'clin_shared:days_to_death'}{'content'};
          };
          if($@){
              for my $j(0..5){                       #假设最多有6次随访
                  my $followup_for=$userxml->{$patient_key}{$follow_key}{$i}[$j]{'clin_shared:days_to_last_followup'}{'content'};
                  my $vital_status_for=$userxml->{$patient_key}{$follow_key}{$i}[$j]{'clin_shared:vital_status'}{'content'};
                  my $death_for=$userxml->{$patient_key}{$follow_key}{$i}[$j]{'clin_shared:days_to_death'}{'content'};
                  if( ($followup_for =~ /\d+/) || ($death_for  =~ /\d+/) ){
                          $followup=$followup_for;
                          $vital_status=$vital_status_for;
                          $death=$death_for;
                          my @survivalArr=split(/\t/,$survivalTime);
                          if($vital_status eq 'Alive'){
                            if($followup>$survivalArr[0]){
                              $survivalTime="$followup\t0";
                            }
                          }
                          else{
                            if($death>$survivalArr[0]){
                              $survivalTime="$death\t1";
                            }
                          }
                  }
              }
          }
          my @survivalArr=split(/\t/,$survivalTime);
          if($vital_status eq 'Alive'){
            if($followup>$survivalArr[0]){
              $survivalTime="$followup\t0";
            }
          }
          else{
            if($death>$survivalArr[0]){
              $survivalTime="$death\t1";
            }
          }
        }
        print WF "$patient_barcode\t$survivalTime\t$age\t$gender\t$grade\t$pathologic_stage\t$pathologic_T\t$pathologic_M\t$pathologic_N\n";
      }
    }
    close(RD);
  }
}
close(WF);


2、使用R语言分析正常组与肿瘤组中HSP90AB1的表达情况


#if (!requireNamespace("BiocManager", quietly = TRUE))
#    install.packages("BiocManager")
#BiocManager::install("limma")
#install.packages("ggplot2")
#install.packages("ggpubr")
#引用包
library(limma)
library(ggplot2)
library(ggpubr)
expFile="symbol.txt"     #表达输入文件
gene="VCAN"              #基因的名称
setwd("C:\\Users\\lexb4\\Desktop\\geneImmune\\07.diff")     #设置工作目录
#读取基因表达文件,并对数据进行处理
rt=read.table(expFile, header=T, sep="\t", check.names=F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)), nrow=nrow(exp), dimnames=dimnames)
data=avereps(data)
data=t(data[gene,,drop=F])
#正常和肿瘤数目
group=sapply(strsplit(rownames(data),"\\-"), "[", 4)
group=sapply(strsplit(group,""), "[", 1)
group=gsub("2", "1", group)
conNum=length(group[group==1])       #正常组样品数目
treatNum=length(group[group==0])     #肿瘤组样品数目
Type=c(rep(1,conNum), rep(2,treatNum))
#差异分析
exp=cbind(data, Type)
exp=as.data.frame(exp)
colnames(exp)=c("gene", "Type")
exp$Type=ifelse(exp$Type==1, "Normal", "Tumor")
exp$gene=log2(exp$gene+1)
#设置比较组
group=levels(factor(exp$Type))
exp$Type=factor(exp$Type, levels=group)
comp=combn(group,2)
my_comparisons=list()
for(i in 1:ncol(comp)){my_comparisons[[i]]<-comp[,i]}
#绘制boxplot
boxplot=ggboxplot(exp, x="Type", y="gene", color="Type",
              xlab="",
              ylab=paste0(gene, " expression"),
              legend.title="Type",
              palette = c("blue","red"),
              add = "jitter")+ 
  stat_compare_means(comparisons=my_comparisons,symnum.args=list(cutpoints = c(0, 0.001, 0.01, 0.05, 1), symbols = c("***", "**", "*", "ns")),label = "p.signif")
#输出图片
pdf(file=paste0(gene,".diff.pdf"), width=5, height=4.5)
print(boxplot)
dev.off()


image.png


3、使用R语言分析不同类型免疫细胞在肝细胞癌中的表达水平及相关关系


#install.packages("corrplot")
library(corrplot)                   #引用包
immFile="CIBERSORT-Results.txt"     #免疫细胞浸润的结果文件
pFilter=0.05                        #免疫细胞浸润结果过滤条件
setwd("C:\\Users\\Administrator\\Desktop\\geneimmune\\10immunePlot")    #设置工作目录
#读取免疫细胞浸润的结果文件,并对数据进行整理
immune=read.table(immFile, header=T, sep="\t", check.names=F, row.names=1)
immune=immune[immune[,"P-value"]<pFilter,]
immune=as.matrix(immune[,1:(ncol(immune)-3)])
data=t(immune)
#绘制柱状图
col=rainbow(nrow(data), s=0.7, v=0.7)
pdf(file="barplot.pdf", width=22, height=10)
par(las=1,mar=c(8,5,4,16),mgp=c(3,0.1,0),cex.axis=1.5)
a1=barplot(data,col=col,yaxt="n",ylab="Relative Percent",xaxt="n",cex.lab=1.8)
a2=axis(2,tick=F,labels=F)
axis(2,a2,paste0(a2*100,"%"))
axis(1,a1,labels=F)
par(srt=60,xpd=T);text(a1,-0.02,colnames(data),adj=1,cex=0.6);par(srt=0)
ytick2=cumsum(data[,ncol(data)]);ytick1=c(0,ytick2[-length(ytick2)])
legend(par('usr')[2]*0.98,par('usr')[4],legend=rownames(data),col=col,pch=15,bty="n",cex=1.3)
dev.off()
#删除正常样品
group=sapply(strsplit(colnames(data),"\\-"), "[", 4)
group=sapply(strsplit(group,""), "[", 1)
group=gsub("2", "1", group)
data=data[,group==0,drop=F]
#绘制免疫细胞相关性的图形
pdf(file="corrplot.pdf", width=13, height=13)
par(oma=c(0.5,1,1,1.2))
immune=immune[,colMeans(immune)>0]
M=cor(immune)
corrplot(M,
         method = "color",
         order = "hclust",
         tl.col="black",
         addCoef.col = "black",
         number.cex = 0.8,
         col=colorRampPalette(c("blue", "white", "red"))(50)
         )
dev.off()


image.png

image.png

4、使用R语言分析正常组及肝癌组中不同免疫细胞浸润水平


#install.packages("pheatmap")
#install.packages("vioplot")
#引用包
library(vioplot)
library(pheatmap)
input="CIBERSORT-Results.txt"      #免疫细胞浸润文件
pFilter=0.05                       #免疫细胞浸润结果过滤条件
setwd("C:\\Users\\Administrator\\Desktop\\生信文章\\geneimmune\\11heatmap\\vioplot-high")      #设置工作目录
#读取免疫结果文件,并对数据进行整理
immune=read.table("CIBERSORT-Results.txt", header=T, sep="\t", check.names=F, row.names=1)
immune=immune[immune[,"P-value"]<pFilter,,drop=F]
immune=as.matrix(immune[,1:(ncol(immune)-3)])
data=t(immune)
#正常和肿瘤数目
group=sapply(strsplit(colnames(data),"\\-"), "[", 4)
group=sapply(strsplit(group,""), "[", 1)
group=gsub("2", "1", group)
conNum=length(group[group==1])       #正常组样品数目
treatNum=length(group[group==0])     #肿瘤组样品数目
#定义热图的注释文件
Type=c(rep("Normal",conNum), rep("Tumor",treatNum))
names(Type)=colnames(data)
Type=as.data.frame(Type)
#绘制热图
pdf(file="heatmap.pdf", width=12, height=6)
pheatmap(data, 
         annotation=Type, 
         color = colorRampPalette(c(rep("green",1), rep("black",1), rep("red",3)))(100),
         cluster_cols =F,
         show_colnames=F,
         fontsize = 8,
         fontsize_row=7,
         fontsize_col=5)
dev.off()
#绘制小提琴图
data=t(data)
outTab=data.frame()
pdf(file="vioplot.pdf", width=13, height=8)
par(las=1, mar=c(10,6,3,3))
x=c(1:ncol(data))
y=c(1:ncol(data))
xMax=ncol(data)*3-2
plot(x,y,
     xlim=c(0,xMax),ylim=c(min(data),max(data)+0.02),
     main="", xlab="", ylab="Fraction",
     pch=21,
     col="white",
     xaxt="n")
#对每个免疫细胞循环,绘制小提琴图,正常样品用绿色表示,肿瘤样品用红色表示
for(i in 1:ncol(data)){
  if(sd(data[1:conNum,i])==0){
      data[1,i]=0.00001
  }
  if(sd(data[(conNum+1):(conNum+treatNum),i])==0){
      data[(conNum+1),i]=0.00001
  }
  conData=data[1:conNum,i]
  treatData=data[(conNum+1):(conNum+treatNum),i]
  vioplot(conData,at=3*(i-1),lty=1,add = T,col = 'green')
  vioplot(treatData,at=3*(i-1)+1,lty=1,add = T,col = 'red')
  wilcoxTest=wilcox.test(conData, treatData)
  p=wilcoxTest$p.value
  if(p<pFilter){
      cellPvalue=cbind(Cell=colnames(data)[i], pvalue=p)
    outTab=rbind(outTab, cellPvalue)
  }
  mx=max(c(conData,treatData))
  lines(c(x=3*(i-1)+0.2,x=3*(i-1)+0.8),c(mx,mx))
  text(x=3*(i-1)+0.5, y=mx+0.02, labels=ifelse(p<0.001, paste0("p<0.001"), paste0("p=",sprintf("%.03f",p))), cex = 0.8)
}
legend("topright", 
       c("Normal", "Tumor"),
       lwd=5,bty="n",cex=1.2,
       col=c("green","red"))
text(seq(1,xMax,3),-0.05,xpd = NA,labels=colnames(data),cex = 1,srt = 45,pos=2)
dev.off()
#输出免疫细胞和p值表格文件
write.table(outTab,file="diff.result.txt",sep="\t",row.names=F,quote=F)


image.png

image.png

5、不同拷贝子数目的HSP90β对中性粒细胞和CD8阳性T细胞在肝癌局部浸润的影响


image.png


6、HSP90β基因的表达水平、拷贝子水平及甲基化水平与不同淋巴细胞数量之间的关系

image.png

image.png

image.png

相关文章
|
3月前
|
文字识别 算法 数据挖掘
视觉智能开放平台产品使用合集之对于统计研究和数据分析,有哪些比较好的工具推荐
视觉智能开放平台是指提供一系列基于视觉识别技术的API和服务的平台,这些服务通常包括图像识别、人脸识别、物体检测、文字识别、场景理解等。企业或开发者可以通过调用这些API,快速将视觉智能功能集成到自己的应用或服务中,而无需从零开始研发相关算法和技术。以下是一些常见的视觉智能开放平台产品及其应用场景的概览。
|
1月前
|
存储 数据可视化 数据挖掘
大数据环境下的房地产数据分析与预测研究的设计与实现
本文介绍了一个基于Python大数据环境下的昆明房地产市场分析与预测系统,通过数据采集、清洗、分析、机器学习建模和数据可视化技术,为房地产行业提供决策支持和市场洞察,探讨了模型的可行性、功能需求、数据库设计及实现过程,并展望了未来研究方向。
大数据环境下的房地产数据分析与预测研究的设计与实现
|
1月前
|
供应链 数据可视化 数据挖掘
【2023年第十一届泰迪杯数据挖掘挑战赛】B题:产品订单的数据分析与需求预测 建模及python代码详解 问题一
本文详细介绍了第十一届泰迪杯数据挖掘挑战赛B题的解决方案,涵盖了对产品订单数据的深入分析、多种因素对需求量影响的探讨,并建立了数学模型进行未来需求量的预测,同时提供了Python代码实现和结果可视化的方法。
70 3
【2023年第十一届泰迪杯数据挖掘挑战赛】B题:产品订单的数据分析与需求预测 建模及python代码详解 问题一
|
1月前
|
机器学习/深度学习 数据采集 数据挖掘
【2023年第十一届泰迪杯数据挖掘挑战赛】B题:产品订单的数据分析与需求预测 建模及python代码详解 问题二
本文提供了第十一届泰迪杯数据挖掘挑战赛B题问题二的详细解题步骤,包括时间序列预测模型的建立、多元输入时间预测问题的分析、时间序列预测的建模步骤、改进模型的方法,以及使用Python进行SARIMA模型拟合和预测的具体实现过程。
46 1
|
1月前
|
供应链 算法 数据挖掘
【2023年第十一届泰迪杯数据挖掘挑战赛】B题:产品订单的数据分析与需求预测 23页论文及实现代码
本文介绍了2023年第十一届泰迪杯数据挖掘挑战赛B题的解决方案,深入分析了产品订单数据,并使用Arimax和Var模型进行了需求预测,旨在为企业供应链管理提供科学依据,论文共23页并包含实现代码。
50 0
【2023年第十一届泰迪杯数据挖掘挑战赛】B题:产品订单的数据分析与需求预测 23页论文及实现代码
|
1月前
|
机器学习/深度学习 安全 算法
【2023年第十一届泰迪杯数据挖掘挑战赛】A题:新冠疫情防控数据的分析 32页和40页论文及实现代码
本文总结了2023年第十一届泰迪杯数据挖掘挑战赛A题的新冠疫情防控数据分析,提供了32页和40页的论文以及实现代码,涉及密接者追踪、疫苗接种影响分析、重点场所管控以及疫情趋势研判等多个方面,运用了机器学习算法和SEIR传染病模型等方法。
33 0
【2023年第十一届泰迪杯数据挖掘挑战赛】A题:新冠疫情防控数据的分析 32页和40页论文及实现代码
|
1月前
|
机器学习/深度学习 安全 算法
【2023年第十一届泰迪杯数据挖掘挑战赛】A题:新冠疫情防控数据的分析 建模方案及python代码详解
本文介绍了2023年第十一届泰迪杯数据挖掘挑战赛A题的解题思路和Python代码实现,涵盖了新冠疫情防控数据的分析、建模方案以及数据治理的具体工作。
41 0
【2023年第十一届泰迪杯数据挖掘挑战赛】A题:新冠疫情防控数据的分析 建模方案及python代码详解
|
3月前
|
数据采集 数据可视化 数据挖掘
数据挖掘实战:使用Python进行数据分析与可视化
在大数据时代,Python因其强大库支持和易学性成为数据挖掘的首选语言。本文通过一个电商销售数据案例,演示如何使用Python进行数据预处理(如处理缺失值)、分析(如销售额时间趋势)和可视化(如商品类别销售条形图),揭示数据背后的模式。安装`pandas`, `numpy`, `matplotlib`, `seaborn`后,可以按照提供的代码步骤,从读取CSV到数据探索,体验Python在数据分析中的威力。这只是数据科学的入门,更多高级技术等待发掘。【6月更文挑战第14天】
184 11
|
3月前
|
数据采集 机器学习/深度学习 数据可视化
数据挖掘实战:Python在金融数据分析中的应用案例
Python在金融数据分析中扮演关键角色,用于预测市场趋势和风险管理。本文通过案例展示了使用Python库(如pandas、numpy、matplotlib等)进行数据获取、清洗、分析和建立预测模型,例如计算苹果公司(AAPL)股票的简单移动平均线,以展示基本流程。此示例为更复杂的金融建模奠定了基础。【6月更文挑战第13天】
870 3
|
4月前
|
算法 搜索推荐 数据挖掘
数据挖掘实战 —— 抖音用户浏览行为数据分析与挖掘(续)
数据挖掘实战 —— 抖音用户浏览行为数据分析与挖掘(续)
203 1

热门文章

最新文章