转录组（七）：提取差异基因

合并reads矩阵

#字符串不会转换成factor
options(stringsAsFactors = FALSE)

#读入四个文件
control1 <- read.table("~/rna_seq_AKAP95/data/rna_seq_test/matrix/SRR3589959.matrix",
                       sep="\t", col.names = c("gene_id","control1"))
control2 <- read.table("~/rna_seq_AKAP95/data/rna_seq_test/matrix/SRR3589961.matrix",
                       sep="\t", col.names = c("gene_id","control2"))
rep1 <- read.table("~/rna_seq_AKAP95/data/rna_seq_test/matrix/SRR3589960.matrix",
                       sep="\t", col.names = c("gene_id","apak951"))
rep2 <- read.table("~/rna_seq_AKAP95/data/rna_seq_test/matrix/SRR3589962.matrix",
                       sep="\t", col.names = c("gene_id","apak952"))

#合并四个文件
raw_count <- merge(merge(control1,control2,by = "gene_id"),
                   merge(rep1,rep2,by = "gene_id"))

#过滤前五行reads_count的结果
raw_count_filter <- raw_count[-1:-5,]

#处理gene_id小数,在EBI上无法查到带小数的gene_id
#从表面看两个反斜杠表示特殊含义
row.names(raw_count_filter) <- gsub("\\.\\d*","",
                                    raw_count_filter$gene_id)

#删除原来的rawnames
raw_count_filter <- raw_count_filter[,-1]

#查看AKAP95表达量
AKAP95 <- raw_count_filter[rownames(raw_count_filter)=="ENSMUSG00000024045",]

构建dds对象

#分组信息，为因子类型
condition <- factor(c(rep("control",2),rep("AKAP95",2)),
                    levels = c("control","AKAP95"))

#与raw_count_filter一样
countdata <- raw_count_filter[,1:4]

colData <- data.frame(row.names = colnames(raw_count_filter), 
                      condition)
dds <- DESeqDataSetFromMatrix(countdata, colData, design= ~ condition)
## design 其实也是一个对象，还可以通过design(dds)来继续修改分组信息，但是一般没有必要。

normalization——消除量纲对数据结构的影响

#标准化一步即可
dds2 <- DESeq(dds)

#查看名称
resultsNames(dds2)

#获取结果
res <- results(dds2)

#查看概要信息
summary(res)

由于只选取了125000条的reads，所以在默认adj_p_value < 0.1的情况下，没有上调和下调的基因。

提取差异基因

#取padj(由p_value校正后得到)小于0.05，log2FoldChange的绝对值大于1作为帅选差异基因的条件
#但是少量reads得到的结果可信度不高
table(res$padj<0.05)
#FALSE 
#13397

#所以为了下游分析的进行，我调整了筛选条件，但这样就失去了生物学意义。
table(res$pvalue<0.1 & (res$log2FoldChange > 1 | res$log2FoldChange < -1))
#FALSE  TRUE 
#13289   108 

#按pvalue从小到大排列，这一步不拍也没有事
res <- res[order(res$pvalue),]
#参考[https://blog.csdn.net/u012543538/article/details/16340907](https://blog.csdn.net/u012543538/article/details/16340907)

#取子集
diff_gene_deseq2 <-subset(res,pvalue < 0.1 & (log2FoldChange > 1 | log2FoldChange < -1))
diff_gene <- as.data.frame(diff_gene_deseq2,sort=FALSE)

参考
http://www.biotrainee.com/thread-1748-1-1.html

https://www.jianshu.com/p/3bfb21d24b74