精华内容
下载资源
问答
  • 差异基因分析

    2013-03-08 11:14:47
    美国斯坦福大学开发的基因芯片差异基因分析EXCEL插件,简单方便
  • 生信(一)——DESeq2差异基因分析 文章目录生信(一)——DESeq2差异基因分析一、差异基因分析原理二、代码实现1、前提:安装DESeq2包2.代码实现三、小结 记录学习过程,共勉。 一、差异基因分析原理 详见 二、...

    生信(一)——DESeq2差异基因分析


    记录学习过程,共勉。

    一、差异基因分析原理

    详见

    二、代码实现

    1、前提:安装DESeq2包

    在这里插入图片描述

    2.代码实现
    setwd("D:\\RData");#设置编码位置
    rt<-read.table("GSE149549_mRNA_Expression_Summary.txt",head=TRUE,row.names = 1)
    head(rt)
    #准备
    rt<-as.matrix(rt) #将数据转换为矩阵格式
    condition<-factor(c(rep("case",2),rep("control",3)),levels = c("case","control")) ## 设置分组信息,建立环境(5个样本,2组处理)
    condition
    coldata<-data.frame(row.names=colnames(rt),condition) 
    #剔除无意义数据(不存在的基因read)
    #nrow(rt)
    i=1
    while(i<=nrow(rt)){
      if(mean(rt[i,])<10)
        rt<-rt[-i,]
      i<-i+1;
    };
    #nrow(rt)
    
    coldata #展示coldata内容
    
    condition  #展示环境
    #构建dds矩阵
    dds<-DESeqDataSetFromMatrix(rt,coldata,design=~condition)
    #标准化,对原始数据进行normalize
    dds<-DESeq(dds)
    dds
    resultsNames(dds)
    #使用deseq2包中的result()函数提取deseq2差异分析结果
    res=results(dds,contrast = c("condition","case","control"))
    head(res)
    summary(res)
    write.csv(res,file = "results.csv")
    #分析结果可视化与导出
    table(res$padj<0.05)
    plotMA(res,ylim=c(-2,2))
    mcols(res,use.names = TRUE)
    #提取差异基因
    res <- res[order(res$padj),]
    resdata <-merge(as.data.frame(res),as.data.frame(counts(dds,normalize=TRUE)),by="row.names",sort=FALSE)
    deseq_res<-data.frame(resdata)
    up_diff_result<-subset(deseq_res,padj < 0.05 & (log2FoldChange > 1)) #提取上调差异表达基因
    down_diff_result<-subset(deseq_res,padj < 0.05 & (log2FoldChange < -1)) #提取下调差异表达基因
    write.csv(up_diff_result,"up_case_VS_control_diff_results.csv") #输出上调基因
    write.csv(down_diff_result,"down_case_VS_control_diff_results.csv") #输出下调基因
    
    #plot(res$log2FoldChange,res$padj)#绘制火山图
    
    

    注:DESeq2的表达数据必须是read counts,不能使用标准化后的FPKM 或TPM
    在这里插入图片描述
    由图可见,表达量(基因的样本平均count)数值越大,离散度越小,差异基因倍数越少

    #提取两种算法标准化后的表达量(VST和rlog,大样本通常使用VST,推荐使用)
    vsd<-vst(dds,blind = FALSE)
    rld<-rlog(dds,blind = FALSE)
    #对照使用 log 2(n+1)的标准化方法
    ntd<-normTransform(dds)
    head(assay(vsd),3)
    head(assay(rld),3)
    head(assay(ntd),3)
    

    在这里插入图片描述

    #导出标准化后的数据(可用于GSEA等下游分析)
    write.csv(as.data.frame(assay(vsd)),file="count_transformation_VST.csv")
    write.csv(as.data.frame(assay(rld)),file="count_transformation_rlog.csv")
    
    #绘制离散度散点图,离群值未做收缩
    plotDispEsts(dds)
    

    在这里插入图片描述

    #对标准化后的数据进行可视化处理
    #安装vsn包
    #library(BiocManager)
    #BiocManager::install("vsn")
    #load package
    #library("vsn")
    #library(ggplot2)
    #library(cowplot)
    
    ntdl<-meanSdPlot(assay(ntd))
    vsdl<-meanSdPlot(assay(vsd))
    rldl<-meanSdPlot(assay(rld))
    
    p1<-ntdl$gg+ggtitle("log2(n+1)")
    p2<-vsdl$gg+ggtitle("VST")
    p3<-rldl$gg+ggtitle("rlog")
    
    plot_grid(p1,p2,p3,ncol = 3)
    

    在这里插入图片描述
    可见使用log2(n+1)标准化方法在低counts区域标准差较大,相反,rlog算法较大;而使用VST标准化后数据的标准差比较稳定;注意,虽然VST方法看起来效果最好,但图中的纵轴为方差的平方根(标准差),需要注意可能为由实验处理产生的真实方差。

    #使用标准化后平均表达量top100的基因绘制热图
    library("pheatmap")
    select<-order(rowMeans(counts(dds,normalized=TRUE)),decreasing = TRUE)[1:100]
    colData(dds)
    df<-as.data.frame(colData(dds)[,c("condition","sizeFactor")])
    mycolors=colorRampPalette(c("orange","white","purple"))(100)
    pheatmap(assay(vsd)[select,],cluster_rows=FALSE,show_rownames=FALSE,border_color=NA,color=mycolors,cluster_cols=FALSE,annotation_col=df)
    

    在这里插入图片描述

    #绘制样本距离聚类热图
    #计算样本距离矩阵
    sampleDists<-dist(t(assay(vsd)))
    #由向量转成矩阵
    sampleDistMatrix<-as.matrix(sampleDists)
    rownames(sampleDistMatrix)<-paste(vsd$condition,vsd$sizeFactor,seq="-")
    colnames(sampleDistMatrix)<-NULL
    
    library("RColorBrewer")
    colors<-colorRampPalette(rev(brewer.pal(8,"YlGnBu")))(255)
    pheatmap(sampleDistMatrix,clustering_distance_rows = sampleDists,clustering_distance_cols = sampleDists,col=colors)
    

    在这里插入图片描述

    三、小结

    刚入门,继续学习,加油!

    展开全文
  • GEO多芯片差异基因分析时候出现以下错误:Error in rt[as.vector(diffSig[, 1]), ] : only 0’s may be mixed with negative subscripts。合并了两个数据,数据数值都小于10,应该是已经取了log的。但是为什么会出现...

    GEO多芯片差异基因分析时候出现以下错误:Error in rt[as.vector(diffSig[, 1]), ] :
    only 0’s may be mixed with negative subscripts。合并了两个数据,数据数值都小于10,应该是已经取了log的。但是为什么会出现这个情况,求大神告知!!!!!!!!

    logFoldChange=1
    > adjustP=0.05
    > 
    > library(limma)
    > setwd("C:\\Users\\acer\\Desktop\\DEGs\\07.diff")
    > rt=read.table("normalize.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)
    > 
    > #differential
    > modType=c(rep("normal",6),rep("Endometrium",6),rep("normal",4),rep("Endometrium",4),rep("normal",9),rep("Endometrium",9))
    > design <- model.matrix(~0+factor(modType))
    > colnames(design) <- c("con","treat")
    > fit <- lmFit(rt,design)
    > cont.matrix<-makeContrasts(treat-con,levels=design)
    > fit2 <- contrasts.fit(fit, cont.matrix)
    > fit2 <- eBayes(fit2)
    > 
    > allDiff=topTable(fit2,adjust='fdr',number=200000)
    > write.table(allDiff,file="limmaTab.xls",sep="\t",quote=F,row.names=F)
    > 
    > #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)
    > 
    > #write expression level of diff gene
    > hmExp=rt[as.vector(diffSig[,1]),]
    **Error in rt[as.vector(diffSig[, 1]), ] : 
      only 0's may be mixed with negative subscripts**
    > diffExp=rbind(id=colnames(hmExp),hmExp)  找不到对象
    > write.table(diffExp,file="diffExp.txt",sep="\t",quote=F,col.names=F)   找不到对象
    > 
    > #volcano
    > 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)
    > abline(h=0,lty=2,lwd=3)
    > dev.off()
    
    展开全文
  • setwd("C:/Users/linqingquan/Desktop/GSE_121787"); Sys.setlocale('LC_ALL','Chinese'); GPL_table = read.table('GPL21185-21174.txt',sep="\t",comment.char="#", stringsAsFactors=F,...
    setwd("C:/Users/linqingquan/Desktop/GSE_121787");
    Sys.setlocale('LC_ALL','Chinese');
    
    GPL_table = read.table('GPL21185-21174.txt',sep="\t",comment.char="#",
                           stringsAsFactors=F,header=T,fill=TRUE,quote="");
    
    GSE_matrix <- read.table('GSE121787_series_matrix.txt',sep="\t",comment.char="!",
                             stringsAsFactors=F,header=T,fill=TRUE);
    ID_Sybmol = GPL_table[,c(1,6)];
    colnames(ID_Sybmol)[2]="Symbol";
    
    Exp=merge(ID_Sybmol,GSE_matrix,by.x="ID",by.y="ID_REF",all=T);
    Exp=Exp[,-1];
    View(Exp)
    
    Exp=Exp[Exp$Symbol != "",];
    Exp=na.omit(Exp);
    
    Exp1=data.frame(Exp[-grep("/",Exp$"Symbol"),]);
    meanfun <- function(x) {
      x1 <- data.frame(unique(x[,1]));
      colnames(x1) <- c("Symbol");
      for (i in 2:ncol(x)) {
        x2 <- data.frame(tapply(x[,i],x[,1],mean));
        x2[,2] <- rownames(x2);
        colnames(x2) <- c(colnames(x)[i],"Symbol");
        x1 <- merge(x1,x2,by.x="Symbol",by.y="Symbol",all=FALSE);
      }
      return(x1);
    }
    
    
    Exp2 <- meanfun(Exp1);
    
    par(cex=0.7);
    n.sample=ncol(Exp2[,-1]);
    if(n.sample>40)par(cex=0.5);
    cols <- rainbow(n.sample*1.2);
    boxplot(Exp2[,-1],col=cols,main="expression value",las=2)
    
    write.table(Exp2,"Exp_Original.txt",row.names=F,quote=F,sep="\t")
    
    row.names(Exp2)=Exp2[,1];
    Exp2=log(Exp2[,-1]);
    
    par(cex=0.7);
    n.sample=ncol(Exp2);
    if(n.sample>40)par(cex=0.5);
    cols <- rainbow(n.sample*1.2);
    boxplot(Exp2,col=cols,main="expression value",las=2);
    
    Symbol=row.names(Exp2);
    Exp_test=cbind(Symbol,Exp2);
    write.table(Exp_test,"Exp.txt",row.names=F,quote=F,sep="\t")
    
    
    top_diff <- read.table('top_diff.txt',sep="\t",comment.char="!", stringsAsFactors=F,header=T,fill=TRUE);
    f2=merge(ID_Sybmol,top_diff,by.x="ID",by.y="ID",all=T)
    f3=merge(f2,Exp_test,by.x = "Symbol",by.y = "Symbol",all=F);
    f4=na.omit(f3)
    f5 <- data.frame(unique(f4[,1]))
    colnames(f5) <- c("ss")
    f6=merge(f4,f5,by.x="Symbol",by.y = "ss",all=F)
    f6 <- f4[!duplicated(f4$Symbol),]
    write.table(f6,"merge_diff_gene.txt",row.names=F,quote=F,sep="\t")

    上述代码参考自:https://www.jianshu.com/p/bc5d885ff4cf

    若有相关疑惑,可点击上述链接进行详细解读。

    若感兴趣 请参见github:https://github.com/Mr-lq7/Bioinformatics_Project2

     

     

     

    展开全文
  • 差异基因分析R代码

    2019-01-06 21:35:09
    资源有用,希望能够对大家有帮助。好资源大家gongxiang
  • 四、DESeq2差异基因分析 获得reads-counts之后,我们就可以开展差异基因分析了。我们以subread中的featureCounts工具得到的counts_id.txt为例,来进行后续的差异基因分析。 目前常见的差异基因分析工具有DESeq2、...

    四、DESeq2差异基因分析

    获得reads-counts之后,我们就可以开展差异基因分析了。我们以subread中的featureCounts工具得到的counts_id.txt为例,来进行后续的差异基因分析。
    目前常见的差异基因分析工具有DESeq2、limma包等等,此处以DESeq2为工具来进行差异基因的筛选。

    First step:获取表达矩阵和分组信息

    进行差异基因分析之前,首先要获取表达矩阵和分组信息。我们的表达矩阵是刚才用featureCounts定量得到的counts_id.txt ,经过格式处理之后,是这样(部分截取):
    表达矩阵

    第一列是基因ID,之后的列都是样本ID
    每一行代表不同的基因在不同样本中的表达量.


    我们选day0和day1做比较,
    为了方便,分组矩阵的制作我在R里面完成,输入如下代码:

    us_count<-read.csv("C:\\Users\\admin\\Desktop\\RNA_seq_Ara_counts_day1_day0.csv",head=T,row.names=1) #输入表达矩阵数据路径
    us_count<-round(us_count,digits=0) #将输入数据取整
    
    #准备
    us_count<-as.matrix(us_count) #将数据转换为矩阵格式
    condition<-factor(c("day0","day0","day0","day0","day1","day1","day1","day1")) ## 设置分组信息,建立环境(8个样本,2组处理)
    coldata<-data.frame(row.names=colnames(us_count),condition)  
    
    coldata #展示coldata内容
    condition  #展示环境

    上面的代码运行之后,我们的分组信息就是这样哒:
    分组信息

    现在,就可以正式使用DESeq2做差异基因分析了,总共其实只有三步:

    • 构建dds矩阵
    • 对dds矩阵进行标准化
    • 提取结果并绘制火山图

    代码如下:

    library(DESeq2) #使用library函数加载DEseq2包
    ##构建dds矩阵 
    dds<-DESeqDataSetFromMatrix(us_count,coldata,design=~condition)
    head(dds) #查看构建好的矩阵
    
    ##进行差异分析
    dds<-DESeq(dds) #对原始的dds进行标准化
    resultsNames(dds)  #查看结果名称
    res<-results(dds) #用results函数提取结果,并赋值给res变量
    summary(res) #查看结果
    plotMA(res,ylim=c(-2,2)) 
    mcols(res,use.names=TRUE)
    
    plot(res$log2FoldChange,res$pvalue) #绘制火山图
    
    #提取差异基因
    res <- res[order(res$padj),]
    resdata <-merge(as.data.frame(res),as.data.frame(counts(dds,normalize=TRUE)),by="row.names",sort=FALSE)
    deseq_res<-data.frame(resdata)
    up_diff_result<-subset(deseq_res,padj < 0.05 & (log2FoldChange > 1)) #提取上调差异表达基因
    down_diff_result<-subset(deseq_res,padj < 0.05 & (log2FoldChange < -1)) #提取下调差异表达基因
    
    write.csv(up_diff_result,"C:\\Users\\admin\\Desktop\\上调_day0_VS_day1_diff_results.csv") #输出上调基因
    write.csv(down_diff_result,"C:\\Users\\admin\\Desktop\\下调_day0_VS_day1_diff_results.csv") #输出下调基因

    至此,差异基因就成功提取了,看看火山图(火山图绘制是这句代码起效果:plot(res$log2FoldChange,res$pvalue) #其实就是对log2foldchange和pvalue作图):
    火山图

    五、 总结

    自己之前处理线虫数据主要是用的RSEM(加bowtie2)来做RNA-seq的序列比对和生物学定量,但是没有接触过salmon和subread。这次使用这两个工具完成了植物的RNA-seq实战,对这些工具有了更深的理解,以后自己如果要开发相关软件,也要多用用,比较工具之间的优劣。

    展开全文
  • DESeq2差异基因分析和批次效应移除

    万次阅读 多人点赞 2018-07-26 15:40:35
    差异基因鉴定 基因表达标准化 不同样品的测序量会有差异,最简单的标准化方式是计算 counts per million (CPM),即原始reads count除以总reads数乘以1,000,000。 这种计算方式的缺点是容易受到极高表达且在...
  • 使用limma包进行差异基因分析时,做最多的是两分类的,例如control组和disease组,但也会碰到按照序列进行的分组。这时,如果逐一使用两两比较求差异基因则略显复杂。其实开发limma包的大神们已经替我们考虑到。我...
  • 基因表达分析时必然会要做差异分析(DE) DE的方法主要有两种: Fold change t-test fold change的意思是样本质检表达量的差异倍数,log2 fold change的意思是取log2,这样可以可以让差异特别大的和差异比较小...
  • R语言进行TCGA配对样本差异基因分析

    千次阅读 2021-03-31 00:37:20
    关于DEseq2配对差异分析,很少有帖子涉及 (生信宝典注:不是很少有帖子涉及,而是看有没有发现,配对的个体信息可以视作一个批次因素,在DESeq2差异基因分析和批次效应移除, 高通量数据中批次效应的鉴定和处理 - ...
  • 高颜值绘图工具ImageGP不只可以绘制常见的生信图形,还提供了不少分析功能,如之前提到的在线WGCNA分析,操作简单,参数丰富,结果交互,绝对是科研必备神器。访问也很简单,百度或谷歌搜...
  • 转录组入门7-用DESeq2进行差异表达分析 Analyzing RNA-seq data with DESeq2 1.关于DESeq2的概述 A basic task in the analysis of count data from RNA-seq is the detection of differentially expressed genes. ...
  • 消除文库大小不同,测序深度对差异分析结果的影响 怎样标准化? 找到一个能反映文库大小的因子,利用这个因子对rawdata进行标准化 3 根据模型检验求p value:泊松分布(poisson distribution)、负二项式分布(NB)等.....
  • 转录组数据中差异基因的相关性分析
  • 手把手教学差异表达基因分析

    千次阅读 2020-09-19 23:18:00
    文章目录引言安装并导入DESeq2包数据要求制作dds对象,进行差异分析筛选差异基因完整代码 引言 对于组学分析来说,常常会寻找组间的差异,例如差异基因(转录组)、差异菌(宏基因组)以及差异通路(宏基因组),而...
  • 基因芯片分析急性心肌梗死患者外周血差异基因表达谱.pdf

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 12,324
精华内容 4,929
关键字:

差异基因分析