精华内容
下载资源
问答
  • 文章目录论文地址四GSE数据集差异表达基因(按logFC值排序)并为一list,正序倒序各一list所有差异基因在四GSE数据集中logFC矩阵筛选共同上调基因筛选共同下调基因合并共同上下调基因logFC.tiff 论文地址 四...

    论文地址

    四个GSE数据集差异表达基因(按logFC值排序)并为一个list,正序倒序各一个list

    setwd("./2.RobustRankAggreg_analysis")
    padj=0.05
    logFC=1
    
    files=c("GSE7476_limmaTab.txt","GSE13507_limmaTab.txt","GSE37815_limmaTab.txt","GSE65635_limmaTab.txt")
    upList=list()
    downList=list()
    allFCList=list()
    for(i in 1:length(files)){
        inputFile=files[i]
        rt=read.table(inputFile,header=T,sep = '\t',quote = '') # 注意文件读取
        header=unlist(strsplit(inputFile,"_"))
        downList[[header[1]]]=as.vector(rt[,1])
        upList[[header[1]]]=rev(as.vector(rt[,1]))
        fcCol=rt[,1:2]
        colnames(fcCol)=c("Gene",header[[1]])
        allFCList[[header[1]]]=fcCol
    }
    
    
    

    所有差异基因在四个GSE数据集中logFC矩阵

    mergeLe=function(x,y){
      merge(x,y,by="Gene",all=T)}
    newTab=Reduce(mergeLe,allFCList)
    rownames(newTab)=newTab[,1]
    newTab=newTab[,2:ncol(newTab)]
    newTab[is.na(newTab)]=0
    

    筛选共同上调基因

    library(RobustRankAggreg)
    upMatrix = rankMatrix(upList)
    upAR = aggregateRanks(rmat=upMatrix)
    colnames(upAR)=c("Name","Pvalue")
    upAdj=p.adjust(upAR$Pvalue,method="bonferroni")
    upXls=cbind(upAR,adjPvalue=upAdj)
    upFC=newTab[as.vector(upXls[,1]),]
    upXls=cbind(upXls,logFC=rowMeans(upFC))
    write.table(upXls,file="up.xls",sep="\t",quote=F,row.names=F)
    upSig=upXls[(upXls$adjPvalue<padj & upXls$logFC>logFC),]
    write.table(upSig,file="upSig.xls",sep="\t",quote=F,row.names=F)
    

    筛选共同下调基因

    downMatrix = rankMatrix(downList)
    downAR = aggregateRanks(rmat=downMatrix)
    colnames(downAR)=c("Name","Pvalue")
    downAdj=p.adjust(downAR$Pvalue,method="bonferroni")
    downXls=cbind(downAR,adjPvalue=downAdj)
    downFC=newTab[as.vector(downXls[,1]),]
    downXls=cbind(downXls,logFC=rowMeans(downFC))
    write.table(downXls,file="down.xls",sep="\t",quote=F,row.names=F)
    downSig=downXls[(downXls$adjPvalue<padj & downXls$logFC< -logFC),]
    write.table(downSig,file="downSig.xls",sep="\t",quote=F,row.names=F)
    

    合并共同上下调基因

    allSig = rbind(upSig,downSig)
    colnames(allSig)
    allSig = allSig[,c("Name","logFC")]
    write.table(allSig,file = 'allSign.xls',sep = '\t',quote = F)
    

    logFC.tiff

    hminput=newTab[c(as.vector(upSig[1:20,1]),as.vector(downSig[1:20,1])),]
    library(pheatmap)
    tiff(file="logFC.tiff",width = 15,height = 20,units ="cm",compression="lzw",bg="white",res=400)
    pheatmap(hminput,display_numbers = TRUE,
             fontsize_row=10,
             fontsize_col=12,
             color = colorRampPalette(c("green", "white", "red"))(50),
             cluster_cols = FALSE,cluster_rows = FALSE, )
    dev.off()
    

    在这里插入图片描述

    本博客内容将同步更新到个人微信公众号生信玩家。欢迎大家关注~~~
    在这里插入图片描述

    展开全文
  • GEO数据挖掘全流程分析

    万次阅读 多人点赞 2020-02-29 19:45:43
    声明:以下学习资料根据“生信技能树”网络系列免费教学材料整理而成,代码来自“生信技能树”校长jimmy的github。...首先便是要广泛阅读,在读文献时,提炼脉络,读懂文献使用了哪个或哪些GSE数据集,对数据做了...

    声明:以下学习资料根据“生信技能树”网络系列免费教学材料整理而成,代码来自“生信技能树”校长jimmy的github。GEO数据库挖掘系列知识分享课程,于2016年首发于生信菜鸟团博客配套教学视频在B站,特此声明。

    前言:关于GEO数据

    我们的目标是要从读懂文献到复刻文献实验,再到掌握GEO数据挖掘的能力。首先便是要广泛阅读,在读文献时,提炼脉络,读懂文献使用了哪个或哪些GSE数据集,对数据做了哪些处理。了解清楚后,便可下载相应的数据集,得到表达矩阵,作差异分析,注释等一系列下游分析。
    一篇文章可以有一个或多个GSE数据集,一个GSE里可以有一个或多个GSM样本。多个研究的GSM样本可以根据研究目的整合为一个GDS,每个数据集有着自己对应的芯片平台(GPL),一个GSE里可能有多个平台测出的数据。
    本分析是基于R语言的平台,所以需要一些R语言的基础知识。
    了解GEO,表达芯片与R

    第一部分:GEO芯片数据下载及整理

    GEO官网
    本例以GSE42872数据集为例,学习GEO数据挖掘分析,分析文献
    在这里插入图片描述
    通过阅读文献找到相应的GSE数据集,并在官网可以下相应的数据集信息及背景知识
    https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE42872
    在这里插入图片描述
    这个数据集使用的是GPL6244这个芯片平台,由6个样本组成,前三个为对照组,后三个为处理组。了解了这个数据集的相关背景后,接下来就需要进行数据下载了。数据下载的方式有很多种,这里我们用R包GEOquery来下载,具体下载代码如下:

    rm(list = ls())  
    ## 魔幻操作,一键清空~当前环境中对象全部删除
    options(stringsAsFactors = F)
    #在调用as.data.frame的时,将stringsAsFactors设置为FALSE可以避免character类型自动转化为factor类型
    f='GSE42872_eSet.Rdata'
    #把GSE42872_eSet.Rdata赋值给f,方便后面流程化处理
    ##根据数据集不同修改相应的GSE号
    
    library(GEOquery)
    # 这个包需要注意两个配置,一般来说自动化的配置是足够的。
    #Setting options('download.file.method.GEOquery'='auto')
    #Setting options('GEOquery.inmemory.gpl'=FALSE)
    
    if(!file.exists(f)){
      gset <- getGEO('GSE42872', destdir=".",
                     AnnotGPL = F,     ## 注释文件
                     getGPL = F)       ## 平台文件
      save(gset,file=f)   ## 保存到本地
    }
    ##这是一个函数,利用包将数据集的表达信息下载下来,赋值给了gset,而不下载注释信息和平台信息,病保存到本地,文件名为f。
    
    load('GSE42872_eSet.Rdata') 
     ## 载入数据
    class(gset) 
     #查看数据类型
    length(gset) 
     ##看一下有几个元素
    gset[[1]]
    #取第一个元素
    class(gset[[1]])
     #查看改元素的数据类型
    # 因为这个GEO数据集只有一个GPL平台,所以下载到的是一个含有一个元素的list
    
    a=gset[[1]] 
    ##取出第一个元素赋值给一个对象a
    dat=exprs(a) 
    #a现在是一个对象,取a这个对象通过看说明书知道要用exprs这个函数,该函数得到表达矩阵
    #现在 得到的dat就是一个表达矩阵,只不过基因的ID是探针名
    dim(dat)
    #看一下dat这个矩阵的维度
    dat[1:5,1:5] 
    #查看dat这个矩阵的1至4行和1至4列,逗号前为行,逗号后为列
    #这个表达矩阵是已经log之后的,表达量一般是0-10左右,如果是原始芯片表达的信号值一般是几千到一万,则需要log处理。
    
    boxplot(dat,las=2) 
    #画个图看一下各样本之间有没有批次效应,一般中位数都差不多,las是将横坐标样本信息竖着排列
    
    pd=pData(a) 
    #通过查看说明书知道取对象a里的临床信息用pData
    View(pd)
    ## 查看一下,挑选一些感兴趣的临床表型,这里我们欲得到其分组title信息。
    library(stringr)
    #运行一个字符分割包
    group_list=str_split(pd$title,' ',simplify = T)[,4]
    #抽取title一列,按照空格分割,取第四个元素即Control和Vemurafenib
    table(group_list)
    #看一下两个分组各有几个
    
    

    也可以使用中国镜像下载:

    library(devtools)
    install_github("jmzeng1314/AnnoProbe")
    library(AnnoProbe)
    gset<-geoChina("GSE12345")
    

    GPL6244这个是该数据的芯片平台,首先在网页查一下有没有对应的包,有的话直接下载相应的包,没有用下面的函数直接下载。
    ##网址:http://www.bio-info-trainee.com/1399.html
    或者在博客查找对应关系
    ##可以查到对应的包为 hugene10sttranscriptcluster
    ##因此我们优先使用包来寻找探针和基因名的对应关系,在bioconductor下载这个包。
    如果下载了相应的芯片平台的包,运行以下程序:

    library(hugene10sttranscriptcluster.db)
    ##运行这个包
    ls("package:hugene10sttranscriptcluster.db")
    #对这个包进行探索,看一下有多少元素
    ##找到有SYMBOL的那个元素就是我们需要的对应关系
    ##例如:[34] "hugene10sttranscriptclusterSYMBOL"    
    ids=toTable(hugene10sttranscriptclusterSYMBOL) 
    #toTable这个函数:通过看hgu133plus2.db这个包的说明书知道提取probe_id(探针名)和symbol(基因名)的对应关系的表达矩阵的函数为toTable。并赋值给ids
    ##这时候我们得到19825个探针对应的基因名。
    ###刚才我们的表达矩阵中是33297个基因探针,这就意味着刚才的表达矩阵中可能存在多个探针重复对应一个基因名。这就需要我们对数据进行进一步筛选、处理。
    head(ids) 
    #head为查看前六行
    

    如果没有查到芯片平台对应的R包就需要下载平台注释信息,按以下操作

    if(F){
      library(GEOquery)
      #Download GPL file, put it in the current directory, and load it:
      gpl <- getGEO('GPL6244', destdir=".")
    ##需要修改相应的平台号,把平台信息赋值给gpl
      colnames(Table(gpl))  
      head(Table(gpl)[,c(1,15)]) 
    ## you need to check this , which column do you need
      probe2gene=Table(gpl)[,c(1,15)]
      head(probe2gene)
      library(stringr)  
      save(probe2gene,file='probe2gene.Rdata')
    }
    load(file='probe2gene.Rdata')
    ids=probe2gene 
    
    

    现在我们就通过两种方法得到了GPL探针和基因的对应关系。

    head(ids)
    #为查看前六行
    colnames(ids)=c('probe_id','symbol')  
    #将列名统一改为'probe_id','symbol'方便后续统一操作。
    
    length(unique(ids$symbol)) 
    #[1] 18832个独特的基因探针,意味着本来19825个里面有一部分是重复的
    tail(sort(table(ids$symbol)))
    table(sort(table(ids$symbol)))
    #每个对象出现的个数
    plot(table(sort(table(ids$symbol))))
    #画图观察
    
    ids=ids[ids$symbol != '',]
    ids=ids[ids$probe_id %in%  rownames(dat),]
    ##%in%用于判断是否匹配,然后取匹配的几行,去掉无法匹配的信息。
    
    dat[1:5,1:5]   
    dat=dat[ids$probe_id,] 
    #取表达矩阵中可以与探针名匹配的那些,去掉无法匹配的表达数据,这时只剩下19825个探针及表达信息,其余已被剔除。
    
    ids$median=apply(dat,1,median) 
    #ids新建median这一列,列名为median,同时对dat这个矩阵按行操作,取每一行的中位数,将结果给到median这一列的每一行
    ids=ids[order(ids$symbol,ids$median,decreasing = T),]
    #对ids$symbol按照ids$median中位数从大到小排列的顺序排序
    ##即先按symbol排序,相同的symbol再按照中位数从大到小排列,方便后续保留第一个值。
    ##将对应的行赋值为一个新的ids,这样order()就相当于sort()
    ids=ids[!duplicated(ids$symbol),]
    #将symbol这一列取取出重复项,'!'为否,即取出不重复的项,去除重复的gene ,保留每个基因最大表达量结果.最后得到18832个基因。
    dat=dat[ids$probe_id,] 
    #新的ids取出probe_id这一列,将dat按照取出的这一列中的每一行组成一个新的dat
    rownames(dat)=ids$symbol
    #把ids的symbol这一列中的每一行给dat作为dat的行名
    head(dat)
    ##至此我们就得到了该数据集的表达矩阵,最后将结果保存。
    save(dat,group_list,file = 'step1-output.Rdata')
    write.csv(dat,file="expressionmetrix_GSE.csv")
    

    第二部分:数据分布判断及分析

    检查样本分布

    exprSet=dat
    exprSet['GAPDH',]
    boxplot(exprSet[,1])
    boxplot(exprSet['GAPDH',])
    exprSet['ACTB',]
    library(reshape2)
    ?melt()
    class(exprSet)
    exprSet_L=melt(exprSet,)
    colnames(exprSet_L)=c('symbol','sample','value')
    exprSet_L$group=rep(group_list,each=nrow(exprSet))
    head(exprSet_L)
    

    关于melt函数的用法及介绍可以参考博客
    下面进行绘图

    ### ggplot2 
    library(ggplot2)
    p=ggplot(exprSet_L,
             aes(x=sample,y=value,fill=group))+geom_boxplot()
    print(p)
    pdf('result1.pdf')#width = 5,height = 10)
    p
    dev.off()
    
    p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_violin()
    print(p)
    pdf('result2.pdf')#width = 5,height = 10)
    p
    dev.off()
    p=ggplot(exprSet_L,aes(value,fill=group))+geom_histogram(bins = 200)+facet_wrap(~sample, nrow = 4)
    print(p)
    pdf('result3.pdf')#width = 5,height = 10)
    p
    dev.off()
    p=ggplot(exprSet_L,aes(value,col=group))+geom_density()+facet_wrap(~sample, nrow = 4)
    pdf('result4.pdf')#width = 5,height = 10)
    p
    dev.off()
    p=ggplot(exprSet_L,aes(value,col=group))+geom_density() 
    pdf('result5.pdf')#width = 5,height = 10)
    p
    dev.off()
    p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_boxplot()
    p=p+stat_summary(fun.y="mean",geom="point",shape=23,size=3,fill="red")
    p=p+theme_set(theme_set(theme_bw(base_size=20)))
    p=p+theme(text=element_text(face='bold'),axis.text.x=element_text(angle=30,hjust=1),axis.title=element_blank())
    
    pdf('result6.pdf')#width = 5,height = 10)
    p
    dev.off()
    

    PCA图

    hc=hclust(dist(t(exprSet)))
    plot(hc)
    colnames(exprSet)=paste(group_list,1:78,sep = " ")
    hc=hclust(dist(t(exprSet)))
    p<-plot(hc)
    pdf('result7.pdf')#width = 5,height = 10)
    p
    dev.off()
    library(ggfortify)
    df=as.data.frame(t(exprSet))
    df$group<-group_list
    p<-autoplot(prcomp(df[,1:(ncol(df)-1)]),
             data=df,
             colour="group")
    pdf('result8.pdf')#width = 5,height = 10)
    p
    dev.off()
    

    第三部分:利用limma包差异表达分析

    rm(list = ls())  ## 魔幻操作,一键清空~
    options(stringsAsFactors = F)
    load(file = 'step1-output.Rdata')
    # 每次都要检测数据
    dat[1:4,1:4] 
    table(group_list) #table函数,查看group_list中的分组个数
    #通过为每个数据集绘制箱形图,比较数据集中的数据分布
    boxplot(dat[1,]~group_list) #按照group_list分组画箱线图
    
    bp=function(g){         #定义一个函数g,函数为{}里的内容
      library(ggpubr)
      df=data.frame(gene=g,stage=group_list)
      p <- ggboxplot(df, x = "stage", y = "gene",
                     color = "stage", palette = "jco",
                     add = "jitter")
      #  Add p-value
      p + stat_compare_means()
    }
    bp(dat[1,]) ## 调用上面定义好的函数,避免同样的绘图代码重复多次敲。
    bp(dat[2,])
    bp(dat[3,])
    bp(dat[4,])
    dim(dat)
    

    用limma包做GEO芯片合并后的差异表达分析:

    #source("http://bioconductor.org/biocLite.R")
    #biocLite("limma")
    
    logFoldChange=1
    adjustP=0.05
    
    library(limma)
    setwd("C:\\Users\\lexb4\\Desktop\\geoBatch\\07.diff")
    rt=read.csv("namalized.csv")
    rt=as.matrix(rt)
    ###解决基因名重复的问题
    rownames(rt)=rt[,1]
    exp=rt[,2:ncol(rt)]
    dimnames=list(rownames(exp),colnames(exp))
    rt=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
    ##删除基因名重复
    
    #differential定义样品类型
    modType=c(rep("normal",25),rep("tumor",25),rep("normal",5),rep("tumor",5))
    #定义肿瘤类型
    design <- model.matrix(~0+factor(modType))
    colnames(design) <- c("normal","tumor")
    
    fit <- lmFit(rt,design)
    allDiff[1:5,1:5]
    cont.matrix<-makeContrasts(tumor-normal,levels=design)
    #处理组减去对照组;
    #去了log之后的值做减法就是除法,得到的值就是logFC
    fit2 <- contrasts.fit(fit, cont.matrix)
    fit2 <- eBayes(fit2)
    
    allDiff=topTable(fit2,adjust='fdr',number=200000)
    #校正方法FDR,定义200000就保证了所有的基因的差异情况
    write.csv(alldiff,file="alldiff.csv")
    
    #write table
    diffSig <- allDiff[with(allDiff, (abs(logFC)>logFoldChange & adj.P.Val < adjustP )), ]
    write.table(diffSig,file="diff.xls",sep="\t",quote=F,row.names=F)
    diffUp <- allDiff[with(allDiff, (logFC>logFoldChange & adj.P.Val < adjustP )), ]
    write.table(diffUp,file="up.xls",sep="\t",quote=F,row.names=F)
    diffDown <- allDiff[with(allDiff, (logFC<(-logFoldChange) & adj.P.Val < adjustP )), ]
    write.table(diffDown,file="down.xls",sep="\t",quote=F,row.names=F)
    #提取差异基因表达矩阵
    hmExp=rt[as.vector(diffSig[,1]),]
    diffExp=rbind(id=colnames(hmExp),hmExp)
    write.table(diffExp,file="diffExp.txt",sep="\t",quote=F,col.names=F)
    
    ##画火山图
    pdf(file="vol.pdf")
    xMax=max(-log10(allDiff$adj.P.Val))
    yMax=max(abs(allDiff$logFC))
    ##最大值可以手动输入
    plot(-log10(allDiff$adj.P.Val), allDiff$logFC, xlab="-log10(adj.P.Val)",ylab="logFC",
         main="Volcano", xlim=c(0,xMax),ylim=c(-yMax,yMax),yaxs="i",pch=20, cex=0.8)
    diffSub=subset(allDiff, adj.P.Val<adjustP & logFC>logFoldChange)
    points(-log10(diffSub$adj.P.Val), diffSub$logFC, pch=20, col="red",cex=0.8)
    diffSub=subset(allDiff, adj.P.Val<adjustP & logFC<(-logFoldChange))
    points(-log10(diffSub$adj.P.Val), diffSub$logFC, pch=20, col="green",cex=0.8)
    ##cex为点的大小,可以自己设置
    abline(h=0,lty=2,lwd=3)
    dev.off()
    

    接着来画个热图(并且定义样本类型)

    #install.packages("pheatmap")
    
    setwd("C:\\....pheatmap")      #设置工作目录
    rt=read.table("diffExp.txt",sep="\t",header=T,check.names=F)
    rt=as.matrix(rt)
    ##检查是否有重复的基因名
    rownames(rt)=rt[,1]
    exp=rt[,2:ncol(rt)]
    dimnames=list(rownames(exp),colnames(exp))
    rt=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
    
    #rt=log2(rt+1)
    #rt[rt>15]=15
    
    library(pheatmap)
    ##注释文件
    Geo=c(rep("GSE33335",50),rep("GSE56807",10))
    ##前面10个来源于GSE33335,后面10个来源于GSE56807
    Type=c(rep("N",25),rep("T",25),rep("N",5),rep("T",5))    #修改正常和癌症样品数目
    
    names(Geo)=colnames(rt)
    ann=cbind(Geo,Type)
    ann=as.data.frame(ann)
    
    tiff(file="heatmap.tiff",
           width = 20,            #图片的宽度
           height =25,            #图片的高度
           units ="cm",
           compression="lzw",
           bg="white",
           res=500)  ##图片分辨率
    pheatmap(rt, annotation=ann, 
             color = colorRampPalette(c("green", "black", "red"))(50),  ##地,中,高表达颜色
             cluster_cols =F,
    #         scale="row",   校正,视图形而定
             fontsize_row=3,
             fontsize_col=5)
    dev.off()
    

    第四部分 功能富集分析

    在做富集分析之前,要先把基因的symble转换成基因的entrezID
    首先准备基因symble文件,一列为symble,一列为logFC

    #install.packages("RSQLite")
    #source("https://bioconductor.org/biocLite.R")
    #biocLite("org.Hs.eg.db")
    
    setwd("C:\\Users\\lexb4\\Desktop\\geoBatch\\09.symbol2id")
    
    library("org.Hs.eg.db")
    rt=read.table("symbol.txt",sep="\t",check.names=F,header=T)
    genes=as.vector(rt[,1])
    entrezIDs <- mget(genes, org.Hs.egSYMBOL2EG, ifnotfound=NA)
    entrezIDs <- as.character(entrezIDs)
    out=cbind(rt,entrezID=entrezIDs)
    write.table(out,file="id.txt",sep="\t",quote=F,row.names=F)
    

    运行后就多了一列为基因的entrezID。接下来就可以进行基因富集分析,结果会以柱状图和气泡图的类型呈现。

    #install.packages("colorspace")
    #install.packages("stringi")
    #source("http://bioconductor.org/biocLite.R")
    #biocLite("DOSE")
    #biocLite("clusterProfiler")
    #biocLite("pathview")
    
    setwd("C:\\...GO")
    library("clusterProfiler")
    library("org.Hs.eg.db")
    rt=read.table("id.txt",sep="\t",header=T,check.names=F)
    rt=rt[is.na(rt[,"entrezID"])==F,]
    ##保留ID不是NA的基因列表
    geneFC=rt$logFC
    gene=rt$entrezID
    names(geneFC)=gene
    
    #GO富集分析
    kk <- enrichGO(gene = gene,OrgDb = org.Hs.eg.db, pvalueCutoff =0.05, qvalueCutoff = 0.05)
    write.table(kk,file="GO.txt",sep="\t",quote=F,row.names = F)
    
    #柱状图
    tiff(file="barplot.tiff",width = 20,height = 30,units ="cm",compression="lzw",bg="white",res=300)
    barplot(kk, drop = TRUE, showCategory = 47)
    dev.off()
    
    #点图
    tiff(file="dotplot.tiff",width = 20,height = 30,units ="cm",compression="lzw",bg="white",res=300)
    dotplot(kk,showCategory = 47)
    dev.off()
    

    KEGG富集分析(柱状图,气泡图,通路图)
    通路中有颜色的就是差异基因
    输入文件:ID.txt

    #install.packages("colorspace")
    #install.packages("stringi")
    #source("http://bioconductor.org/biocLite.R")
    #biocLite("DOSE")
    #biocLite("clusterProfiler")
    #biocLite("pathview")
    
    setwd("C:\\.....KEGG")
    library("clusterProfiler")
    rt=read.table("id.txt",sep="\t",header=T,check.names=F)
    rt=rt[is.na(rt[,"entrezID"])==F,]
    
    geneFC=rt$logFC
    gene=rt$entrezID
    names(geneFC)=gene
    
    #kegg富集分析
    kk <- enrichKEGG(gene = gene, organism = "hsa", pvalueCutoff =0.05, qvalueCutoff =0.05)
    ##然后看一下有多少个通路有意义
    ##太少了可以把p值或者q值适当调整
    write.table(kk,file="KEGG.txt",sep="\t",quote=F,row.names = F)
    
    #柱状图
    tiff(file="barplot.tiff",width = 20,height = 20,units ="cm",compression="lzw",bg="white",res=300)
    barplot(kk, drop = TRUE, showCategory = 12)
    dev.off()
    
    #点图
    tiff(file="dotplot.tiff",width = 20,height = 20,units ="cm",compression="lzw",bg="white",res=300)
    dotplot(kk)
    dev.off()
    
    #通路图
    library("pathview")
    keggxls=read.table("KEGG.txt",sep="\t",header=T)
    for(i in keggxls$ID){
      pv.out <- pathview(gene.data = geneFC, pathway.id = i, species = "hsa", out.suffix = "pathview")
    }
    
    
    展开全文
  • 合并shp文件 用到的工具 UnionFeatureCollection 需要的依赖 <dependency> <groupId>org.geotools</groupId> <artifactId>gt-process</artifactId> <version>${geo...

    项目需要对多个shp文件做合并和求交,而geotools常用的几何工具都是针对单个Geometry的,若手动获取shp里的要素集合并遍历每个要素的geometry,可能效率会很低,所以找到了geotools的process工具,测试了一下

    合并两个shp文件

    用到的工具
    UnionFeatureCollection
    需要的依赖

            <dependency>
                <groupId>org.geotools</groupId>
                <artifactId>gt-process</artifactId>
                <version>${geotools.version}</version>
            </dependency>
            <dependency>
                <groupId>org.geotools</groupId>
                <artifactId>gt-process-feature</artifactId>
                <version>${geotools.version}</version>
            </dependency>
    

    代码

        // ShapeFileReaderUtils和ShapeFileWriterUtils是自己定义的,内容见下面地址
        // https://blog.csdn.net/Neuromancerr/article/details/119981046
        public static void mergeTest() throws Exception {
            // 待合并的两个文件地址
            String pathHead = "D:\\mapresources\\test\\";
            String file1 = pathHead + "test1.shp";
            String file2 = pathHead + "test2.shp";
            // 获取要素集
            ContentFeatureCollection features1 = ShapeFileReaderUtils.getFeatureSourceFromShapeFile(file1).getFeatures();
            ContentFeatureCollection features2 = ShapeFileReaderUtils.getFeatureSourceFromShapeFile(file2).getFeatures();
            // 定义输出路径
            String outPut = "D:\\mapresources\\test\\testOut.shp";
            // 初始化
            UnionFeatureCollection collection2 = new UnionFeatureCollection();
            // 执行合并
            SimpleFeatureCollection result = collection2.execute(
                    features1,
                    features2
            );
            // 输出结果生成shp文件
            ShapeFileWriterUtils.buildShp(result, outPut);
        }
    
        public static void main(String[] args) throws Exception {
            mergeTest();
        }
    

    结果成功
    合并前
    合并后

    两个shp求交

    用到的是 IntersectionFeatureCollection
    以下是测试代码

    // ShapeFileReaderUtils和ShapeFileWriterUtils是自己定义的,内容见下面地址
    // https://blog.csdn.net/Neuromancerr/article/details/119981046
    public static void mergeTest() throws Exception {
        // 待合并的两个文件地址
        String pathHead = "D:\\mapresources\\test\\";
        String file1 = pathHead + "test1.shp";
        String file2 = pathHead + "test2.shp";
        // 获取要素集
        ContentFeatureCollection features1 = ShapeFileReaderUtils.getFeatureSourceFromShapeFile(file1).getFeatures();
        ContentFeatureCollection features2 = ShapeFileReaderUtils.getFeatureSourceFromShapeFile(file2).getFeatures();
        // 定义输出路径
        String outPut = "D:\\mapresources\\test\\testOut3.shp";
        // 初始化
        IntersectionFeatureCollection collection = new IntersectionFeatureCollection();
        List<String> attr1 = new ArrayList<>();
        features1.features().next().getAttributes().forEach(attr -> {
            attr1.add(attr.toString());
        });
        List<String> attr2 = new ArrayList<>();
        features2.features().next().getAttributes().forEach(attr -> {
            attr2.add(attr.toString());
        });
        // 执行合并
        SimpleFeatureCollection result = collection.execute(
                features1,
                features2,
                attr1,
                attr2,
                IntersectionFeatureCollection.IntersectionMode.FIRST,
                true,
                true
        );
    
        // 这一步是为了消除多边形自相交的几何错误
        SimpleFeature[] resArr = (SimpleFeature[]) result.toArray();
        DefaultFeatureCollection res2 = new DefaultFeatureCollection();
        for (SimpleFeature simpleFeature : resArr) {
            simpleFeature.setDefaultGeometry(parseGeom((Geometry) simpleFeature.getDefaultGeometry()));
            res2.add(simpleFeature);
        }
    
        // 输出结果生成shp文件
        ShapeFileWriterUtils.buildShp(res2, outPut);
    }
    
    // Geometry转multiPolygon
    private static MultiPolygon parseGeom(Geometry geometry) {
    
        MultiPolygon multiPolygon = null;
        // 防止自相交
        geometry = geometry.intersection(geometry);
    
        if (geometry.getGeometryType().equals("MultiPolygon")) {
            multiPolygon = (MultiPolygon) geometry;
        } else if (geometry.getGeometryType().equals("Polygon")) {
            GeometryFactory gf = new GeometryFactory();
            Polygon[] polys = new Polygon[1];
            polys[0] = (Polygon) geometry;
            multiPolygon = gf.createMultiPolygon(polys);
        }
        return multiPolygon;
    }
    

    但结果并不是期望中两个要素的相交部分,而是IntersectionFeatureCollection.execute()方法传入的第二个要素集的整个复制。这把我搞懵了,经过各种参数调整测试,都没啥用,只能暂时放弃。

    求交的另一种思路

    用的方法比较笨,先用UnionFeatureCollection把需要求交的shp都合成一个巨大的数据集,然后遍历其中的要素,分别对每个geometry求交。我是这么写的

    // 传入一整个要素集
    private void mergeShp(DefaultFeatureCollection features) throws Exception {
        String outPut = "D:\\mapresources\\jiangsu\\output\\test.shp";
        // 输出要素集
        DefaultFeatureCollection featureCollection = new DefaultFeatureCollection();
        // 输入要素集转数组
        SimpleFeature[] featuresArray = new SimpleFeature[features.size()];
        featuresArray = features.toArray(featuresArray);
        for (int i = 0; i < featuresArray.length; i++) {
            // 取出单个要素和他的geometry
            SimpleFeature feature1 = featuresArray[i];
            MultiPolygon polygon1 = parseGeom((Geometry) feature1.getDefaultGeometry());
            boolean intersected = false;
            for (int j = 0; j < i; j++) {
                SimpleFeature feature2 = featuresArray[j];
                MultiPolygon polygon2 = parseGeom((Geometry) feature2.getDefaultGeometry());
                assert polygon1 != null;
                assert polygon2 != null;
                if (polygon1.intersects(polygon2)) {
                    try {
                        intersected = true;
                        SimpleFeature newFeature;
                        // 按照定义的mdate字段,较新的属性覆盖旧的
                        if (Long.parseLong((String) feature1.getAttribute("mdate")) > Long.parseLong(((String) feature2.getAttribute("mdate")))) {
                            newFeature = SimpleFeatureBuilder.copy(feature1);
                        } else {
                            newFeature = SimpleFeatureBuilder.copy(feature2);
                        }
                        newFeature.setDefaultGeometry(polygon1.union(polygon2));
                        featureCollection.add(newFeature);
                        System.out.println("featureCollection add:" + newFeature.getID() + " intersected " + feature1.getID() + "," + feature2.getID());
                    } catch (Exception e) {
                        e.printStackTrace();
                    }
    
                }
    
            }
    
            if (!intersected) {
                featureCollection.add(feature1);
                System.out.println("featureCollection add:" + feature1.getID());
            }
        }
    
        ShapeFileWriterUtils.buildShp(featureCollection, outPut);
    }
    

    问题首先是效率,对两个各有799个要素的shp求交,耗时数分钟。同时若有多个几何同时相交,会重复产生新要素,比如A、B、C相交,会产生A-B和B-C两个要素,B的部分重负了。
    虽然算法上可以并应当优化,但负责python的同事可以接手这部分工作,遂放弃。在这里记录一下。

    展开全文
  • GEO数据库数据下载方法总结

    千次阅读 2021-03-10 10:10:30
    数据集GSE1001为例, 可以直接点击“Series Matrix Files”获取该样本txt格式的表达谱数据,一般我们认为这种处理过的表达谱数据是没有问题的,当然,具体情况具体分析。 打开下载的文件可以看到许多“#”开头的...

    GEO数据下载

    GEO是生信分析经常用到的数据库。经常需要从中获取表达矩阵,平台信息,meta信息等,本博文总结了几种下载GEO数据的方法,各有优劣,实际应用过程中自行选择适合自己的。

    方法一:直接从浏览器中下载,手动

    以数据集GSE1001为例,
    在这里插入图片描述
    可以直接点击“Series Matrix Files”获取该样本txt格式的表达谱数据,一般我们认为这种处理过的表达谱数据是没有问题的,当然,具体情况具体分析。
    打开下载的文件可以看到许多“#”开头的行,这些是注释信息,一般关注这些注释信息中的“data processing”,这行中可以看到数据是如何归一化和标准化的,以及是否已经经过log转化等。
    之后我们需要下载平台数据以注释表达谱中的探针,点击“GPL85”,点击下图中的红框,
    在这里插入图片描述
    之后可以用R进行表达谱和平台数据的合并。

    注意!到这里我们获得的只是初步的表达数据,还没有经过预处理,需要用R处理多个探针对应一个表达值,无对应symbol,以及合并多个探针对应一个symbol的情况后才可进行后续分析。

    优缺点:

    优点:该方法的优点是下载较快,数据也比较完整。
    缺点:无法直接得到metadata,需要手动通过GEO2R获取,或在表达谱的注释信息中寻找。

    方法二:代码下载,自动化

    常用的R包有GEOquery,直接copy了“生信技能树”的代码,以下代码直接将表达谱封装成expression set对象,这有一个非常显著的优点即可以直接导入limma进行差异分析。

    downGSE <- function(studyID = "GSE1009", destdir = ".") {
    
        library(GEOquery)
        eSet <- getGEO(studyID, destdir = destdir, getGPL = F)
    
        exprSet = exprs(eSet[[1]])
        pdata = pData(eSet[[1]])#表型信息,里面包含了需要的分组信息
    
        write.csv(exprSet, paste0(studyID, "_exprSet.csv"))
        write.csv(pdata, paste0(studyID, "_metadata.csv"))
        return(eSet)
    
    }
    
    #导入limma包进行分析
    library(limma)
    design=model.matrix(~factor(sCLLex$Disease))
    fit=lmFit(eSet,design)
    fit=eBayes(fit)
    options(digits = 4)
    topTable(fit,coef=2,adjust='BH')
    

    关于原始芯片数据的处理可参考以下文章link of original microarray

    如果自己手头上已经有了表达谱数据和分组信息,也可以手动构造expression set对象:

    metadata <- data.frame(labelDescription=c('SampleID', 'Disease'),
                           row.names=c('SampleID', 'Disease'))
        phenoData <- new("AnnotatedDataFrame",data=meta,varMetadata=metadata)
        myExpressionSet <- ExpressionSet(assayData=exprMatrix,
                                         phenoData=phenoData,
                                         annotation="hgu95av2")
        > myExpressionSet
        ExpressionSet (storageMode: lockedEnvironment)
        assayData: 12625 features, 22 samples 
          element names: exprs 
        protocolData: none
        phenoData
          sampleNames: CLL11.CEL CLL12.CEL ... CLL9.CEL (22 total)
          varLabels: SampleID Disease
          varMetadata: labelDescription
        featureData: none`在这里插入代码片`
        experimentData: use 'experimentData(object)'
        Annotation: hgu95av2 
        > 
    

    如果只是想获得某一GSE的metadata,可以通过GEOmetadb包获取:

    library(GEOmetadb)
    if(!file.exists('GEOmetadb.sqlite')) getSQLiteFile()
    ## 取决于网速哦
    file.info('/path/GEOmetadb.sqlite')
    con <- dbConnect(SQLite(),'/path/GEOmetadb.sqlite'))
    #dbListTables(con2)
    #dbListFields(con2,'gse')
    GeoList = read.table("diabetes.GEO.list")
    query = paste("select + from gsm where series_id in
                (  ' ", gsub(", ", " ', ' ", paste(Geolist[,1], collapse=",")," '", seq=" ")
    query
    tmp = dbGetQuery(con2, query)
    write.csv(tmp, "diabetes.GEO.meta.csv")
    

    优缺点:

    优点:自动化
    缺点:取决于网速,下载慢

    展开全文
  • 数据挖掘—GEO,TCGA,Oncomine联合数据(一)挖掘概述 数据挖掘—GEO,TCGA,Oncomine联合(二)GEO在线工具的应用 前言 要做分析那肯定要下载数据,这下载数据的过程大家肯定都会,但是下载完的数据真的能直接就...
  • # 以数据集GSE89657为例,芯片平台是GPL6244。 1.下载表达谱数据 # GEO网站手动下载表达谱数据,解压,去注释 gunzip GSE89657_series_matrix.txt.gz cat GSE89657_series_matrix.txt|grep -v "\!" >GSE...
  • 欢迎关注”生信修炼手册”!GEO数据库中的数据是公开的,很的科研工作者会下载其中的数据自己去分析,其中差异表达分析是最常见的分析策略之一,为了方便大家更好的挖掘GEO中的数据,官网提供...
  • GeoJson数据合并

    千次阅读 2020-06-19 20:06:55
    本文主要是基于geojson-merge,实现多个geojson文件合并为一个geojson文件,以便实现基于该文件进行数据分析展示 geojson合并概述 当前在 datav的geoatlas中,可以下载单个地市或区县的数据,例如福建省下面每个...
  • GEO数据下载 GEOquery包安装 注意:需要R 4.0及以上 if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("GEOquery") .cel数据下载 getGEOSuppFiles...
  • ML主要数据集收藏

    2015-03-21 06:22:36
    ML主要数据集收藏 2014-05-20 09:19 462人...AWS(亚马逊网络服务)公共数据集,提供了一集中的资料库,可以无缝集成到基于AWS的云应用程序的公共数据集。生物测定数据,在 虚拟筛选,生物测定数据,对化学信
  • 文章目录论文地址GSE7476GSE7476数据下载到表达矩阵GSE7476数据下载getGEO包下载的探针注释文件不全,需要在GEO网站下载筛选探针分位数标准化预处理分组差异表达表达矩阵分组矩阵差异表达矩阵按照logFC排序保存差异...
  • 数据集汇总 一、免费大数据存储库的网站 1、深度学习数据集收集网站 http://deeplearning.net/datasets/** 收集大量的各深度学习相关的数据集,但并不是所有开源的数据集都能在上面找到相关信息。 2、Tiny ...
  • 数据集汇总

    2019-09-29 17:16:00
    一、免费大数据存储库的...收集大量的各深度学习相关的数据集,但并不是所有开源的数据集都能在上面找到相关信息。 2、Tiny Images Datasethttp://horatio.cs.nyu.edu/mit/tiny/data/index.html 包含8000万的32x3...
  • GEO是当今最大、最全的公共基因数据资源库,包括基因的表达、突变、修饰等信息,涵盖几乎所有的疾病,且单个实验检测样品数目较,是我们分析、学习的很好资源。 实验设计 原始文章对14溃疡性结肠炎病人 ...
  • geohash和geohash聚合

    2020-05-22 15:07:24
    GeoHash将二维的经纬度转换成字符串,比如下图展示了北京9区 域的GeoHash字符串,分别是WX4ER,WX4G2、WX4G3等等,每 一字符串代表了某一矩形区域。也就是说,这矩形区域内所有的 点(经纬度坐标)都共享相同...
  • 目前移动终端相当普及,一部手机就是一移动位置源,很应用都基于LBS功能,附近的某某(餐馆、银行、KTV等等)。一般地,基础信息数据中都会维护了标记位置的经纬度,利用用户提供的经纬度,进行对比,从而获得...
  • AWS(亚马逊网络服务)公共数据集,提供了一集中的资料库,可以无缝集成到基于AWS的云应用程序的公共数据集。生物测定数据,在 虚拟筛选,生物测定数据,对化学信息学,J.由阿曼达Schierz的,有21生物测定数据...
  • 机器学习相关数据集

    2015-08-06 09:38:00
    KDD杯的中心,所有的数据,任务和结果。 UCI机器学习和知识发现研究中使用的大型...AWS(亚马逊网络服务)公共数据集,提供了一集中的资料库,可以无缝集成到基于AWS的云应用程序的公共数据集。 生物测定数...
  • Geohash原理

    2020-06-12 15:46:19
    GeoHash方式建立空间索引,可以提高对空间poi数据进行经纬度检索的效率。 2.认识GeoHash GeoHash将二维的经纬度转换成字符串,比如下图展示了北京9区域的GeoHash字符串,分别是WX4ER,WX4G2...
  • Geohash算法

    千次阅读 2018-11-23 15:54:08
    1. 引言 ...以GeoHash方式建立空间索引,可以提高对空间poi数据进行经纬度检索的效率。 2. 认识Geohash GeoHash将二维的经纬度转换成字符串,比如下图展示了北京9区域的GeoHash字符串,分别是WX...
  • 100款机器学习数据集

    千次阅读 2018-06-18 22:12:08
    Kaggle   ...书籍推荐数据集(goodreads/上万图书/百万评价)【Kaggle】 https://www.kaggle.com/zygmunt/goodbooks-10k   带有预期点数和获胜概率的NFL比赛详情数据集(2009-2016)【Kaggle】 ...
  • 免费数据集整理

    千次阅读 2017-09-13 13:50:06
    Wikipedia:Database :向感兴趣的用户提供所有可用的内容的免费...这个数据保存在亚马逊s3bucket中,请求者可能花费一些钱来访问它。 Common crawl :建立并维护一个开放的网络,向所有人开放。 EDRM File Formats D
  • 实体对齐(Entity Alignment)相关论文与数据集整理

    千次阅读 多人点赞 2021-03-21 21:46:13
    实体对齐(Entity Alignment)、知识图谱融合方法总结整理 年份 模型 ...实体对齐数据集整理 名称 — DBpe-dia(DBP) LinkedGeoData(LGD) Geonames(GEO) YAGO ...
  • 我想同时利用GeoLayout和Maps of countries layouts两工具。我已经下载了航空公司的sample.gexf文件,并更改经纬度名为lat 和lng。下面是我试过的步骤: Try1 1)打开.gexf文件。 2)应用Maps of countries ...
  • bitmap是一种伪数据类型,是基于String实现的。因为redis的key和value本身就支持二进制的存储方式,所以bitmaps只是一独特的扩展。因为是面向字节操作,所以他的最大长度就是512M,最适合设置成2^32不同字节。 ...
  • 不同平台的数据,同一平台的不同时期的数据,同一样品不同试剂的数据,以及同一样品不同时间的数据等等都会产生一种 batch effect。 校正批次效应的目的是减少batch之间的差异,从而更好的获取不同生物学状态...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 2,370
精华内容 948
关键字:

geo多个数据集合并