转录组练习（6）
原贴地址：https://www.jianshu.com/p/3bfb21d24b74
https://www.jianshu.com/p/4910d7cec5c8
1.读取自己表达矩阵
# 构建自己的表达矩阵并读取
> control1 <- read.table("~/disk2/data/rna-seq/matrix/SRR3589959.count", sep="\t", col.names = c("gene_id","control1"))
> control2 <- read.table("~/disk2/data/rna-seq/matrix/SRR3589961.count", sep="\t", col.names = c("gene_id","control2"))
> rep1 <- read.table("~/disk2/data/rna-seq/matrix/SRR3589960.count", sep="\t", col.names = c("gene_id","akap951"))
> rep2 <- read.table("~/disk2/data/rna-seq/matrix/SRR3589962.count", sep="\t",col.names = c("gene_id","akap952"))
> raw_count <- merge(merge(control1, control2,by="gene_id"),merge(rep1,rep2, by="gene_id"))
> raw_count_filt <- raw_count[-48823:-48825,]
> raw_count_filter <- raw_count_filt[-1:-2,]
> ENSEMBL <- gsub("\\.\\d*","", raw_count_filter$gene_id)
> row.names(raw_count_filter) <- ENSEMBL
> raw_count_filter <- raw_count_filter[ ,-1]
这里有两个问题，第一个就是删除行列，自己根据情况选择，理论上把名字差别太大的删了。
矩阵数据结构
2.构建dds对象
# 这一步很关键，要明白condition这里是因子，不是样本名称；小鼠数据有对照组和处理组，各两个重复
> condition <- factor(c(rep("control",2),rep("akap95",2)), levels = c("control","akap95"))# 获取count数据
> countData <- raw_count_filter[,1:4]
> colData <- data.frame(row.names=colnames(raw_count_filter), condition)
> dds <- DESeqDataSetFromMatrix(countData, colData, design= ~ condition)# 查看一下dds的内容
> head(dds)
adds概要信息
3.DESeq标准化dds
# normalize 数据
> dds2 <- DESeq(dds)# 查看结果的名称，本次实验中是"Intercept"，"condition_akap95_vs_control"
> resultsNames(dds2)# 将结果用results()函数来获取，赋值给res变量
res <- results(dds2)# summary一下，看一下结果的概要信息summary(res)
result结果可以看到一些基本的信息，p值默认小于0.1，上调基因有625个，下调基因有445个。
res的概要信息
4.提取差异分析结果
# 获取padj（p值经过多重校验校正后的值）小于0.05，表达倍数取以2为对数后大于1或者小于-1的差异表达基因。
> table(res$padj<0.05)
> res <- res[order(res$padj),]
> diff_gene_deseq2 <-subset(res,padj <0.05& (log2FoldChange >1| log2FoldChange <-1))
> diff_gene_deseq2 <- row.names(diff_gene_deseq2)
resdata<-merge(as.data.frame(res),as.data.frame(counts(dds2,normalize=TRUE)),by="row.names",sort=FALSE)
# 得到csv格式的差异表达分析结果
> write.csv(resdata,file="control_vs_akap95.cvs",row.names = F)
把名字变成上图，gsub那个不会整，那就导出csv，然后手工删除在导入进行后续分析
后续就可以转换基因名 各种折腾了
require(DOSE)
require(clusterProfiler)
ekk <- enrichKEGG(gene=gene,organism="human",pvalueCutoff=0.01)
ego <- enrichGO(gene=gene,OrgDb="org.Hs.eg.db",ont="CC",pvalueCutoff=0.01,readable=TRUE)
write.csv(summary(ekk),"KEGG-enrich.csv",row.names =F)
write.csv(summary(ego),"GO-enrich.csv",row.names =F)
ego <- enrichGO(gene = gene,
universe = names(geneList),
OrgDb = org.Hs.eg.db,
ont = "CC",
pAdjustMethod = "BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05)
ego2 <- enrichGO(gene = gene.df$ENSEMBL,
OrgDb = org.Hs.eg.db,
keytype = 'ENSEMBL',
ont = "CC",
pAdjustMethod = "BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05)
ego3 <- enrichGO(gene = gene.df$SYMBOL,
OrgDb = “org.Hs.eg.db”,
keytype = 'SYMBOL',
ont = "CC",
pAdjustMethod = "BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05)
ego <- enrichGO(gene = DEG$,
OrgDb = "org.Hs.eg.db",
ont = "CC",
pAdjustMethod = "BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05)
ego <- enrichGO(gene=gene,OrgDb="org.Hs.eg.db",ont="CC",pvalueCutoff=0.01,readable=TRUE)
DO<-enrichDO(gene=DEG$Gene.ID, ont = "DO", pvalueCutoff = 0.01, pAdjustMethod = "BH",qvalueCutoff = 0.05)
eg=bitr(geneID = PP$geneID, "ENTREZID", "SYMBOL", "org.Hs.eg.db")
df.id<-bitr(df$SYMBOL, fromType ="SYMBOL", toType ="ENTREZID",OrgDb ="org.Hs.eg.db")
easy.df<-merge(df,df.id,by="SYMBOL",all=F)
sortdf<-easy.df[order(easy.df$foldChange, decreasing =T),]
gene.expr = sortdf$foldChange
names(gene.expr) <- sortdf$ENTREZID
