DEseq2篩選差異表達基因并注釋(bioMart)
DESeq2對于輸入數(shù)據(jù)的要求
1.DEseq2要求輸入數(shù)據(jù)是由整數(shù)組成的矩陣。
2.DESeq2要求矩陣是沒有標準化的省艳。
DESeq2包分析差異表達基因簡單來說只有三步：構建dds矩陣，標準化浙巫，以及進行差異分析物邑。
（1）構建dds矩陣
構建dds矩陣需要：
a)表達矩陣：即上述代碼中的count寂殉，就是通過read count計算后并融合生成的矩陣，行為各個基因，列為各個樣品那先，中間為計算reads或者fragment得到的整數(shù)。我們后面要用的是這個文件（readcount.cvs）
b)樣品信息矩陣：即上述代碼中的colData赡艰，它的類型是一個dataframe（數(shù)據(jù)框）售淡，第一列是樣品名稱，第二列是樣品的處理情況（對照還是處理等）慷垮，即condition揖闸，condition的類型是一個factor。這些信息可以從readcount.cvs中導出料身，也可以自己單獨建立汤纸。
c)差異比較矩陣：即design。 差異比較矩陣就是告訴差異分析函數(shù)是要從要分析哪些變量間的差異芹血，簡單說就是說明哪些是對照哪些是處理贮泞。
構建dds矩陣的基本代碼：
dds <- DESeqDataSetFromMatrix(countData, DataFrame(condition), design= ~ condition )
> library(DESeq2)
> mycounts <- read.csv("/media/yanfang/FYWD/RNA_seq/matrix/readcount.csv")
> head(mycounts)
X control1 control2 treat1 treat2
1 ENSG00000000003 807 357 800 963
2 ENSG00000000005 0 0 0 1
3 ENSG00000000419 389 174 405 509
4 ENSG00000000457 288 108 218 283
5 ENSG00000000460 505 208 451 543
6 ENSG00000000938 0 0 0 0
> rownames(mycounts)<-mycounts[,1] #把第一列當作行名來處理
> mycounts<-mycounts[,-1] #把帶X的列刪除
> head(mycounts)
control1 control2 treat1 treat2
ENSG00000000003 807 357 800 963
ENSG00000000005 0 0 0 1
ENSG00000000419 389 174 405 509
ENSG00000000457 288 108 218 283
ENSG00000000460 505 208 451 543
ENSG00000000938 0 0 0 0
> condition <- factor(c(rep("control",2),rep("treat",2)), levels = c("control","treat"))
#condition是因子，不是樣本名稱;對照組和處理組幔烛，各兩個重復啃擦。conditions必須是cloData中的一列
> condition #看一下condition
[1] control control treat treat
Levels: control treat
> colData <- data.frame(row.names=colnames(mycounts), condition)
> colData #看一下colData
condition
control1 control
control2 control
treat1 treat
treat2 treat
> dds <- DESeqDataSetFromMatrix(mycounts, colData, design= ~ condition)
#正式構建dds矩陣
DESeq標準化過程：有這些顯示說明是正常運行
DESeq標準化過程
> dds_norm <- DESeq(dds) #dds標準化
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
> dds_norm #查看標準化后的數(shù)據(jù)
class: DESeqDataSet
dim: 62292 4
metadata(1): version
assays(4): counts mu H cooks
rownames(62292): ENSG00000000003 ENSG00000000005 ... ENSG00000288110 ENSG00000288111
rowData names(22): baseMean baseVar ... deviance maxCooks
colnames(4): control1 control2 treat1 treat2
colData names(2): condition sizeFactor
獲取標準化后的數(shù)據(jù)：
#獲取標準化后的數(shù)據(jù)
> normalized_counts <- counts(dds_norm, normalized=TRUE)
> head(normalized_counts)
control1 control2 treat1 treat2
ENSG00000000003 680.3985 696.3627 670.5922 668.9699158
ENSG00000000005 0.0000 0.0000 0.0000 0.6946728
ENSG00000000419 327.9740 339.4037 339.4873 353.5884602
ENSG00000000457 242.8188 210.6644 182.7364 196.5924052
ENSG00000000460 425.7760 405.7239 378.0464 377.2073357
ENSG00000000938 0.0000 0.0000 0.0000 0.0000000
#根據(jù)基因在不同的樣本中表達變化的差異程度mad值對數(shù)據(jù)排序，差異越大的基因排位越前饿悬。
> normalized_counts_mad <- apply(normalized_counts, 1, mad)
> normalized_counts <- normalized_counts[order(normalized_counts_mad, decreasing=T), ]
# 標準化后的數(shù)據(jù)輸出
> write.table(normalized_counts, file="dds_normalized_counts.xls",
+ quote=F, sep="\t", row.names=T, col.names=T)
# 只在Linux下能運行议惰，目的是填補表格左上角內容
> system(paste("sed -i '1 s/^/ID\t/'", "dds_normalized_counts.xls"))
# 把輸出的標準化后的count取一個log值，log轉換后的結果
> rld <- rlog(dds_norm, blind=FALSE)
> rlogMat <- assay(rld)
> rlogMat <- rlogMat[order(normalized_counts_mad, decreasing=T), ]
#把log轉化后的數(shù)據(jù)再輸出一個表格
> write.table(rlogMat, file="dds_normalized_counts_rlog.xls",
+ quote=F, sep="\t", row.names=T, col.names=T)
#填補表格左上角內容
> system(paste("sed -i '1 s/^/ID\t/'", "dds_normalized_counts_rlog.xls"))
差異基因分析
> res = results(dds_norm, contrast=c("condition", "control", "treat"))
#使用results()函數(shù)獲取結果乡恕，并賦值給res
# contrast: 定義誰和誰比較
> res = res[order(res$pvalue),]
> head(res)
log2 fold change (MLE): condition control vs treat
Wald test p-value: condition control vs treat
DataFrame with 6 rows and 6 columns
baseMean log2FoldChange lfcSE stat pvalue padj
<numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
ENSG00000178691 531.97208257559 2.89986543009505 0.234227480912105 12.3805516705499 3.33045617928295e-35 2.57277739849608e-31
ENSG00000135535 1195.85832961519 1.23459771315993 0.195266300782888 6.32263584760924 2.57138788753732e-10 9.93198571561288e-07
ENSG00000196504 1798.5833532804 1.12058684067728 0.200292976167903 5.59473857804131 2.20954442255512e-08 5.68957688807944e-05
ENSG00000141425 1314.98437636948 0.969266426507209 0.176627372789244 5.4876342845441 4.07352327905737e-08 7.86699183267955e-05
ENSG00000164172 256.736460977222 1.27926938242781 0.240962413113634 5.3089997145095 1.1022850607365e-07 0.000170303041883789
ENSG00000173905 551.3081531682 1.17836225683544 0.223950630573659 5.26170546524926 1.42725274324417e-07 0.000183758790692687
> summary(res)
# 對res矩陣進行總結，利用summary命令統(tǒng)計顯示一共多少個genes上調和下調
out of 31349 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up) : 42, 0.13% #42個上調基因
LFC < 0 (down) : 8, 0.026% #8個下調基因
outliers [1] : 0, 0% #沒有離群值
low counts [2] : 23624, 75% #75%基因表達量不高俯萎，屬于low count
(mean count < 166)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results
> write.csv(res,file="All_results.csv") #將所有的數(shù)據(jù)輸出
> table(res$padj<0.05)#取padj小于0.05的數(shù)據(jù)傲宜，得到31行
FALSE TRUE
7694 31
image.png
提取差異基因分析的結果
在利用數(shù)據(jù)比較分析兩個樣品中同一個基因是否存在差異表達的時候，一般選取Foldchange值和經過FDR矯正過后的p值夫啊，取padj值(p值經過多重校驗校正后的值)小于0.05函卒，log2FoldChange大于1的基因作為差異基因集
> diff_gene_deseq2 <-subset(res,padj < 0.05 & (log2FoldChange > 1 | log2FoldChange < -1))
#也可以寫成diff_gene_deseq2 <-subset(res, padj < 0.05 & abs(log2FoldChange) > 1)
> dim(diff_gene_deseq2) #訪問數(shù)組
[1] 15 6
> head(diff_gene_deseq2)
log2 fold change (MLE): condition control vs treat
Wald test p-value: condition control vs treat
DataFrame with 6 rows and 6 columns
baseMean log2FoldChange lfcSE stat pvalue padj
<numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
ENSG00000178691 531.97208257559 2.89986543009505 0.234227480912105 12.3805516705499 3.33045617928295e-35 2.57277739849608e-31
ENSG00000135535 1195.85832961519 1.23459771315993 0.195266300782888 6.32263584760924 2.57138788753732e-10 9.93198571561288e-07
ENSG00000196504 1798.5833532804 1.12058684067728 0.200292976167903 5.59473857804131 2.20954442255512e-08 5.68957688807944e-05
ENSG00000164172 256.736460977222 1.27926938242781 0.240962413113634 5.3089997145095 1.1022850607365e-07 0.000170303041883789
ENSG00000173905 551.3081531682 1.17836225683544 0.223950630573659 5.26170546524926 1.42725274324417e-07 0.000183758790692687
ENSG00000172239 238.384608006302 1.30055924781765 0.252531313724155 5.15009100708309 2.60360094361487e-07 0.000287325961277498
> write.csv(diff_gene_deseq2,file= "DEG_treat_vs_control.csv")
#把所有的差異表達基因輸出一個文件
用bioMart對差異表達基因進行注釋
bioMart包是一個連接bioMart數(shù)據(jù)庫的R語言接口，能通過這個軟件包自由連接到bioMart數(shù)據(jù)庫撇眯。這個包可以做以下幾個工作：
1.查找某個基因在染色體上的位置报嵌。反之虱咧，給定染色體每一區(qū)間，返回該區(qū)間的基因s锚国；
2.通過EntrezGene的ID查找到相關序列的GO注釋腕巡。反之，給定相關的GO注釋血筑，獲取相關的EntrezGene的ID绘沉；
3.通過EntrezGene的ID查找到相關序列的上游100bp序列（可能包含啟動子等調控元件）；
4.查找人類染色體上每一段區(qū)域中已知的SNPs豺总；
5.給定一組的序列ID车伞，獲得其中具體的序列
> library('biomaRt')#載入R包，要用bioMart包從ensembl數(shù)據(jù)庫獲得基因的注釋
> library("curl")
> mart <- useDataset("hsapiens_gene_ensembl", useMart("ensembl"))
#用useMart函數(shù)選定數(shù)據(jù)庫喻喳，用useDataset()函數(shù)選定數(shù)據(jù)庫中的基因組.
#這條語句的意思是：選定ensembl數(shù)據(jù)庫中的hsapiens_gene_ensembl基因組
> my_ensembl_gene_id<-row.names(diff_gene_deseq2)
> hg_symbols<- getBM(attributes=c('ensembl_gene_id','external_gene_name',"description"),filters = 'ensembl_gene_id', values = my_ensembl_gene_id, mart = mart)
#用getBM()函數(shù)獲取注釋另玖。這個函數(shù)有4個參數(shù)：
#attributers()里面的值為我們要獲取的注釋類型，filters（）里面的值為我們已知的注釋類型
#values= 這個值就是我們已知的注釋類型的數(shù)據(jù)表伦，把上面我們通過數(shù)據(jù)處理得到的ensembl基因序號作為ensembl_gene_id 的值
#mart= 這個值是我們所選定的數(shù)據(jù)庫的基因組
#mart <- useDataset("hsapiens_gene_ensembl", useMart("ensembl"))
> head(hg_symbols)
ensembl_gene_id external_gene_name description
1 ENSG00000064666 CNN2 calponin 2 [Source:HGNC Symbol;Acc:HGNC:2156]
2 ENSG00000135535 CD164 CD164 molecule [Source:HGNC Symbol;Acc:HGNC:1632]
3 ENSG00000140526 ABHD2 abhydrolase domain containing 2 [Source:HGNC Symbol;Acc:HGNC:18717]
4 ENSG00000144357 UBR3 ubiquitin protein ligase E3 component n-recognin 3 [Source:HGNC Symbol;Acc:HGNC:30467]
5 ENSG00000152818 UTRN utrophin [Source:HGNC Symbol;Acc:HGNC:12635]
6 ENSG00000163848 ZNF148 zinc finger protein 148 [Source:HGNC Symbol;Acc:HGNC:12933]
> ensembl_gene_id<-rownames(diff_gene_deseq2)
> diff_gene_deseq2<-cbind(ensembl_gene_id,diff_gene_deseq2)
> colnames(diff_gene_deseq2)[1]<-c("ensembl_gene_id")
> diff_name<-merge(diff_gene_deseq2,hg_symbols,by="ensembl_gene_id")
#把注釋文件和基因表達量文件合并起來
> head(diff_name)
DataFrame with 6 rows and 9 columns
ensembl_gene_id baseMean log2FoldChange lfcSE stat pvalue padj
<character> <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
1 ENSG00000064666 233.649167762814 -1.05516841956132 0.261201368199475 -4.03967416723294 5.3525508533232e-05 0.0179775892790964
2 ENSG00000135535 1195.85832961519 1.23459771315993 0.195266300782888 6.32263584760924 2.57138788753732e-10 9.93198571561288e-07
3 ENSG00000140526 454.248659298761 1.09955045620162 0.244298072590034 4.50085604255593 6.76803336861755e-06 0.00326769111078566
4 ENSG00000144357 411.010386096953 1.00761195016462 0.244841355759736 4.11536665053185 3.86564433730833e-05 0.0149310512528534
5 ENSG00000152818 568.00627589957 1.22331171400082 0.254288749757859 4.81071897661886 1.50388340150822e-06 0.0011617499276651
6 ENSG00000163848 315.923451036876 1.12786190454491 0.232224503174498 4.85677389391337 1.19313704152794e-06 0.00102410929397815
external_gene_name description
<character> <character>
1 CNN2 calponin 2 [Source:HGNC Symbol;Acc:HGNC:2156]
2 CD164 CD164 molecule [Source:HGNC Symbol;Acc:HGNC:1632]
3 ABHD2 abhydrolase domain containing 2 [Source:HGNC Symbol;Acc:HGNC:18717]
4 UBR3 ubiquitin protein ligase E3 component n-recognin 3 [Source:HGNC Symbol;Acc:HGNC:30467]
5 UTRN utrophin [Source:HGNC Symbol;Acc:HGNC:12635]
6 ZNF148 zinc finger protein 148 [Source:HGNC Symbol;Acc:HGNC:12933]
查看CD164的情況：
> CD164 <- diff_name[diff_name$external_gene_name=="CD164",]
> CD164
DataFrame with 1 row and 9 columns
ensembl_gene_id baseMean log2FoldChange lfcSE stat pvalue padj external_gene_name
<character> <numeric> <numeric> <numeric> <numeric> <numeric> <numeric> <character>
1 ENSG00000135535 1195.85832961519 1.23459771315993 0.195266300782888 6.32263584760924 2.57138788753732e-10 9.93198571561288e-07 CD164
description
<character>
1 CD164 molecule [Source:HGNC Symbol;Acc:HGNC:1632]
數(shù)據(jù)可視化
1.MA plot
MA-plot 谦去，也叫 mean-difference plot或者Bland-Altman plot，用來估計模型中系數(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窃诉，衡量基因表達量變化，上調還是下調赤套。A表示每個基因的count的均值飘痛。根據(jù)summary可知，low count的比率很高容握，所以大部分基因表達量不高宣脉，也就是集中在0的附近（log2(1)=0，也就是變化1倍）.提供了模型預測系數(shù)的分布總覽剔氏。
library(ggplot2)
library(BiocGenerics)
plotMA(res,ylim=c(-2,2))
topGene <- rownames(res)[which.min(res$padj)]
with(res[topGene, ], {
points(baseMean, log2FoldChange, col="dodgerblue", cex=6, lwd=2)
text(baseMean, log2FoldChange, topGene, pos=2, col="dodgerblue")
})
image.png
lfcShrink 收縮log2 fold change
> res_order<-res[order(row.names(res)),]
#前面res結果已經按padj排序了塑猖，所以這次要按照行名升序再排列回來，否則和dds不一致
> res = res_order
> res.shrink <- lfcShrink(dds, contrast = c("condition","treat","control"), res=res)
#將DESeq函數(shù)處理后生成的的dds對象傳遞給lfcShrink函數(shù)即可
> plotMA(res.shrink, main="DESeq2", 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")
+ })
效果好的話谈跛，收縮的圖形不會出現(xiàn)左寬右窄羊苟。
image.png
火山圖 原始數(shù)據(jù)來源：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE50177
>BiocManager::install("EnhancedVolcano")
>library(EnhancedVolcano)
>library(airway)
>EnhancedVolcano(res,
lab = rownames(res),
x = 'log2FoldChange',
y = 'pvalue',
xlim = c(-8, 8),
title = 'resistant versus control',
pCutoff = 10e-17,
FCcutoff = 2,
transcriptPointSize = 1.5,
transcriptLabSize = 3.0,
col=c('black', 'blue', 'green', 'red1'),
colAlpha = 1,
legend=c('NS','Log2 fold-change','P value',
'P value & Log2 fold-change'),
legendPosition = 'right',
legendLabSize = 14,
legendIconSize = 5.0,
)
image.png
箱圖
> plotCounts(dds, gene="ENSG00000064666", intgroup="condition", returnData=TRUE) %>%
+ ggplot(aes(condition, count)) + geom_boxplot(aes(fill=condition)) + scale_y_log10() + ggtitle("CNN2")
image.png
PCA(principal components analysis)
熱圖(heatmap), PCA或聚類(clustering)我們需要data的轉換后的格式。
在DESeq2包中專門由一個PCA分析的函數(shù)感憾，即plotPCA蜡励。在DESeq2包中實際上已經有了歸一化的方法，rlog和vst。在使用的需要根據(jù)樣本量的多少來選擇方法凉倚。樣本量少于30的話兼都，選擇rlog，多于30的話稽寒，建議選擇vst扮碧。
#歸一化，因為樣本量超過了30瓦胎，因此用vst
> vsdata <- vst(dds, blind=FALSE)
> plotPCA(vsdata, intgroup="condition")
#intgroup:也就是在構建dds的時候的分組芬萍，實際上這個分組最后影響的是PCA圖中的顏色，但是并不影響PCA圖中各個樣本的位置搔啊。
#這個參數(shù)的作用就是當你的PCA圖需要加樣本名的時候柬祠，也就是你希望在PCA圖上知道哪個點是哪個樣本，你就需要導出這個负芋。
image.png
熱圖
1.count matrix 熱圖
library("pheatmap")
有兩種方式莽龟，畫出來的圖是一樣的毯盈。
第一種：
> select<-order(rowMeans(counts(dds_norm, normalized = TRUE)),decreasing = TRUE)[1:20]
> df <- as.data.frame(colData(dds_norm)[,c("condition","sizeFactor")])
> pheatmap(assay(ntd)[select,], cluster_rows=FALSE, show_rownames=FALSE,cluster_cols=FALSE, annotation_col=df)
第二種：
> ntd <- normTransform(dds_norm)
> pheatmap(assay(ntd)[select,], cluster_rows=FALSE, show_rownames=FALSE,
+ cluster_cols=FALSE, annotation_col=df)
畫出來的圖都是這樣：
image.png
sample-to-sample distances熱圖
轉換數(shù)據(jù)還可以做出樣本聚類熱圖。用dist函數(shù)來獲得sample-to-sample距離脑奠。距離矩陣熱圖中可以清楚看到samples之間的相似與否的總概宋欺。需要給heatmap函數(shù)基于sample距離提供等級聚類hc齿诞。
首先調用RColorBrewer包掌挚，這個包的說明請見：http://www.reibang.com/p/a8856757a0d2
> library("RColorBrewer")
> sampleDists <- dist(t(assay(vsdata)))
> sampleDistMatrix <- as.matrix(sampleDists)
> rownames(sampleDistMatrix) <- paste(vsdata$condition, vsdata$type, sep="-")
> colnames(sampleDistMatrix) <- NULL
> colors <- colorRampPalette( rev(brewer.pal(9, "Oranges")) )(255)
> pheatmap(sampleDistMatrix,
+ clustering_distance_rows=sampleDists,
+ clustering_distance_cols=sampleDists,
+ col=colors)
image.png
功能富集分析
A：差異基因富集分析（不需要表達值抽米，只需要gene name)
最常用的基因注釋工具是GO和KEGG注釋云茸，這基本上是差異基因分析一定做的兩件事懊纳。GO可以在GO:BP(生物過程）嗤疯，GO:MF(分子功能）茂缚，GO:CC(細胞組分)三個方面分別進行注釋脚囊，用的比較多的是GO:BP悔耘，但其他兩方面也很重要衬以。得到的差異表達基因列表就可以，也就是說不需要其他的值（http://www.reibang.com/p/5d82acad6008）
拿什么來富集备籽？ 得到的差異表達基因列表就可以车猬。
B: 基因集(gene set)富集分析（不管有無差異珠闰，需要全部genes表達值)
好處：可以發(fā)現(xiàn)被差異基因舍棄的genes可能參與了某重要生理過程或信號通路伏嗜。工具：GSEA
library(clusterProfiler)
library(DOSE)
library(org.Hs.eg.db)
setwd("/media/yanfang/FYWD/RNA_seq/matrix/")
sig.gene<-read.csv(file="readcount_all.csv")#用DEG_treat_vs_control.csv后續(xù)分析沒有富集到承绸，所以換用這個
head(sig.gene)
gene<-sig.gene[,2] #ID所在的列
head(gene)
gene.df<-bitr(gene, fromType = "ENSEMBL",
toType = c("SYMBOL","ENTREZID"),
OrgDb = org.Hs.eg.db)
# 人類樣品用org.Hs.eg.db军熏；小鼠用org.Mm.eg.db
head(gene.df)
library("stringr")
GO_CC富集
ego_cc <- enrichGO(gene = gene.df$ENSEMBL,
OrgDb = org.Hs.eg.db,
keyType = 'ENSEMBL',
ont = "CC",
pAdjustMethod = "BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05)
barplot(ego_cc,showCategory = 18,title="The GO_CC enrichment analysis of all DEGs ")+
scale_size(range=c(2, 12))+
scale_x_discrete(labels=function(ego_cc) str_wrap(ego_cc,width = 25))
GO_BP富集
ego_bp <- enrichGO(gene = gene.df$ENSEMBL,
OrgDb = org.Hs.eg.db,
keyType = 'ENSEMBL',
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05)
barplot(ego_bp,showCategory = 18,title="The GO_BP enrichment analysis of all DEGs ")+
scale_size(range=c(2, 12))+
scale_x_discrete(labels=function(ego_bp) str_wrap(ego_bp,width = 25))
image.png
image.png
#功能富集分析氣泡圖
#這個結果是另一篇文獻中的數(shù)據(jù)得來的
dotplot(ego_bp, font.size=10)
image.png
#GO圖
plotGOgraph(ego_bp)
image.png
KEGG富集分析
KEGG將基因組信息和高一級的功能信息有機地結合起來彤委，通過對細胞內已知生物學過程的計算機化處理和將現(xiàn)有的基因功能解釋標準化葫慎，對基因的功能進行系統(tǒng)化的分析偷办。
KEGG的另一個任務是一個將基因組中的一系列基因用一個細胞內的分子相互作用的網絡連接起來的過程椒涯，如一個通路或是一個復合物，通過它們來展現(xiàn)更高一級的生物學功能湖苞。
參考文章：http://blog.sciencenet.cn/blog-364884-779116.html
KEGG物種縮寫：http://www.genome.jp/kegg/catalog/org_list.html
GO和KEGG輸出文件解讀：http://www.bio-info-trainee.com/370.html
#這里我用另外的一篇文獻里的原始數(shù)據(jù)做的圖
>library(stringr)
>kk<-enrichKEGG(gene =gene.df$ENTREZID,
organism = 'hsa',
pvalueCutoff = 0.05)
>kk[1:30]
>barplot(kk,showCategory = 25, title="The KEGG enrichment analysis of all DEGs")+
scale_size(range=c(2, 12))+
scale_x_discrete(labels=function(kk) str_wrap(kk,width = 25))
>dotplot(kk,showCategory = 25, title="The KEGG enrichment analysis of all DEGs")+
scale_size(range=c(2, 12))+
scale_x_discrete(labels=function(kk) str_wrap(kk,width = 25))
image.png
Gene Set Enrichment Analysis（GSEA）
基因集富集分析 (Gene Set Enrichment Analysis, GSEA) 的基本思想是使用預定義的基因集（通常來自功能注釋或先前實驗的結果），將基因按照在兩類樣本中的差異表達程度排序，然后檢驗預先設定的基因集合是否在這個排序表的頂端或者底端富集捌臊±砼欤基因集合富集分析檢測基因集合而不是單個基因的表達變化矾端，因此可以包含這些細微的表達變化秩铆，預期得到更為理想的結果。（http://www.reibang.com/p/4910d7cec5c8）
參考資料：GSEA分析是什么鬼(上)和GSEA分析是什么鬼(下)添祸。
setwd("/media/yanfang/FYWD/RNA_seq/matrix/")
sig.gene <- read.csv(file="DEG_resistant_vs_control.csv")
genelist <- sig.gene$log2FoldChange
names(genelist) <- sig.gene[,1]
genelist <- sort(genelist, decreasing = TRUE)
gsemf <- gseGO(genelist, OrgDb = org.Hs.eg.db, keyType = "ENSEMBL", ont="BP")#OrgDb:指定物種注釋數(shù)據(jù)滚粟。keytype:ID類型
head(gsemf)
gseaplot(gsemf, geneSetID="GO:0030198")
image.png
