转录组练习(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

推荐阅读更多精彩内容

  • 一、读文章获取下载数据 1、读文章 一般我都从NCBI上面下载文章,找到数据号 2、下载数据 进入NCBI的GEO...
    黄思源_3a22阅读 4,268评论 0 0
  • pyspark.sql模块 模块上下文 Spark SQL和DataFrames的重要类: pyspark.sql...
    mpro阅读 7,942评论 0 11
  • library(Biobase)library(GEOquery)library(limma)## load se...
    苏慕晨枫阅读 1,126评论 1 1
  • mean to add the formatted="false" attribute?.[ 46% 47325/...
    ProZoom阅读 1,465评论 0 3
  • 打开终端 sudo vi ~/.bash_profile 然后黏贴 export GOPATH=$HOME/wo...
    炳良哥哥阅读 379评论 0 49