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

相关文章

来自专栏HansBug's Lab

2005: [Noi2010]能量采集

2005: [Noi2010]能量采集 Time Limit: 10 Sec  Memory Limit: 552 MB Submit: 1831  Solve...

2726
来自专栏机器之心

教程 | 无需复杂深度学习算法,基于计算机视觉使用Python和OpenCV计算道路交通

4428
来自专栏社区的朋友们

Steering Behaviors 详解

Steering Behaviors 意在使游戏中的AI个体具备真实的运动行为,通过对力的施加与整合,使游戏个体具备类生命体般的运动特征。

5541
来自专栏余林丰

初学数据挖掘——相似性度量(一)

  好久没有写这个了。也就是在去年到今年这个时间段里,同时决定好几件事情。第一:考研。第二:以后方向就是大数据或者是叫数据挖掘。这两件事当然是有联系的,第一件事...

1748
来自专栏潇涧技术专栏

Things of Math

因为近期换了博客主题,对Latex的支持较弱,而且以后可能会很少写和数学有关的内容,所以下线了之前数学专题下的所有文章,但竟然有网友评论希望重新上线,我还以为那...

491
来自专栏编程

基于micropython的数字滤波器

刚参加了昨天的硕士研究生考试,专业课考的信号与系统,报考学校今年出题出的有点偏,不是题型偏而是考察知识点有明显的偏重,简单说考纲里所罗列的考点最多只考了百分之三...

1626
来自专栏WOLFRAM

Mathematica 谜中智 | 赏九美图 戏九连环【谜底篇】

1756
来自专栏生信宝典

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

2972
来自专栏机器之心

教程 | 遗传算法的基本概念和实现(附Java实现案例)

选自Medium 作者:MallawaarachchiFollow 机器之心编译 参与:俞云开、蒋思源 基因遗传算法是一种灵感源于达尔文自然进化理论的启发式搜索...

2748
来自专栏技术随笔

LIDC-IDRI肺结节Dicom数据集解析与总结Reference:

7427

扫码关注云+社区