DESeq2에서 heatmap, PCA, MA, volcano plot 그리기

DESeq2에서 heatmap, PCA, MA, volcano plot 그리기

분석을 위한 DataSet이 필요하다. 예시에서는 여러 cell line에 몇 가지 drug를 각기 다른 dose로 처리했을 때의 결과를 분석하고자 한다.

DEG분석은 dose에만 초점이 맞춰지도록 디자인 하였다.

res변수는 결과가 adjust p-value 낮은순서대로 sorting하여 저장하였다.

sampleTable<-data.frame(sampleName=sampleFiles, fileName=sampleFiles, dose=I(dose), drug=drug, cell_line=cell_line)

ddsHTSeq<-DESeqDataSetFromHTSeqCount(sampleTable=sampleTable, directory=directory, design=~dose)

dds<-DESeq(ddsHTSeq)

res <- results(dds)

res <- res[order(res$padj),]

heatmap을 그리기 위한 color를 지정해 준다. 지정하지 않으면 그림을 그릴 때 마다 다른 color가 나오기 때문에 결과를 정리하는데 매 번 수정하기가 상당히 귀찮기 때문이다. 각 항목에 맞는 값을 넣어주고 뒤에 16진법 RGB코드를 넣어주면 된다. 

16진법 RGB코드는 아래 포스팅한 방법으로 찾으면 쉽게 찾을 수 있다.

2018/08/31 - [etc.] - 16진수 RGB코드 알아내는법

## annotation color

ann_colors = list(

        dose = c(".." = "#F9EBD1"),

        cell_line = c(".." = "#FFFFFF"),

        drug = c(".." = "#BDBE32")

)

heatmap을 그리기 위해서 pheatmap 라이브러리를 로딩한다. order뒤에 [1:1000]이라고 되어 있는 부분은 유전자가 너무 많아서 발현량 기준으로 상위 1000개만 가져온 것이다. 전부 넣으려다가는 메모리 덤프가 생길 것이다. 1000개라는 갯수는 임의로 설정한 것이다.

## heatmap

library("pheatmap")

ntd <- normTransform(dds)

select <- order(rowMeans(counts(dds,normalized=TRUE)),

                decreasing=TRUE)[1:1000]

df <- as.data.frame(colData(dds)[,c("cell_line","drug","dose")])

rownames(df) = sampleFiles

pdf("heatmap_normalized_readcount.pdf")

pheatmap(assay(ntd)[select,], cluster_rows=T, show_rownames=F,

         cluster_cols=T, annotation_col=df, annotation_colors = ann_colors)

dev.off()

plotPCA는 함수 자체를 다시 정의한 것인데 이유는 DESeq2에서 제공하는 plotPCA함수는 PC1과 PC2밖에 고려하지 않기 때문이다.

마지막 ggplot을 불러오는 함수에서 x와 y를 PC3, PC4로 바꾸면 결과가 달라진다. PCA 분석에서 항상 PC1과 PC2의 조합이 정답이 아니기 때문에 필요한 방법이다.

## PCA

plotPCA <- function (object, intgroup = "condition", ntop = 500, returnData = FALSE)

{

    rv <- rowVars(assay(object))

    select <- order(rv, decreasing = TRUE)[seq_len(min(ntop,

        length(rv)))]

    pca <- prcomp(t(assay(object)[select, ]))

    percentVar <- pca$sdev^2/sum(pca$sdev^2)

    if (!all(intgroup %in% names(colData(object)))) {

        stop("the argument 'intgroup' should specify columns of colData(dds)")

    }

    intgroup.df <- as.data.frame(colData(object)[, intgroup,

        drop = FALSE])

    group <- if (length(intgroup) > 1) {

        factor(apply(intgroup.df, 1, paste, collapse = " : "))

    }

    else {

        colData(object)[[intgroup]]

    }

    d <- data.frame(PC1 = pca$x[, 1], PC2 = pca$x[, 2], PC3 = pca$x[, 3], PC4 = pca$x[, 4], group = group,

        intgroup.df, name = colnames(object))

    if (returnData) {

        attr(d, "percentVar") <- percentVar[1:2]

        return(d)

    }

    ggplot(data = d, aes_string(x = "PC1", y = "PC2", color = "group")) +

        geom_point(size = 3) + xlab(paste0("PC1: ", round(percentVar[1] *

        100), "% variance")) + ylab(paste0("PC2: ", round(percentVar[2] *

        100), "% variance")) + coord_fixed()

}

se <- SummarizedExperiment(log2(counts(dds, normalized=TRUE) + 1), colData=colData(dds,levels=c('0','0.1','1')))

pdf("PCA_normlaized_readcount.pdf")

plotPCA(DESeqTransform(se),intgroup=c("cell_line","drug","dose"),ntop=500)

dev.off()

MAplot은 유전자의 발현 값과 p-value의 상관관계를 보기 위해서 사용된다. 

#MAplot

pdf("MA_plot.pdf")

plotMA(res, ylim=c(-2,2))

dev.off()

Volcano plot은 각 유전자의 log2FoldChange와 p-value를 보기 위해 사용된다. plotting을 하는건 위의 4번째 줄 까지만 입력하면 되지만 특정 조건에 맞는 유전자일 때 색을 다르게 주거나 cutoff을 선으로 표시하기 위한 줄이 추가되어 있다.

#volcano plot

pdf("volcano_plot.pdf")

par(mar=c(5,5,5,5), cex=1.0, cex.main=1.4, cex.axis=1.4, cex.lab=1.4)

topT <- as.data.frame(res)

with(topT, plot(log2FoldChange, -log10(padj), pch=20, main="Volcano plot", cex=1.0, xlab=bquote(~Log[2]~fold~change), ylab=bquote(~-log[10]~Q~value)))

with(subset(topT, padj<0.05 & abs(log2FoldChange)>2), points(log2FoldChange, -log10(padj), pch=20, col="red", cex=0.5))

#Add lines for absolute FC>2 and P-value cut-off at FDR Q<0.05

abline(v=0, col="black", lty=3, lwd=1.0)

abline(v=-2, col="black", lty=4, lwd=2.0)

abline(v=2, col="black", lty=4, lwd=2.0)

abline(h=-log10(max(topT$pvalue[topT$padj<0.05], na.rm=TRUE)), col="black", lty=4, lwd=2.0)

dev.off()