RNA-seq之差异表达分析： 使用DESeq2进行差异基因归类

RNA-seq之差异表达分析： 使用DESeq2进行差异基因归类

这是本人进行课题的时候一个实例，本篇博客仅用于个人的代码记录：
希望对有缘点进来的你有一些作用。

1. 安装过程所需要的包

使用Biomanager包安装Bioconductor里面的DESeq2包

BiocManager::install("DESeq2")

如果过程中遇到某些包的版本不对，将其手动删除即可。

2.读入数据-构建dds矩阵-输出

getwd()
library("DESeq2")
countData <- as.matrix(read.csv('gene_count.DH_DXH.csv',header=T,row.names=1,sep=","))
condition<-factor(c("DH","DH","DH","DXH","DXH","DXH"))
coldata<-data.frame(row.names=colnames(countData),condition)  
dds<-DESeqDataSetFromMatrix(countData,coldata,design=~condition)
dds <- DESeq(dds)
res <- results(dds) 
out <- cbind(res$log2FoldChange,res$pvalue,res$padj) 
rownames(out) <- row.names(res) 
colnames(out) <- c('log2(FoldChange)','Pvalue','FDR') 
FCcut <- 2 
FDRcut <- 0.05 
diff <- out[(!is.na(res$padj) & res$padj<FDRcut)&abs(res$log2FoldChange)>abs(log2(FCcut)),]
write.table(out,"DH_DXH-DESeq2_AllResult.table", sep = "\t", quote=F,row.names = T,col.names=T) 
write.table(diff,"DH_DXH-DESeq2_DEG.table", sep = "\t", row.names = TRUE,quote=F) 
write.table(rownames(diff),"DH_DXH-DESeq2_DEGlist.txt", sep = "\t",quote=F,row.names=F,col.names=F)

PS：dds是表达矩阵
coldata是用来指定样本分组的数据集

为了让自己的那么多组对比迅速我写了个破perl脚本，
你问我为什么不用R自带的循环做，因为我R语言循环写的烂啊。
反正perl脚本只需要改几个参数就ok了不是吗？

奇怪的perl脚本

use strict;
use warnings;

my @files=`ls ./DEapp/*.csv`;

for my $file(@files){
    chomp($file);
    #print "$file";exit;
    $file =~ /(.*)gene_count\.(\w+)_(\w+)\.csv/;
    #print "$2\t$3\n";exit;
    print "
countData <- as.matrix(read.csv(\'gene_count.$2_$3.csv\',header=T,row.names=1,sep=\",\"))
condition<-factor(c(\"$2\",\"$2\",\"$2\",\"$3\",\"$3\",\"$3\"))
coldata<-data.frame(row.names=colnames(countData),condition)
dds<-DESeqDataSetFromMatrix(countData,coldata,design=~condition)
dds <- DESeq(dds)
res <- results(dds)
out <- cbind(res\$log2FoldChange,res\$pvalue,res\$padj)
rownames(out) <- row.names(res)
colnames(out) <- c('log2(FoldChange)','Pvalue','FDR')
FCcut <- 2
FDRcut <- 0.05
diff <- out[(!is.na(res\$padj) & res\$padj<FDRcut)&abs(res\$log2FoldChange)>abs(log2(FCcut)),]
write.table(out,\"$2_$3-DESeq2_AllResult.table\", sep = \"\\t\", quote=F,row.names = T,col.names=T)
write.table(diff,\"$2_$3-DESeq2_DEG.table\", sep = \"\\t\", row.names = TRUE,quote=F)
write.table(rownames(diff),\"$2_$3-DESeq2_DEGlist.txt\", sep = \"\\t\",quote=F,row.names=F,col.names=F)
"
}

至于后续的火山图和venn图等我什么时候想做了就去弄吧，现在确实有点懒得弄了。
又是个批量循环操作。啊这，有点难受的。