上一步驟得到了表達矩陣，兩個樣本分別是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