載入表達(dá)矩陣，然后設(shè)置好分組信息
用DEseq2進(jìn)行差異分析，也可以走走edgeR或者limma的voom流程
基本任務(wù)是得到差異分析結(jié)果，進(jìn)階任務(wù)是比較多個(gè)差異分析結(jié)果的異同點(diǎn)。
了解差異基因分析背后的統(tǒng)計(jì)學(xué)基礎(chǔ)（數(shù)據(jù)標(biāo)準(zhǔn)化，假設(shè)檢驗(yàn)）
簡(jiǎn)單回顧之前genek-轉(zhuǎn)錄組原理篇的知識(shí)點(diǎn)，參考學(xué)習(xí)徐洲更的轉(zhuǎn)錄組分析流程七（差異基因分析），及相關(guān)統(tǒng)計(jì)學(xué)書(shū)本知識(shí)。
(處理樣品中，有一個(gè)基因的表達(dá)量異常的高，大量的reads被占，導(dǎo)致其他gene的相對(duì)表達(dá)量下降。)TPM的標(biāo)準(zhǔn)化是在 各自的==樣品內(nèi)部==的 差異基因==相對(duì)表達(dá)量==對(duì)于組間來(lái)說(shuō)，也就是處理組和對(duì)照組來(lái)說(shuō)就需要進(jìn)一步的標(biāo)準(zhǔn)化。目的就是能夠得到==樣品間的絕對(duì)表達(dá)量==是否有差異。
標(biāo)準(zhǔn)化方法：
內(nèi)參基因，看家基因。某一基因在樣本中的表達(dá)肯定是一樣的。(有限制，依賴基因的功能注釋，數(shù)目少準(zhǔn)確不高。)
通用方法：假設(shè)大多數(shù)基因是沒(méi)有差異表達(dá)的。統(tǒng)計(jì)學(xué)方法找到標(biāo)準(zhǔn)化的因子。
DESeq2(estimateSizeFactors/sizeFactors)、edgeR（TMM,calcNormFactors）。差異表達(dá)分析鑒定軟件。
輸入文件是reads count矩陣。
可用Trinity軟件包里的run_DE_analysis.pl
FPKM/TPM/TMM的選擇
組內(nèi)比較 TPM（不同基因在同一樣品中的表達(dá)差異）
組間比較 TMM（同一基因在不同樣品間的表達(dá)差異）
低表達(dá)基因的過(guò)濾TMM-normalized FPKM(用于數(shù)據(jù)過(guò)濾、heatmap、共表達(dá)等)
TMM-TPM還未結(jié)合起來(lái)（可能TMM矯正TPM后意義就變了？或者TMM還未流行起來(lái)）
建庫(kù)測(cè)序是一個(gè)隨機(jī)抽樣的過(guò)程?？v使是兩次測(cè)序，即抽樣，同一基因的reads_count也會(huì)有差異。
當(dāng)觀察到的基因在兩組reads_count 中不一樣的時(shí)候，可能有兩種原因：可能是==表達(dá)豐度上的差異，也可能是隨機(jī)抽樣過(guò)程中的波動(dòng)==。假設(shè)是隨機(jī)波動(dòng)導(dǎo)致的，然后計(jì)算在隨機(jī)波動(dòng)的前提下，出現(xiàn)當(dāng)前事件的概率，如果概率低。 反證法。
假設(shè)檢驗(yàn)：需要根據(jù)實(shí)際情況（總體均值、方差是否已知，單樣本還是多樣本，重復(fù)個(gè)數(shù)的多少/大樣本量還是小樣本量）
差異表達(dá)基因通常用t檢驗(yàn)（組間差異除以/組內(nèi)差異，組間差異更大，而組內(nèi)差異小則更可能有差異【此即是需要測(cè)更多生物學(xué)重復(fù)，能更準(zhǔn)確的估計(jì)組內(nèi)差異】）
t值的雙側(cè)檢驗(yàn)和單側(cè)檢驗(yàn)。
單側(cè)檢驗(yàn)：G1在處理組中表達(dá) 大于 對(duì)照組。右側(cè)面積->則單側(cè)檢驗(yàn)，即upregulate。
雙側(cè)檢驗(yàn)：G1處理組與對(duì)照組有明顯差異。則雙側(cè)檢驗(yàn)
一類錯(cuò)誤(即為假陽(yáng)性，與P-value對(duì)應(yīng)，0.05 or 0.01)，二類錯(cuò)誤(指漏診率)
降低一類錯(cuò)誤：對(duì)于一條gene可能有5%可能誤診，而在差異基因鑒定時(shí)都是大規(guī)模的，如果一萬(wàn)條，50條則太高了。p-value 進(jìn)一步的矯正 FDR方法，q-value為升級(jí)版。
降低二類錯(cuò)誤：G3增加生物學(xué)重復(fù)，G2增加測(cè)序量來(lái)提高抽樣次數(shù)。
source("https://bioconductor.org/biocLite.R") 載入安裝工具bioconductorbiocLite("DESeq2") bioconductor安裝DESEQ2library(DESeq2)對(duì)于DESeq2需要輸入的三個(gè)數(shù)據(jù)：表達(dá)矩陣、樣品信息矩陣、差異比較矩陣
而對(duì)于DESeq2的差異分析步驟，總結(jié)起來(lái)就是三步：
<font color=darkred>構(gòu)建一個(gè)dds(DESeqDataSet)的對(duì)象</font>；
<font color=darkred>利用DESeq函數(shù)進(jìn)行標(biāo)準(zhǔn)化處理</font>；
<font color=darkred>用result函數(shù)來(lái)提取差異比較的結(jié)果</font>。
構(gòu)建dds矩陣的基本代碼
dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design= ~ batch + condition)
其輸入的三個(gè)文件：
表達(dá)矩陣：即代碼中的countData，就是通過(guò)之前的HTSeq-count生成的reads-count計(jì)算融合的矩陣。行為基因名，列為樣品名，值為reads或fragment的整數(shù)。
樣品信息矩陣：即上述代碼中的colData，它的類型是一個(gè)dataframe（數(shù)據(jù)框），第一列是樣品名稱，第二列是樣品的處理情況（對(duì)照還是處理等）?？梢詮谋磉_(dá)矩陣中導(dǎo)出或是自己?jiǎn)为?dú)建立。
差異比較矩陣：即代碼中的design，差異比較矩陣就是告訴差異分析函數(shù)哪些是對(duì)照，哪些是處理。
正式構(gòu)建dds的矩陣：
輸入HTSeq-count的各個(gè)結(jié)果文件
KD1 <- read.table("mus_E14_cell_Akap95_shRNA_rep1_SRR3589959.count",sep = "\t",col.names = c("gene_id","rep1"))KD2 <- read.table("mus_E14_cell_Akap95_shRNA_rep2_SRR3589960_matrix.count",sep = "\t",col.names = c("gene_id","rep2"))control1 <- read.table("mus_E14_cell_control_shRNA_rep1_SRR3589961_matrix.count",sep = "\t",col.names = c("gene_id","control1"))control2 <- read.table("mus_E14_cell_control_shRNA_rep2_SRR3589962_matrix.count",sep = "\t",col.names = c("gene_id","control2"))merge_file
raw_count <- merge(merge(merge(KD1,KD2,by="gene_id"),control1,by="gene_id"),control2,by="gene_id")
構(gòu)建DESeq_data_set基本信息
count_data <- raw_count[,2:5]row.names(count_data) <- raw_count[,1]condition <- factor(c("rep","rep","condition","condition"))col_data <- data.frame(row.names = colnames(count_data), condition)dds <- DESeqDataSetFromMatrix(countData = count_data,colData = col_data,design = ~ condition)過(guò)濾低質(zhì)量的低count數(shù)據(jù)：
nrow(dds) ## 顯示52636行dds_filter <- dds[ rowSums(counts(dds))>1, ] ##將所有樣本基因表達(dá)量之和小于1的基因過(guò)濾掉nrow(dds_filter) ##顯示27101行,過(guò)濾掉了近50%
使用DESeq進(jìn)行差異表達(dá)分析：(一堆統(tǒng)計(jì)統(tǒng)計(jì)學(xué)知識(shí)的坑待填啊)
estimation of size factors（estimateSizeFactors)，
estimation of dispersion（estimateDispersons)，
Negative Binomial GLM fitting and Wald statistics（nbinomWaldTest）
DESeq自動(dòng)化的基本過(guò)程如圖：
dds_out <- DESeq(dds_filter)res <- results(dds_out)summary(res)
<font color=darkred>結(jié)果顯示:</font> 2.6%的基因：717個(gè)上調(diào)，2.1%的基因579下調(diào)，另外還有41%的基因表達(dá)量是不高的，屬于low count。
在利用數(shù)據(jù)比較分析兩個(gè)樣品中同一個(gè)基因是否存在差異表達(dá)的時(shí)候，一般選取Foldchange值和經(jīng)過(guò)FDR矯正過(guò)后的p值，取padj值(p值經(jīng)過(guò)多重校驗(yàn)校正后的值)小于0.05，log2FoldChange大于1的基因作為差異基因集
table(res$padj<0.05)res_deseq <- res[order(res$padj),]diff_gene_deseq2 <- subset(res_deseq, padj<0.05 & (log2FoldChange > 1 | log2FoldChange < -1))## or diff_gene_deseq2 <- subset(res_desq, padj<0.05 & abs(log2FoldChange) >=1)diff_gene_deseq2 <- row.names(diff_gene_deseq2)res_diff_data <- merge(as.data.frame(res),as.data.frame(counts(dds_out,normalize=TRUE)),by="row.names",sort=FALSE)write.csv(res_diff_data,file = "11_30_mouse_data.csv",row.names = F)
提取結(jié)果：
MA-plot (R. Dudoit et al. 2002) ，也叫 mean-difference plot或者Bland-Altman plot，用來(lái)估計(jì)模型中系數(shù)的分布。 X軸, the “A” （ “average”）；Y軸，the “M” （“minus”） – subtraction of log values is equivalent to the log of the ratio。
M表示log fold change，衡量基因表達(dá)量變化，上調(diào)還是下調(diào)。A表示每個(gè)基因的count的均值。根據(jù)summary可知，low count的比率很高，所以大部分基因表達(dá)量不高，也就是集中在0的附近（log2(1)=0，也就是變化1倍）.提供了模型預(yù)測(cè)系數(shù)的分布總覽。
library(ggplot2)plotMA(res,ylim=c(-5,5))topGene <- rownames(res)[which.min(res$padj)]with(res[topGene, ], { points(baseMean, log2FoldChange, col="dodgerblue", cex=2, lwd=2) text(baseMean, log2FoldChange, topGene, pos=2, col="dodgerblue")})lfcShrink 收縮log2 fold change
resLFC <- lfcShrink(dds_out,coef = 2,res=res)plotMA(resLFC, ylim=c(-5,5))topGene <- rownames(resLFC)[which.min(res$padj)]with(resLFC[topGene, ], { points(baseMean, log2FoldChange, col="dodgerblue", cex=2, lwd=2) text(baseMean, log2FoldChange, topGene, pos=2, col="dodgerblue")})idx <- identify(res$baseMean, res$log2FoldChange)基本上是按照PANDA姐的分析的代碼走一遍，很多參數(shù)的含義還不是很清楚，等接下來(lái)的時(shí)間再學(xué)習(xí)作圖相關(guān)。
library(genefilter)library(pheatmap)rld <- rlogTransformation(dds_out,blind = F)write.csv(assay(rld),file="mm.DESeq2.pseudo.counts.csv")topVarGene <- head(order(rowVars(assay(rld)),decreasing = TRUE),20)mat <- assay(rld)[ topVarGene, ]### mat <- mat - rowMeans(mat) 減去一個(gè)平均值，讓數(shù)值更加集中。第二個(gè)圖anno <- as.data.frame(colData(rld)[,c("condition","sizeFactor")])pheatmap(mat, annotation_col = anno)結(jié)果兩個(gè)熱圖：
#### perl 單行：第三列l(wèi)og2foldchange大于1為上調(diào)，小于-1下調(diào)，其他的不顯著。perl -F'\t' -alne '$F[5]=~s/\r//;if(/baseMean/){print "gene\tlog2FoldChange\tpadj\tsignificant"} else{unless(/NA/){if($F[2] >=1 && $F[5]<0.05){print "$F[0]\t$F[2]\t$F[5]\tup"} elsif($F[2]<=-1 && $F[5]<0.05){print "$F[0]\t$F[2]\t$F[5]\tdown"} else{print "$F[0]\t$F[2]\t$F[5]\tno"}}}' mm_deseq2_resLFS.csv >volcano.txt### 火山圖data <-read.table(file="volcano.txt",header = TRUE, row.names =1,sep = "\t")volcano <-ggplot(data = volcano_data,aes(x=log2FoldChange,y= -1*log10(padj)))+geom_point(aes(color=significant))+scale_color_manual(values = c("red","grey","blue")) + labs(title="Volcano_Plot",x=expression((log[2](FC)), y=expression(-log[10](padj)) ))+geom_hline(yintercept=1.3,linetype=4)+geom_vline(xintercept=c(-1,1),linetype=4)volcano11/30下午小結(jié)：本次差異基因分析的學(xué)習(xí)，代碼部分基本上是參考學(xué)習(xí)著各個(gè)大神們的筆記一步步操作來(lái)的。好的方面是總算對(duì)差異基因分析有了一定系統(tǒng)化的了解，對(duì)DESeq2差異分析軟件了解了基本使用方法，對(duì)一些簡(jiǎn)單可視化圖形的展示與解讀有了一定的認(rèn)識(shí)。但仍然有許多不足，許多坑待填比如：R語(yǔ)言，ggplot2可視化方面，還有最大的坑：統(tǒng)計(jì)學(xué)相關(guān)知識(shí)。這些將會(huì)是轉(zhuǎn)錄組學(xué)習(xí)完之后新的學(xué)習(xí)方向。
【PANDA姐的轉(zhuǎn)錄組入門-差異基因分析-2】 https://mp.weixin.qq.com/s/05Gv1pnL0GWZwjz893xcOw
【用DESeq2包來(lái)對(duì)RNA-seq數(shù)據(jù)進(jìn)行差異分析】http://www.bio-info-trainee.com/bioconductor_China/software/DESeq2.html
【（偽）從零開(kāi)始學(xué)轉(zhuǎn)錄組（7）：差異基因表達(dá)分析】https://mp.weixin.qq.com/s?__biz=MzI1MjU5MjMzNA==&mid=2247484514&idx=1&sn=1c6d227c6d7ac432d8baffb93b0b9072&chksm=e9e02dc3de97a4d59d918ee37655153fa4ccc8122e3259c0201ff441c922548f22f84311ad91&scene=21#wechat_redirect
【青山屋主-轉(zhuǎn)錄組入門】https://zhuanlan.zhihu.com/p/30350531?group_id=905527998389362688
【Analyzing RNA-seq data with DESeq2
】http://www.bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html
【轉(zhuǎn)錄組差異表達(dá)篩選的真相
】https://mp.weixin.qq.com/s?__biz=MzI4NjMxOTA3OA==&mid=2247484080&idx=1&sn=7de397add158271aa51763d0d3598deb&chksm=ebdf8a39dca8032f1cb8b8154a024ee7f564c2e33b2f824efce7e34a016ef7baf92d52d64c79&scene=21#wechat_redirect
【FPKM / RPKM與TPM哪個(gè)用來(lái)篩選差異表達(dá)基因更準(zhǔn)確？你可別逗了！】https://mp.weixin.qq.com/s?__biz=MzI4NjMxOTA3OA==&mid=2247483987&idx=1&sn=aa2ca81e7fe128edaaedc47479c517c9&chksm=ebdf8adadca803cc31261a1ccabf8a6bdb835ce12b670cc01969fcfde51cfdb991d0619e7695&scene=21#wechat_redirect
聯(lián)系客服
