lncRNA实战项目-第五步-差异表达的mRNA和lncRNA

上一步骤得到了表达矩阵,两个样本分别是F_1yr.OC和M_1yr.OC, 所以接下来的差异分析就是比较1岁猕猴脑OC区域女性和男性的差别,差异分析的分析方法很多,主要根据前面标准化的方法,有基于counts的差异分析,也有基于标准化后的FPKM,TPM等的差异分析。 常见的R包有(摘自https://github.com/jmzeng1314/my-R/tree/master/DEG_scripts): edgeR (Robinson et al., 2010) DESeq / qDESeq2 (Anders and Huber, 2010, 2014) DEXSeq (Anders et al., 2012) limmaVoom Cuffdiff / Cuffdiff2 (Trapnell et al., 2013) PoissonSeq baySeq 作业里给的参考是一步法差异分析,是对常见的R包做了下封装,包括了对转录组的raw counts数据分析DEseq2包和edgeR包,及对于芯片等normalization好的表达矩阵数据的limma和t.test等。用的时候只要设置好表达矩阵和分组矩阵,然后选择特定的方法,一步就可以进行差异分析。

但是这里的样本是无生物学重复的,无重复的数据做差异分析是一件很麻烦的事,可靠性都不能保证。。。但是目前由于测序的价格,还有样本自身的珍贵稀缺性,部分实验设计仍然是没有生物学重复的。对于无重复样本的差异分析有几种方法可以选择,如edgeR,DEGseq和GFOLD等。下面分别尝试edgeR,DEGseq及GFOLD:

edgeR做无重复样本的差异分析

edgeR针对无重复样本给出了四条建议,第一条建议是仅分析MDS plot和fold changes,不做显著性分析;第二条建议是设置合适的离散度值,然后做个exactTest 或glmFit;第三条看不懂;第四条建议是基于大量的稳定的参照转录本。

###下载安装edgeR包
#source("http://bioconductor.org/biocLite.R")
#biocLite("edgeR")
library("edgeR")
library('ggplot2')

###读取数据
setwd("G:/My_exercise/edgeR")
rawdata <- read.table("hisat_matrix.out",header=TRUE,row.names=1,check.names = FALSE)
head(rawdata)
#重命名列名
names(rawdata) <- c("F.1yr.OC.count","M.1yr.OC.count")
#进行分组
group <- factor(c("F.1yr.OC.count","M.1yr.OC.count"))

###过滤与标准化
y <- DGEList(counts=rawdata,genes=rownames(rawdata),group = group)
###TMM标准化
y<-calcNormFactors(y)
y$samples
###推测离散度,根据经验设置,若样本是人,设置bcv = 0.4,模式生物设置0.1.(这里没有经验,我就多试几个)
#bcv <- 0.1
bcv <- 0.2
#bcv <- 0.4
et <- exactTest(y, dispersion=bcv^2)
topTags(et)
summary(de <- decideTestsDGE(et))
###图形展示检验结果
png('F_1yr.OC-vs-M_yrM.OC.png')
detags <- rownames(y)[as.logical(de)];
plotSmear(et, de.tags=detags)
abline(h=c(-4, 4), col="blue");
dev.off()

###导出数据
DE <- et$table
DE$significant <- as.factor(DE$PValue<0.05 & abs(DE$logFC) >1)
write.table(DE,file="edgeR_all2",sep="\t",na="NA",quote=FALSE)
write.csv(DE, "edgeR.F-vs-M.OC2.csv")

#DE2 <- topTags(et,20000)$table
#DE2$significant <- as.factor(DE2$PValue<0.05 & DE2$FDR<0.05 & abs(DE2$logFC) >1)
#write.csv(DE2, "F_1yr.OC-vs-M_1yr.OC3.csv")

edgeR

DEGseq对无重复样本差异分析

也有推荐DEGSeq 中MARS方法的(MARS: MA-plot-based method with Random Sampling model)。

## 1.安装DEGseq
source("https://bioconductor.org/biocLite.R")
biocLite("DEGseq")
library("DEGseq")
setwd("G:/My_exercise/DEG/")
#读入数据,每组样本构建单独一个矩阵
matrix1 <- readGeneExp(file="hisat_matrix.out", geneCol=1, valCol=2)
matrix2 <- readGeneExp(file="hisat_matrix.out", geneCol=1, valCol=3)
DEGexp(geneExpMatrix1=matrix1, geneCol1=1, expCol1=2, groupLabel1="F_1yr.OC",
               geneExpMatrix2=matrix2, geneCol2=1, expCol2=2, groupLabel2="M_1yr.OC",
               method="MARS", outputDir="G:/My_exercise/DEG/")

MA.plot

GFOLD对无重复样本进行差异分析

该软件称尤其适合做无重复样本的差异分析,他对foldchange 的计算考虑到posterior distribution,即克服了pvalue评估显著性的缺点,同时也克服了 fold change 在评估低counts 数的gene时的缺点。

下载软件:
wget https://bitbucket.org/feeldead/gfold/get/e78560195469.zip
unzip e78560195469.zip
cd feeldead-gfold-e78560195469 
#查看REDEME安装说明

安装GFOLD时,需要先安装gsl,然后再编译安装gfold。

#安装gsl
wget ftp://ftp.gnu.org/gnu/gsl/gsl-2.4.tar.gz
tar zxf gsl-2.4.tar.gz
cd gsl-2.4
./configure
make 
make install
#查看帮助文档
cd doc/
firefox gfold.html &
该软件的功能包括5部分:

1)Count reads and rank genes; 2)Count reads; 3)Identify differentially expressed genes without replicates; 4)Identify differentially expressed genes with replicates; 5)Identify differentially expressed genes with replicates only in one condition. 下面是无重复样本计算差异的例子:

对于前面得到的counts列表(hisat_matrix.out)每个样本单独分开,并命名为samplename.read_cnt(一定要加后缀.read_cnt).

awk '{print $1,$2}' OFS='\t' hisat_matrix.out >F.OC.read_cnt
awk '{print $1,$3}' OFS='\t' hisat_matrix.out >M.OC.read_cnt

这里查看下F.OC.read_cnt是否有头文件,若有最好注释掉,否则后面差异结果有错位。然后用gfold diff 一步就可以求出差异基因。输出文件包含4列,第一列GeneID, 第二列是gfold值,gfold值的正负对应着基因的上调和下调,gfold=0认为是无差异的,E-FDR对无重复样本总是1,第四列是log2fold change。

gfold diff -s1 F.OC -s2 M.OC -suf .read_cnt -o F_M.OC.diff
awk '{if($2>0 && $3=1) print $0}' F_M.OC.diff OFS='\t' > up_diff.gene
awk '{if($2<0 && $3=1) print $0}' F_M.OC.diff  OFS='\t' > down_diff.gene
#筛选差异倍数为2
awk '{if($2>1 && $3=1) print $0}' F_M.OC.diff OFS='\t' > up_diff.gene_2
awk '{if($2<-1 && $3=1) print $0}' F_M.OC.diff  OFS='\t' > down_diff.gene_2

上调基因:4324,下调基因:4240,差异变化阈值设置gfold为1时,上调的基因有83个,下调有97个。

差异基因初步统计

用edgeR共筛选到1322个差异显著基因(筛选条件:PValue<0.05 & abs(logFC) >1); 用DEGseq共筛选到743个差异显著基因(筛选条件:abs(log2(Fold_change) normalized ) >1 & p-value < 0.05 & q-value(Storey et al. 2003) <0.05 & Signature(p-value < 0.001)=TRUE), 用GFOLD共筛选到180个差异基因(gfold>1 && gfold<-1,E- FDR=1)。其中gfold筛选到的180个基因全部包含在edgeR和DEGSeq中,edgeR和DEGseq筛选到显著差异基因共有720个基因重合。


参考资料:

一步法差异分析:https://github.com/jmzeng1314/my-R/tree/master/DEG_scripts 从零开始学转录组(7):差异基因表达分析 从零开始学转录组(8):富集分析 RNA-seq项目设计:生物学重复和单个样本测序量对结果的影响 clusterProfiler参考文档 差异基因分析 文献:Efficient experimental design and analysis strategies for the detection of differential expression using RNA-Sequencing

编辑:jimmy

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2018-02-16

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏机器之心

深度 | 生物神经网络基础:从尖峰放电神经元谈起

1822
来自专栏大数据挖掘DT机器学习

车辆目标检测

4394
来自专栏PPV课数据科学社区

【学习】SPSS探索分析实践操作

SPSS为我们提供了探索分析,所谓探索分析之所以是探索,是因为有时候我们对于变量的分布特点不是很清楚,探索的目的在于帮助我们完成以下的工作:识别...

3618
来自专栏专知

机器翻译新时代:Facebook 开源无监督机器翻译模型和大规模训练语料

【导读】基于深度学习的机器翻译往往需要数量非常庞大的平行语料,这一前提使得当前最先进的技术无法被有效地用于那些平行语料比较匮乏的语言之间。为了解决这一问题,Fa...

48011
来自专栏CreateAMind

Keras和DDPG玩赛车游戏(自动驾驶)

这里,s是状态,a是行为/动作,θ是策略网络的模型参数,π是常见的表示策略的符号。我们可以设想策略是我们行为的代理人,即一个从状态到动作的映射函数。

2492
来自专栏大数据挖掘DT机器学习

详细步骤:用R语言做文本挖掘

目录 Part1 安装依赖包 Part2 分词处理 Part3文本聚类 Part4 文本分类 Part5情感分析 Part1 安装依赖包 R语言中中文分析的...

64112
来自专栏章鱼的慢慢技术路

解救小哈——DFS算法举例

2398
来自专栏生信技能树

第41周生信文献分享:肝癌复发的CpG甲基化信号特征

前面我们讲解了一篇2013年多组学数据探索乳腺癌细胞系药物敏感性使用的也是两个机器学习算法,不过是LS-SVM和RF,但是也有借鉴意义。

1182
来自专栏生信宝典

WGCNA分析,简单全面的最新教程

6692
来自专栏机器学习人工学weekly

机器学习人工学weekly-2018/3/17

1. PyTorch构架分析 PyTorch – Internal Architecture Tour 链接:http://blog.christianper...

3037

扫码关注云+社区