DESeq2

DESeq2에서 heatmap, PCA, MA, volcano plot 그리기 2018.08.31
DESeq2에서 연속적인 값을 condtion으로 받기 2018.08.10
DESeq2 에서 multiple condition 수행하기 2018.07.27 1

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

바닐라스카이 2018. 8. 31. 16:48

2018. 8. 31. 16:48

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()

저작자표시 비영리 변경금지

'bioinformatics' 카테고리의 다른 글

oncotator 설치 및 실행하기 (0)	2018.10.04
liftover하기 (0)	2018.09.28
Optical duplicate와 Library duplicates (0)	2018.08.27
Genome에 존재하는 Variant란? (0)	2018.08.20
DESeq2에서 연속적인 값을 condtion으로 받기 (0)	2018.08.10

DESeq2에서 연속적인 값을 condtion으로 받기

바닐라스카이 2018. 8. 10. 09:03

2018. 8. 10. 09:03

DESeq2에서 연속적인 값을 condtion으로 받기

DESeq2에서는 2가지 이상의 condition에 대한 차이를 보고 Differentially Expressed Genes를 구할 수도 있지만 시간이나 약물 투여 농도에 따른 DEGs도 계산할 수 있다.

아래의 예시처럼 약물 농도는 0nM, 0.1nM, 1nM을 투여하고 이에 따른 유전자의 변화를 가지고 p-value를 구하게 된다.

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

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

dose=I(dose) 이 부분에서 I가 연속적인 value를 의미한다. 물론 dose안의 모든 값은 숫자로만 저장되어야 하며 dose에 기반하여 DEG를 구하여야 하기 때문에 밑에design=~dose가 되는 것이다.

design에 여러가지 factor를 넣어서 같이 계산할 수도 있지만 실험설계를 그렇게 해본 적은 없기 때문에 아래 참조 페이지에서 직접 확인하길 바란다.

그림이 너무 작은데 계속 진행하면 아래와 같은 결과를 얻을 수 있다.

내용은 임의로 채워놓았기때문에 크게 신경쓰지 말고 보면 마지막 6개의 column이 의미하는 것은 농도가 0, 1, 3으로 올라갔을 때 normalized readcount가 급격히 감소하는 것을 볼 수 있으며 log2FoldChange도 -값으로 크고, 당연히 결과적으로 adjust p-value도 낮은것을 확인할 수 있다.

Reference -

https://support.bioconductor.org/p/97262/

저작자표시 비영리 변경금지

'bioinformatics' 카테고리의 다른 글

Optical duplicate와 Library duplicates (0)	2018.08.27
Genome에 존재하는 Variant란? (0)	2018.08.20
RSeQC를 사용하여 stranded 데이터 확인하기 (0)	2018.08.06
RNA-seq 라이브러리 종류와 구별법 (1)	2018.08.03
MultiQC 설치 및 실행하기 (1)	2018.07.31

DESeq2 에서 multiple condition 수행하기

바닐라스카이 2018. 7. 27. 14:15

2018. 7. 27. 14:15

DESeq2 에서 multiple condition 수행하기

DESeq2는 여러 컨디션이 있어도 DEG분석이란 결국 pairwise하게 비교해야 하기 때문에 condtion을 구분 해 놓은뒤 결과를 원하는 두 컨디션을 놓고 비교해서 각각 분석해야 한다.

SRA010153 데이터로 전체적인 DEG analysis workflow는 추후에 정리하기로 하고 여기서는 multiple condition의 비교에만 중점을 두고 설명하겠다.

SRA010153은 여러 컨디션이 있지만 Brain과 UHR phi X 데이터만 받은 상태에서 임의로 test 컨디션을 아래처럼 추가했다. (원래는 7개의 Brain과 7개의 UHR이다.)

sampleCondition<-c('Brain','Brain','Brain','Brain','Brain','Brain','test', 'UHR','UHR','UHR','UHR','test','UHR','test')

sampleTable<-data.frame(sampleName=sampleFiles, fileName=sampleFiles, condition=sampleCondition)

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

write.csv(assay(ddsHTSeq),quote=F,file="raw_count.csv")

colData(ddsHTSeq)$condition<-factor(colData(ddsHTSeq)$condition, levels=c('Brain','UHR','test'))

cts <- counts(ddsHTSeq)

geoMeans <- apply(cts, 1, function(row) if (all(row == 0)) 0 else exp(mean(log(row[row != 0]))))

ddsHTSeq <- estimateSizeFactors(ddsHTSeq, geoMeans=geoMeans)

dds<-DESeq(ddsHTSeq)

res <- results(dds)

쭉 진행하면 나오는 결과는 test와 brain간의 비교이다. DESeq2는 default로 가장 처음과 마지막에 있는 컨디션을 비교하기 때문이다.

이를 아래처럼 수정하면 test와 UHR간의 비교 결과를 보여준다.

results(dds,contrast=c('condition','test','UHR'))

2, 3번째 값을 수정하면 그에 맞게 바뀌는데 당연히 이미 설정해놓은 컨디션 값을 지정해주어야하며

순서를 바꾸면 log2FoldChange값이 +에서 -로 바뀌게 된다.

> results(dds,contrast=c('condition','UHR','test'))

log2 fold change (MLE): condition UHR vs test

Wald test p-value: condition UHR vs test

DataFrame with 60461 rows and 6 columns

baseMean log2FoldChange lfcSE

ENSG00000000003 92.1871316897525 0.479423175690374 0.43657650183317

ENSG00000000005 2.33464385903914 0.243364522893621 0.807778836598724

> results(dds,contrast=c('condition','test','UHR'))

log2 fold change (MLE): condition test vs UHR

Wald test p-value: condition test vs UHR

DataFrame with 60461 rows and 6 columns

baseMean log2FoldChange lfcSE

ENSG00000000003 92.1871316897525 -0.479423175690374 0.43657650183317

ENSG00000000005 2.33464385903914 -0.243364522893621 0.807778836598724

log2FC값 외에 변경되는 값은 없다.

저작자표시 비영리 변경금지

'Computer Science > R' 카테고리의 다른 글

Kegg pathway에 속하는 유전자 정보 가져오기 (2)	2018.09.11
pheatmap으로 heatmap그리기 (0)	2018.09.11
R에서 Dataframe 합치기 (0)	2018.09.05
centos 6에서 R 설치를 위한 라이브러리 설치 (2)	2018.07.25
Arguments in R (0)	2018.07.25

PREV 이전 1 NEXT 다음

Be great

DESeq2