作業(yè)要求
這個步驟推薦在R里面做，載入表達矩陣某饰，然后設置好分組信息，統(tǒng)一用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渡紫。繼續(xù)下一課的內容，最后一課的內容考赛，GO富集分析惕澎。
