DESeq2 - lmigueel/Bioinformatica GitHub Wiki

1. Sobre

DESeq2 é um pacote para expressão diferencial de genes/transcritos. Ele estima a dependência média da variância dos read counts dos sequenciamentos NGS, e teste a expressão diferencial com base em um modelo usando a distribuição binomial negativa.

2. Script

Vamos utilizar o R para realizar a expressão diferencial dos genes.

if (! require(DESeq2)) {
  source("https://bioconductor.org/biocLite.R")
  biocLite("DESeq2")
}

library(DESeq2)

De modo geral temos que passar as matrizes dos COUNTS para o DESeq2. Essa matriz pode ser gerada através do kallisto, salmon, HTSeq count ou RSEM.

data = read.table("kallisto.isoform.counts.matrix", header=T, row.names=1, com='')

Todos os comandos abaixo são para reordenar as colunas da matriz, caso alguma esteja invertida, além de deixar os dados como inteiros (round()) e corta para no mínimo 2 counts por linhas.

col_ordering = c(1,2,3,4,5,6)
rnaseqMatrix = data[,col_ordering]
rnaseqMatrix = round(rnaseqMatrix)
rnaseqMatrix = rnaseqMatrix[rowSums(rnaseqMatrix) >= 2,]

O próximo passo é alterar as colunas e criar o Full Dataframe, que contém todas as informações.

conditions = data.frame(conditions=factor(c(rep("Controle", 3), rep("Tratado", 3))))
rownames(conditions) = colnames(rnaseqMatrix)
ddsFullCountTable <- DESeqDataSetFromMatrix(
  countData = rnaseqMatrix,
  colData = conditions,
  design = ~ conditions)

Agora partimos para análise em si, e geramos o resultado. Se você colocar a variável contrast sendo como indicado abaixo, o valor do log2FC, quando positivo, indicará uma maior expressão no Tratado em relação ao controle.

dds = DESeq(ddsFullCountTable)
contrast=c("conditions","Tratado","Controle")
res = results(dds, contrast)

Melhorando o resultado

baseMeanA <- rowMeans(counts(dds, normalized=TRUE)[,colData(dds)$conditions == "Controle"])
baseMeanB <- rowMeans(counts(dds, normalized=TRUE)[,colData(dds)$conditions == "Tratado"])
res = cbind(baseMeanA, baseMeanB, as.data.frame(res))
res = cbind(sampleA="Controle", sampleB="Tratado", as.data.frame(res))
res$padj[is.na(res$padj)]  <- 1
res = as.data.frame(res[order(res$pvalue),])

Salvando os resultados ...

write.table(res, file='DESeq2.DE_results', sep='\t', quote=FALSE)

Salvando o Volcano plot ...

pdf("DESeq2_Volcano.pdf")
plot(res$log2FoldChange, -1*log10(res$padj), col=ifelse(res$padj<=0.05, "red", "black"),xlab="logCounts", ylab="logFC", title="Volcano plot", pch=20)
dev.off()

Lista de diferenciais ...

Deseq2_list <- NULL
for(i in 1:length(res[,1])){
  if(res[i,]$padj <= 0.05){
    #print("Entrou\n")
    Deseq2_list[i] <- row.names(res)[i]
  }
}

Script geral

if (! require(DESeq2)) {
  source("https://bioconductor.org/biocLite.R")
  biocLite("DESeq2")
  library(DESeq2)
}

data = read.table("kallisto.isoform.counts.matrix", header=T, row.names=1, com='')
col_ordering = c(1,2,3,4,5,6)
rnaseqMatrix = data[,col_ordering]
rnaseqMatrix = round(rnaseqMatrix)
rnaseqMatrix = rnaseqMatrix[rowSums(rnaseqMatrix) >= 2,]
conditions = data.frame(conditions=factor(c(rep("Controle", 3), rep("Tratado", 3))))
rownames(conditions) = colnames(rnaseqMatrix)
ddsFullCountTable <- DESeqDataSetFromMatrix(
  countData = rnaseqMatrix,
  colData = conditions,
  design = ~ conditions)
dds = DESeq(ddsFullCountTable)
contrast=c("conditions","Tratado","Controle")
res = results(dds, contrast)
baseMeanA <- rowMeans(counts(dds, normalized=TRUE)[,colData(dds)$conditions == "Controle"])
baseMeanB <- rowMeans(counts(dds, normalized=TRUE)[,colData(dds)$conditions == "Tratado"])
res = cbind(baseMeanA, baseMeanB, as.data.frame(res))
res = cbind(sampleA="Controle", sampleB="Tratado", as.data.frame(res))
res$padj[is.na(res$padj)]  <- 1
res = as.data.frame(res[order(res$pvalue),])
write.table(res, file='DESeq2.DE_results', sep='     ', quote=FALSE)

pdf("DESeq2_Volcano.pdf")
plot(res$log2FoldChange, -1*log10(res$padj), col=ifelse(res$padj<=0.05, "red", "black"),xlab="logCounts", ylab="logFC", title="Volcano plot", pch=20)
dev.off()

Deseq2_list <- NULL
for(i in 1:length(res[,1])){
  if(res[i,]$padj <= 0.05){
    #print("Entrou\n")
    Deseq2_list[i] <- row.names(res)[i]
  }
}