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 条评论
登录 后参与评论

相关文章

来自专栏WOLFRAM

最新北京高考数学试题 Wolfram 完整版

974
来自专栏斑斓

统计学中的相关性分析

掌握一点儿统计学介绍了统计学中常用到的函数,特别重点介绍了Standard Deviation(标准差)。接下来结合一个案例来谈谈相关性(Correlation...

2697
来自专栏数据派THU

教你用机器学习匹配导师 !(附代码)

作者:Zipporah Polinsky-Nagel, Gregory Brucchieri, Marissa Joy, William Kye, Nan Li...

862
来自专栏CSDN技术头条

Bandit算法与推荐系统

推荐系统里面有两个经典问题:EE和冷启动。前者涉及到平衡准确和多样,后者涉及到产品算法运营等一系列。Bandit算法是一种简单的在线学习算法,常常用于尝试解决这...

2899
来自专栏ACM算法日常

各种博弈问题

(一)巴什博奕(Bash Game):只有一堆n个物品,两个人轮流从这堆物品中取物,规定每次至少取一个,最多取m个。最后取光者得胜。

623
来自专栏机器学习之旅

应用:推荐系统-威尔逊区间法

理论上讲,p越大应该越好,但是n的不同,导致p的可信性有差异。100个人投票,50个人投喜欢;10个人投票,6个人喜欢,我们不能说后者比前者要好。

784
来自专栏机器学习算法与Python学习

tweet情感分析流程

关键字全网搜索最新排名 【机器学习算法】:排名第一 【机器学习】:排名第二 【Python】:排名第三 【算法】:排名第四 前言 自然语言处理(NLP)中一个很...

3228
来自专栏Y大宽

Cytoscape插件3:Enrichment Map(2)

所有的微阵列表达数据下载与GEO数据库。Raw.CEL文件用bioconductor的affy包进行RMA。数据集的选择依据以下几个质量控制标准:可靠的并且高覆...

692
来自专栏机器学习之旅

基于Tensorflow实现FFMFFM理论代码实现论文结论总结

没错,这次登场的是FFM。各大比赛中的“种子”算法,中国台湾大学Yu-Chin Juan荣誉出品,美团技术团队背书,Michael Jahrer的论文的fiel...

722
来自专栏机器之心

教程 | 如何保持运动小车上的旗杆屹立不倒?TensorFlow利用A3C算法训练智能体玩CartPole游戏

本教程面向所有对强化学习感兴趣的人,不会涉及太深的机器学习基础,但主题中涵盖了高级策略网络和价值网络的相关知识。此外,我建议阅读 Voldymyr Mnih 的...

453

扫描关注云+社区