很明显，得到了表达矩阵之后，根据上面的样本信息，可以按照年龄，性别，取样部位来进行分组找差异。
可以参考：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；第三条看不懂；第四条建议是基于大量的稳定的参照转录本。（PS：看不懂原理这里的原理，，看用的多是第二条建议，那就试试第二个吧）
###下载安装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.
下面是无重复样本计算差异的例子：
发现gfold不同版本输入文件格式不同，如果是需要输入文件5列，可以参考这里https://www.jianshu.com/p/50cd51c090eb
对于前面得到的counts列表（hisat_matrix.out）每个样本单独分开，并命名为samplename.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
# -suf:后面加后缀
#也可以不加后缀，以上代码等同于gfold diff -s1 F.OC.read_cnt -s2 M.OC.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个基因重合。
注释
接下来注释这些差异基因分别包含的mRNA和lincRNA，及其GO和KEGG pathway分析。
mRNA和lincRNA的注释
写个脚本从注释文件里提取出来差异基因是什么分类就可以了。
GO和KEGG注释
GO和KEGG注释用Y叔的clusterprofiler。首先下载clusterProfiler包及猕猴的OrgDb，目前Bioconductor共包含20个物种。
##加载clusterProfiler及OrgDb
source("https://bioconductor.org/biocLite.R")
biocLite("clusterProfiler")
biocLite("org.Mmu.eg.db")
library(clusterProfiler)
library("org.Mmu.eg.db")
setwd("G:/My_exercise/cluster_profile/")
##读入差异基因
gene_diff <- read.csv(file="gfold.all.csv",header = TRUE,sep = ",")
dim(gene_diff)
gene <- gene_diff$X.GeneSymbol
clusterProfiler提供了ID转化的函数bitr(), 25种ID格式可以相互转化。对于GO注释不需要进行ID转换，KEGG分析时需要先进行ID转换，ID的类型（fromType & toType）必须有一个是kegg id, ncbi-geneid,ncbi-preoteinid 中的一个。
For GO analysis, user don’t need to convert ID, all ID type provided by OrgDb can be used in groupGO, enrichGO and gseGO by specifying keytype parameter.
bitr_kegg: converting biological IDs using KEGG API,The ID type (both fromType & toType) should be one of ‘kegg’, ‘ncbi-geneid’, ‘ncbi-proteinid’ or ‘uniprot’. The ‘kegg’ is the primary ID used in KEGG database.
## bitr()ID转换
gene.df_ID <- bitr(gene,fromType = "ENSEMBL",toType = c("SYMBOL","UNIPROT","ENTREZID"),OrgDb = org.Mmu.eg.db)
write.table(as.data.frame(gene.df_ID),file="gene_transID",sep = "\t",quote = FALSE)
# 获取按照log2FC大小来排序的基因列表
genelist <-gene_diff$log2fdc
#names(genelist) <- rownames(diff_gene_deseq2)
genelist <- sort(genelist, decreasing = TRUE)
参考clusterProfiler文档。
GO分析
具体参数参考：(伪）从零开始学转录组（8）：富集分析
#groupGO(gene, OrgDb, keytype = "ENTREZID", ont = "CC", level = 2,readable = FALSE)
#groupGO()支持ENTREZID
#If readable is setting to TRUE, the input gene IDs will be converted to gene symbols.
#BP:Biological Process; CC:Cellular Componen;MF:Molecular Function
#GO:https://www.zhihu.com/question/53055375
ggo_CC <- groupGO(gene = gene.df_ID$ENTREZID,
OrgDb = org.Mmu.eg.db,
ont = "CC",
level = 3,
readable = TRUE)
ggo_BP <- groupGO(gene = gene.df_ID$ENTREZID,
OrgDb = org.Mmu.eg.db,
ont = "BP",
level = 3,
readable = TRUE)
ggo_MF <- groupGO(gene = gene.df_ID$ENTREZID,
OrgDb = org.Mmu.eg.db,
ont = "MF",
level = 3,
readable = TRUE)
#可视化
barplot(ggo_CC,showCategory=12,font.size=7,title="groupGO Cellular Component")
barplot(ggo_BP,showCategory=12,font.size=7,title="groupGO Biological Process")
barplot(ggo_MF,showCategory=12,font.size=7,title="groupGO Molecular Function")
GO over-representation test
#enrichGO(gene, OrgDb, keytype = "ENTREZID", ont = "MF",
# pvalueCutoff = 0.05, pAdjustMethod = "BH", universe, qvalueCutoff = 0.2,
# minGSSize = 10, maxGSSize = 500, readable = FALSE, pool = FALSE)
ego_CC <- enrichGO(gene = gene,
OrgDb = org.Mmu.eg.db,
keytype = "ENSEMBL",
ont = "CC",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
readable = TRUE)
ego_BP <- enrichGO(gene = gene,
OrgDb = org.Mmu.eg.db,
keytype = "ENSEMBL",
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
readable = TRUE)
ego_MF <- enrichGO(gene = gene,
OrgDb = org.Mmu.eg.db,
keytype = "ENSEMBL",
ont = "MF",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
readable = TRUE)
ego_all <- enrichGO(gene = gene,
OrgDb = org.Mmu.eg.db,
keytype = "ENSEMBL",
ont = "ALL",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
readable = TRUE,
pool = TRUE)
#可视化barplot
barplot(ego_CC, showCategory=12,font.size=7,title="EnrichmentGO_CC")
barplot(ego_BP, showCategory=12,font.size=7,title="EnrichmentGO_BP")
barplot(ego_MF, showCategory=12,font.size=7,title="EnrichmentGO_MF")
barplot(ego_all, showCategory=12,font.size=7,title="EnrichmentGO_all")
#dotplot
dotplot(ego_CC, showCategory=12,font.size=7,title="EnrichmentGO_CC")
dotplot(ego_BP, showCategory=12,font.size=7,title="EnrichmentGO_BP")
dotplot(ego_MF, showCategory=12,font.size=7,title="EnrichmentGO_MF")
#enrichMap
#plotGOgraph(ego_MF)
enrichMap(ego_CC)
enrichMap(ego_BP)
enrichMap(ego_MF)
#cnetplot
cnetplot(ego_CC, categorySize="pvalue", foldChange=genelist)
cnetplot(ego_BP, categorySize="pvalue", foldChange=genelist)
cnetplot(ego_MF, categorySize="pvalue", foldChange=genelist)
##############################################
####GO Gene Set Enrichment Analysis
##gseGO(geneList, ont = "BP", OrgDb, keyType = "ENTREZID", exponent = 1,
## nPerm = 1000, minGSSize = 10, maxGSSize = 500, pvalueCutoff = 0.05,
## pAdjustMethod = "BH", verbose = TRUE, seed = FALSE, by = "fgsea")
#ego_gsea_CC <- gseGO(geneList = genelist,
# OrgDb = org.Mmu.eg.db,
# keyType = "ENSEMBL",
# ont = "CC")
#dotplot(ego_gsea_CC)
KEGG（pathway）分析
kegg <- gene.df_ID[,4]
ekk <- enrichKEGG(kegg, keyType = "kegg",organism = "mcc", pvalueCutoff = 0.05, pAdjustMethod = "BH", qvalueCutoff = 0.1)
head(summary(ekk))
mkk <- enrichMKEGG(kegg,organism = 'mcc')
##可视化
dotplot(ekk,title="enrichKEGG")
cnetplot(ekk, categorySize="pvalue", foldChange=genelist)
write.table(as.data.frame(ekk),file = "kegg",sep="\t",quote=F,row.names=F)
write.csv(as.data.frame(ekk),file = "kegg.csv")
参考资料：
一步法差异分析：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
