lncRNA实战项目-第五步-差异表达的mRNA和lncRNA
响应生信技能树的号召：lncRNA数据分析传送门 , 一起来一个lncRNA数据分析实战！
往期目录
很明显，得到了表达矩阵之后，根据上面的样本信息，可以按照年龄，性别，取样部位来进行分组找差异。
可以参考：https://github.com/jmzeng1314/my-R/tree/master/DEG_scripts
上一步骤得到了表达矩阵，两个样本分别是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")
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/")
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
后续分析，请大家持续关注~