中文字幕理论片,69视频免费在线观看,亚洲成人app,国产1级毛片,刘涛最大尺度戏视频,欧美亚洲美女视频,2021韩国美女仙女屋vip视频

差異基因表達(dá)分析

我按照前面的流程轉(zhuǎn)錄組差異表達(dá)分析小實(shí)戰(zhàn)（一），將小鼠的4個(gè)樣本又重新跑了一遍，從而獲得一個(gè)新的count文件：mouse_all_count.txt，有需要的話，可以下載下來(lái)進(jìn)行后續(xù)的差異分析。

一般來(lái)說(shuō)，由于普遍認(rèn)為高通量的read count符合泊松分布，所以一些差異分析的R包都是基于負(fù)二項(xiàng)式分布模型的，比如DESeq、DESeq2和edgeR等，所以這3個(gè)R包從整體上來(lái)說(shuō)是類似的（但各自標(biāo)準(zhǔn)化算法是不一樣的）。

當(dāng)然還有一個(gè)常用的R包則是Limma包，其中的limma-trend和limma-voom都能用來(lái)處理RNA-Seq數(shù)據(jù)（但對(duì)應(yīng)適用的情況不一樣）

下面準(zhǔn)備適用DESeq2和edgeR兩個(gè)R包分別對(duì)小鼠的count數(shù)據(jù)進(jìn)行差異表達(dá)分析，然后取兩者的結(jié)果的交集的基因作為差異表達(dá)基因。

1. DEseq2

   library(DESeq2)##數(shù)據(jù)預(yù)處理database <- read.table(file = "mouse_all_count.txt", sep = "\t", header = T, row.names = 1)database <- round(as.matrix(database))##設(shè)置分組信息并構(gòu)建dds對(duì)象condition <- factor(c(rep("control",2),rep("Akap95",2)), levels = c("control", "Akap95"))coldata <- data.frame(row.names = colnames(database), condition)dds <- DESeqDataSetFromMatrix(countData=database, colData=coldata, design=~condition)dds <- dds[ rowSums(counts(dds)) > 1, ]##使用DESeq函數(shù)估計(jì)離散度，然后差異分析獲得res對(duì)象dds <- DESeq(dds)res <- results(dds)#最后設(shè)定閾值，篩選差異基因，導(dǎo)出數(shù)據(jù)(全部數(shù)據(jù)。包括標(biāo)準(zhǔn)化后的count數(shù))res <- res[order(res$padj),]diff_gene <- subset(res, padj < 0.05 & (log2FoldChange > 1 | log2FoldChange < -1))diff_gene_DESeq2 <- row.names(diff_gene)resdata <- merge(as.data.frame(res), as.data.frame(counts(dds, normalized=TRUE)),by="row.names",sort=FALSE)write.csv(resdata,file = "control_vs_Akap95.csv",row.names = F)

   最終獲得572個(gè)差異基因（篩選標(biāo)準(zhǔn)為padj < 0.05, |log2FoldChange| > 1）

2. edgeR

   library(edgeR)##跟DESeq2一樣，導(dǎo)入數(shù)據(jù)，預(yù)處理（用了cpm函數(shù)）exprSet<- read.table(file = "mouse_all_count.txt", sep = "\t", header = TRUE, row.names = 1, stringsAsFactors = FALSE)group_list <- factor(c(rep("control",2),rep("Akap95",2)))exprSet <- exprSet[rowSums(cpm(exprSet) > 1) >= 2,]##設(shè)置分組信息，并做TMM標(biāo)準(zhǔn)化exprSet <- DGEList(counts = exprSet, group = group_list)exprSet <- calcNormFactors(exprSet)##使用qCML（quantile-adjusted conditional maximum likelihood）估計(jì)離散度（只針對(duì)單因素實(shí)驗(yàn)設(shè)計(jì)）exprSet <- estimateCommonDisp(exprSet)exprSet <- estimateTagwiseDisp(exprSet)##尋找差異gene(這里的exactTest函數(shù)還是基于qCML并且只針對(duì)單因素實(shí)驗(yàn)設(shè)計(jì))，然后按照閾值進(jìn)行篩選即可et <- exactTest(exprSet)tTag <- topTags(et, n=nrow(exprSet))diff_gene_edgeR <- subset(tTag$table, FDR < 0.05 & (logFC > 1 | logFC < -1))diff_gene_edgeR <- row.names(diff_gene_edgeR)write.csv(tTag$table,file = "control_vs_Akap95_edgeR.csv")

   最終獲得688個(gè)差異基因（篩選標(biāo)準(zhǔn)為FDR < 0.05, |log2FC| > 1）

3. 取DESeq2和edgeR兩者結(jié)果的交集

   diff_gene <- diff_gene_DESeq2[diff_gene_DESeq2 %in% diff_gene_edgeR]

   最終的差異基因數(shù)目為545個(gè)

   head(diff_gene)[1] "ENSMUSG00000003309.14" "ENSMUSG00000046323.8"  "ENSMUSG00000001123.15"[4] "ENSMUSG00000023906.2"  "ENSMUSG00000044279.15" "ENSMUSG00000018569.12"

4. 其他兩個(gè)R包（DESeq和limma）就不在這嘗試了，我之前做過(guò)對(duì)于這4個(gè)R包的簡(jiǎn)單使用筆記，可以參考下：
   簡(jiǎn)單使用DESeq做差異分析
   簡(jiǎn)單使用DESeq2/EdgeR做差異分析
   簡(jiǎn)單使用limma做差異分析

GO&&KEGG富集分析

以前一直沒(méi)有機(jī)會(huì)用Y叔寫的clusterProfiler包，這次正好看說(shuō)明用一下。

1. GO富集，加載clusterProfiler包和org.Mm.eg.db包（小鼠嘛），然后將ENSEMBL ID后面的版本號(hào)去掉，不然后面不識(shí)別這個(gè)ID，然后按照clusterProfiler包的教程說(shuō)明使用函數(shù)即可。

   library(clusterProfiler)library(org.Mm.eg.db)##去除ID的版本號(hào)diff_gene_ENSEMBL <- unlist(lapply(diff_gene, function(x){strsplit(x, "\\.")[[1]][1]}))##GOid mapping + GO富集ego <- enrichGO(gene = diff_gene_ENSEMBL, OrgDb = org.Mm.eg.db,         keytype = "ENSEMBL", ont = "BP", pAdjustMethod = "BH",         pvalueCutoff = 0.01, qvalueCutoff = 0.05)##查看富集結(jié)果數(shù)據(jù)enrich_go <- as.data.frame(ego)##作圖barplot(ego, showCategory=10)dotplot(ego)enrichMap(ego)plotGOgraph(ego)

2. KEGG富集，首先需要將ENSEMBL ID轉(zhuǎn)化為ENTREZ ID，然后使用ENTREZ ID作為kegg id，從而通過(guò)enrichKEGG函數(shù)從online KEGG上抓取信息，并做富集

   library(clusterProfiler)library(org.Mm.eg.db)##ID轉(zhuǎn)化ids <- bitr(diff_gene_ENSEMBL, fromType = "ENSEMBL", toType = "ENTREZID", OrgDb = "org.Mm.eg.db")kk <- enrichKEGG(gene = ids[,2], organism = "mmu", keyType = "kegg",         pvalueCutoff = 0.05, pAdjustMethod = "BH", qvalueCutoff = 0.1)##查看富集結(jié)果數(shù)據(jù)enrich_kegg <- as.data.frame(kk)##作圖dotplot(kk)

到這里為止，轉(zhuǎn)錄組的差異表達(dá)分析算是做完了，簡(jiǎn)單的來(lái)說(shuō)，這個(gè)過(guò)程就是將reads mapping 到reference上，然后計(jì)數(shù)獲得count數(shù)，然后做差異分析，最后來(lái)個(gè)GO KEGG，over了。。。

對(duì)于mapping和計(jì)數(shù)這兩部還有其實(shí)還有好多軟件，具體可見(jiàn)文獻(xiàn)：Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis，有時(shí)間可以都嘗試下。

至于GO && KEGG這兩步，對(duì)于人、小鼠等模式物種來(lái)說(shuō)，不考慮方便因素來(lái)說(shuō)，完全可以自己寫腳本來(lái)完成，數(shù)據(jù)可以從gene ontology官網(wǎng)下載，然后就是GO id與gene id相互轉(zhuǎn)化了。KEGG 也是一樣，也可以用腳本去KEGG網(wǎng)站上自行抓取，先知道URL規(guī)則，然后爬數(shù)據(jù)即可。