STEP5:差異表達的mRNA和lncRNA

很明顯宏浩，得到了表達矩陣之后柱查，根據(jù)上面的樣本信息廓俭，可以按照年齡，性別唉工，取樣部位來進行分組找差異研乒。
可以參考：https://github.com/jmzeng1314/my-R/tree/master/DEG_scripts

上一步驟得到了表達矩陣，兩個樣本分別是F_1yr.OC和M_1yr.OC, 所以接下來的差異分析就是比較1歲獼猴腦OC區(qū)域女性和男性的差別酵紫，差異分析的分析方法很多，主要根據(jù)前面標(biāo)準(zhǔn)化的方法错维，有基于counts的差異分析奖地，也有基于標(biāo)準(zhǔn)化后的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
作業(yè)里給的參考是一步法差異分析参歹，是對常見的R包做了下封裝，包括了對轉(zhuǎn)錄組的raw counts數(shù)據(jù)分析的DEseq2包和edgeR包隆判，及對于芯片等normalization好的表達矩陣數(shù)據(jù)的limma和t.test等犬庇。用的時候只要設(shè)置好表達矩陣和分組矩陣，然后選擇特定的方法侨嘀，一步就可以進行差異分析臭挽。

但是這里的樣本是無生物學(xué)重復(fù)的，而無生物學(xué)重復(fù)對差異基因的檢出率和結(jié)果的可靠性都有影響咬腕。目前由于測序的價格及樣本自身的珍貴稀缺性欢峰，部分實驗設(shè)計仍然是沒有生物學(xué)重復(fù)的。對于無重復(fù)樣本的差異分析可以選擇；edgeR纽帖，DEGseq和GFOLD等宠漩。

下面分別嘗試edgeR，DEGseq及GFOLD：

edgeR做無重復(fù)樣本的差異分析

edgeR針對無重復(fù)樣本給出了四條建議懊直，第一條建議是僅分析MDS plot和fold changes扒吁，不做顯著性分析；第二條建議是設(shè)置合適的離散度值室囊，然后做個exactTest 或glmFit雕崩；第三條看不懂；第四條建議是基于大量的穩(wěn)定的參照轉(zhuǎn)錄本波俄。（PS：看不懂原理這里的原理晨逝，，看用的多是第二條建議懦铺，那就試試第二個吧）

###下載安裝edgeR包
#source("http://bioconductor.org/biocLite.R")
#biocLite("edgeR")
library("edgeR")
library('ggplot2')

###讀取數(shù)據(jù)
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"))

###過濾與標(biāo)準(zhǔn)化
y <- DGEList(counts=rawdata,genes=rownames(rawdata),group = group)
###TMM標(biāo)準(zhǔn)化
y<-calcNormFactors(y)
y$samples
###推測離散度,根據(jù)經(jīng)驗設(shè)置捉貌，若樣本是人，設(shè)置bcv = 0.4冬念，模式生物設(shè)置0.1.（這里沒有經(jīng)驗趁窃，我就多試幾個）
#bcv <- 0.1
bcv <- 0.2
#bcv <- 0.4
et <- exactTest(y, dispersion=bcv^2)
topTags(et)
summary(de <- decideTestsDGE(et))
###圖形展示檢驗結(jié)果
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()

###導(dǎo)出數(shù)據(jù)
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對無重復(fù)樣本差異分析

也有推薦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/")
#讀入數(shù)據(jù)急前，每組樣本構(gòu)建單獨一個矩陣
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對無重復(fù)樣本進行差異分析

該軟件稱尤其適合做無重復(fù)樣本的差異分析醒陆，他對foldchange 的計算考慮到posterior distribution，即克服了pvalue評估顯著性的缺點裆针，同時也克服了 fold change 在評估低counts 數(shù)的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.
下面是無重復(fù)樣本計算差異的例子：

發(fā)現(xiàn)gfold不同版本輸入文件格式不同耘婚，如果是需要輸入文件5列罢浇，可以參考這里http://www.reibang.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是否有頭文件沐祷，若有最好注釋掉嚷闭，否則后面差異結(jié)果有錯位。然后用gfold diff 一步就可以求出差異基因赖临。輸出文件包含4列胞锰，第一列GeneID, 第二列是gfold值，gfold值的正負對應(yīng)著基因的上調(diào)和下調(diào)兢榨，gfold=0認為是無差異的胜蛉，E-FDR對無重復(fù)樣本總是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

#篩選差異倍數(shù)為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

上調(diào)基因：4324领突，下調(diào)基因：4240，差異變化閾值設(shè)置gfold為1時案怯，上調(diào)的基因有83個君旦，下調(diào)有97個。

差異基因初步統(tǒng)計

用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）金砍。其中g(shù)fold篩選到的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轉(zhuǎn)化的函數(shù)bitr(), 25種ID格式可以相互轉(zhuǎn)化骆捧。對于GO注釋不需要進行ID轉(zhuǎn)換，KEGG分析時需要先進行ID轉(zhuǎn)換髓绽，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轉(zhuǎn)換
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分析

具體參數(shù)參考：(偽）從零開始學(xué)轉(zhuǎn)錄組（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")
groupgo

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)
ego.dotplot

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")
enrichkegg

參考資料：

一步法差異分析：https://github.com/jmzeng1314/my-R/tree/master/DEG_scripts
從零開始學(xué)轉(zhuǎn)錄組（7）：差異基因表達分析
從零開始學(xué)轉(zhuǎn)錄組（8）：富集分析
RNA-seq項目設(shè)計：生物學(xué)重復(fù)和單個樣本測序量對結(jié)果的影響
clusterProfiler參考文檔
差異基因分析
文獻：Efficient experimental design and analysis strategies for the detection of differential expression using RNA-Sequencing

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 人面猴
    序言：七十年代末顺呕，一起剝皮案震驚了整個濱河市枫攀，隨后出現(xiàn)的幾起案子，更是在濱河造成了極大的恐慌株茶，老刑警劉巖来涨，帶你破解...
    沈念sama閱讀 206,013評論 6贊 481
  • 死咒
    序言：濱河連續(xù)發(fā)生了三起死亡事件，死亡現(xiàn)場離奇詭異忌卤，居然都是意外死亡扫夜，警方通過查閱死者的電腦和手機楞泼，發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 88,205評論 2贊 382
  • 救了他兩次的神仙讓他今天三更去死
    文/潘曉璐 我一進店門驰徊，熙熙樓的掌柜王于貴愁眉苦臉地迎上來，“玉大人堕阔，你說我怎么就攤上這事棍厂。” “怎么了超陆？”我有些...
    開封第一講書人閱讀 152,370評論 0贊 342
  • 道士緝兇錄：失蹤的賣姜人
    文/不壞的土叔 我叫張陵牺弹，是天一觀的道長浦马。 經(jīng)常有香客問我，道長张漂，這世上最難降的妖魔是什么晶默？ 我笑而不...
    開封第一講書人閱讀 55,168評論 1贊 278
  • ?港島之戀（遺憾婚禮）
    正文 為了忘掉前任，我火速辦了婚禮航攒，結(jié)果婚禮上磺陡，老公的妹妹穿的比我還像新娘。我一直安慰自己漠畜，他們只是感情好币他，可當(dāng)我...
    茶點故事閱讀 64,153評論 5贊 371
  • 惡毒庶女頂嫁案：這布局不是一般人想出來的
    文/花漫 我一把揭開白布。 她就那樣靜靜地躺著憔狞，像睡著了一般蝴悉。 火紅的嫁衣襯著肌膚如雪。 梳的紋絲不亂的頭發(fā)上瘾敢，一...
    開封第一講書人閱讀 48,954評論 1贊 283
  • 城市分裂傳說
    那天拍冠，我揣著相機與錄音，去河邊找鬼廉丽。 笑死倦微，一個胖子當(dāng)著我的面吹牛，可吹牛的內(nèi)容都是我干的正压。 我是一名探鬼主播欣福，決...
    沈念sama閱讀 38,271評論 3贊 399
  • 雙鴛鴦連環(huán)套：你想象不到人心有多黑
    文/蒼蘭香墨 我猛地睜開眼，長吁一口氣：“原來是場噩夢啊……” “哼焦履！你這毒婦竟也來了拓劝？” 一聲冷哼從身側(cè)響起，我...
    開封第一講書人閱讀 36,916評論 0贊 259
  • 萬榮殺人案實錄
    序言：老撾萬榮一對情侶失蹤嘉裤，失蹤者是張志新（化名）和其女友劉穎郑临，沒想到半個月后，有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體屑宠，經(jīng)...
    沈念sama閱讀 43,382評論 1贊 300
  • ?護林員之死
    正文 獨居荒郊野嶺守林人離奇死亡厢洞，尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點故事閱讀 35,877評論 2贊 323
  • ?白月光啟示錄
    正文 我和宋清朗相戀三年，在試婚紗的時候發(fā)現(xiàn)自己被綠了典奉。 大學(xué)時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片躺翻。...
    茶點故事閱讀 37,989評論 1贊 333
  • 活死人
    序言：一個原本活蹦亂跳的男人離奇死亡，死狀恐怖卫玖，靈堂內(nèi)的尸體忽然破棺而出公你，到底是詐尸還是另有隱情，我是刑警寧澤假瞬，帶...
    沈念sama閱讀 33,624評論 4贊 322
  • ?日本核電站爆炸內(nèi)幕
    正文 年R本政府宣布陕靠，位于F島的核電站迂尝，受9級特大地震影響，放射性物質(zhì)發(fā)生泄漏剪芥。R本人自食惡果不足惜垄开，卻給世界環(huán)境...
    茶點故事閱讀 39,209評論 3贊 307
  • 男人毒藥：我在死后第九天來索命
    文/蒙蒙 一、第九天 我趴在偏房一處隱蔽的房頂上張望税肪。 院中可真熱鬧说榆，春花似錦、人聲如沸寸认。這莊子的主人今日做“春日...
    開封第一講書人閱讀 30,199評論 0贊 19
  • 一樁弒父案，背后竟有這般陰謀
    文/蒼蘭香墨 我抬頭看了看天上的太陽偏塞。三九已至唱蒸，卻和暖如春，著一層夾襖步出監(jiān)牢的瞬間灸叼，已是汗流浹背神汹。 一陣腳步聲響...
    開封第一講書人閱讀 31,418評論 1贊 260
  • 情欲美人皮
    我被黑心中介騙來泰國打工， 沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留古今，地道東北人屁魏。 一個月前我還...
    沈念sama閱讀 45,401評論 2贊 352
  • 代替公主和親
    正文 我出身青樓，卻偏偏與公主長得像捉腥，于是被迫代替她去往敵國和親氓拼。 傳聞我的和親對象是個殘疾皇子，可洞房花燭夜當(dāng)晚...
    茶點故事閱讀 42,700評論 2贊 345

推薦閱讀更多精彩內(nèi)容