作业要求
这个步骤推荐在R里面做，载入表达矩阵，然后设置好分组信息，统一用DEseq2进行差异分析，当然也可以走走edgeR或者limma的voom流程。
基本任务是得到差异分析结果，进阶任务是比较多个差异分析结果的异同点。
来源于生信技能树：http://www.biotrainee.com/forum.php?mod=viewthread&tid=1750#lastpost
实验过程
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)
resdata数据结构
有大神们的帮助，总算是完成了这一课的内容，感谢黯蓝小伙伴的帮助，还有群里小伙伴的帮助，当然还参考了Jimmy大神的博客，此外还参考了张翼翔的贴文：http://www.biotrainee.com/thread-1984-1-2.html。继续下一课的内容，最后一课的内容，GO富集分析。
