【생물정보학】 bulk RNA-seq에서 DEG 획득하기 (DESeq2, t test, ANOVA)

bult RNA-seq에서 DEG 획득하기 (DESeq2, t test, ANOVA)

 

추천글 : 【생물정보학】 생물정보학 분석 목차 

---

1. 개요 [본문]

2. 함수 [본문]

3. 예제 [본문]

---

a. 전사체 분석 파이프라인 

---

 

1. 개요 [목차]

⑴ 두 조건 DEG 분석에서 Wilcoxon rank-sum test가 DESeq2, edgeR, limma-voom, dearseq, NOISeq보다 더 좋은 성능을 보임

⑵ DESeq2, edgeR가 부정확할 수 있는 이유

① 유전자 발현이 음이항 분포를 따른다는 가정은 homogeneous distribution에서만 쓸 수 있음

② expression level은 음이항 분포를 따르고 measurement는 Poisson 분포를 따르기에 실제로는 복잡한 양상으로 나타남

 

 

2. 함수 [목차]

⑴ 2개 집단의 DESeq2 기반 DEG 획득하기 : R 기반

 

a = 'DSS_ST_late'
b = 'DSS_Saline_late'

Idents(pbmc_gene) <- pbmc_gene$orig.ident
pbmc_gene_ = pbmc_gene[, pbmc_gene$orig.ident %in% c(a, b)]
counts <- pbmc_gene_@assays$RNA@counts
sample_info <- data.frame(
  row.names = colnames(counts),
  condition = pbmc_gene_$orig.ident[colnames(counts)]
)
dds <- DESeqDataSetFromMatrix(
  countData = counts,
  colData = sample_info,
  design = ~ condition
)
dds <- DESeq(dds)
res <- results(dds, contrast=c("condition", a, b))
res_ordered <- res[order(res$log2FoldChange), ]
res_significant <- subset(res_ordered, padj < 0.05) 
res_significant

 

⑵ 여러 집단의 DESeq2 기반 DEG 획득하기

 

library(compGenomRData)
library(DESeq2)
library(stats)

counts_file  <- system.file ("extdata/rna-seq/SRP029880.raw_counts.tsv", package = "compGenomRData")
coldata_file  <- system.file ("extdata/rna-seq/SRP029880.colData.tsv", package = "compGenomRData")
counts <- as.matrix(read.table(counts_file, header = T, sep = '\t'))
countData <- as.matrix(subset(counts, select = c(-width)))
colData <- read.table(coldata_file, header = T, sep = '\t', 
                      stringsAsFactors = TRUE)
designFormula <- "~ group"
dds <- DESeqDataSetFromMatrix(countData = countData, 
                              colData = colData, 
                              design = as.formula(designFormula))
dds <- dds[ rowSums(DESeq2::counts(dds)) > 1, ]
dds <- DESeq(dds)
DEresults = results(dds, contrast = c("group", 'CASE', 'CTRL'))
DEresults <- DEresults[order(DEresults$pvalue),]

 

⑶ 두 집단의 t test

 

t.test(v1, v2, paired = FALSE)
# maybe you can activate 'paired' in a special condition

 

⑷ 두 집단의 ANOVA test

 

one_way_2_factor_anova <- function(v1, v2){
  dat <- matrix(0, nrow = ( length(v1) + length(v2) ), ncol = 2 )
  for(i in 1 : length(v1) ){
    dat[i, 1] <- v1[i]
    dat[i, 2] <- 'v1'
  }
  for(i in 1 : length(v2) ){
    dat[i + length(v1), 1] <- v2[i]
    dat[i + length(v1), 2] <- 'v2'
  }
  dat <- as.data.frame(dat)
  
  colnames(dat) <- c('val', 'factor')

  anova_IS <- aov(val ~ factor, data = dat)
  print(summary(anova_IS))

  anova_residuals <- anova_IS$residuals
  print(summary(anova_residuals))
}

 

⑸ 세 집단의 ANOVA test 

 

one_way_3_factor_anova <- function(v1, v2, v3){
  dat <- matrix(0, nrow = ( length(v1) + length(v2) + length(v3) ), ncol = 2 )
  for(i in 1 : length(v1) ){
    dat[i, 1] <- v1[i]
    dat[i, 2] <- 'v1'
  }
  for(i in 1 : length(v2) ){
    dat[i + length(v1), 1] <- v2[i]
    dat[i + length(v1), 2] <- 'v2'
  }
  for(i in 1 : length(v3) ){
    dat[i + length(v1) + length(v2), 1] <- v3[i]
    dat[i + length(v1) + length(v2), 2] <- 'v3'
  }
  dat <- as.data.frame(dat)
  
  colnames(dat) <- c('val', 'factor')

  anova_IS <- aov(val ~ factor, data = dat)
  print(summary(anova_IS))

  anova_residuals <- anova_IS$residuals
  print(summary(anova_residuals))
}

 

⑹ 네 집단의 ANOVA

 

one_way_4_factor_anova <- function(v1, v2, v3, v4){
  dat <- matrix(0, nrow = ( length(v1) + length(v2) + length(v3) + length(v4) ), ncol = 2 )
  for(i in 1 : length(v1) ){
    dat[i, 1] <- v1[i]
    dat[i, 2] <- 'v1'
  }
  for(i in 1 : length(v2) ){
    dat[i + length(v1), 1] <- v2[i]
    dat[i + length(v1), 2] <- 'v2'
  }
  for(i in 1 : length(v3) ){
    dat[i + length(v1) + length(v2), 1] <- v3[i]
    dat[i + length(v1) + length(v2), 2] <- 'v3'
  }
  for(i in 1 : length(v4) ){
    dat[i + length(v1) + length(v2) + length(v3), 1] <- v4[i]
    dat[i + length(v1) + length(v2) + length(v3), 2] <- 'v4'
  }
  dat <- as.data.frame(dat)
  
  colnames(dat) <- c('val', 'factor')

  anova_IS <- aov(val ~ factor, data = dat)
  print(summary(anova_IS))

  anova_residuals <- anova_IS$residuals
  print(summary(anova_residuals))
}

 

 

3. 예제 [목차]

 

# http://compgenomr.github.io/book/gene-expression-analysis-using-high-throughput-sequencing-technologies.html#quantification

### step 1. 데이터 로딩 ###

library(dplyr)
# install.packages("Seurat")
library(Seurat)

# 출처 : https://github.com/compgenomr/compGenomRData
devtools::install_github("compgenomr/compGenomRData")
library(compGenomRData)
counts_file  <- system.file ("extdata/rna-seq/SRP029880.raw_counts.tsv", package = "compGenomRData")
coldata_file  <- system.file ("extdata/rna-seq/SRP029880.colData.tsv", package = "compGenomRData")
counts <- as.matrix(read.table(counts_file, header = T, sep = '\t'))


### step 2. TPM 정규화 ###
library(pheatmap)
library(stats)
library(ggplot2)

cpm <- apply(subset(counts, select = c(-width)), 2, 
             function(x) x/sum(as.numeric(x)) * 10^6)
geneLengths <- as.vector(subset(counts, select = c(width)))
rpkm <- apply(X = subset(counts, select = c(-width)),
              MARGIN = 2, 
              FUN = function(x) {
                10^9 * x / geneLengths / sum(as.numeric(x))
               })
rpk <- apply( subset(counts, select = c(-width)), 2, 
              function(x) x/(geneLengths/1000))
tpm <- apply(rpk, 2, function(x) x / sum(as.numeric(x)) * 10^6)
V <- apply(tpm, 1, var)
selectedGenes <- names(V[order(V, decreasing = T)][1:100])
colData <- read.table(coldata_file, header = T, sep = '\t', 
                      stringsAsFactors = TRUE)
M <- t(tpm[selectedGenes,])
M <- log2(M + 1)
pcaResults <- prcomp(M)
correlationMatrix <- cor(tpm)


### step 3. DESeq2 DEG 획득 ###
library(DESeq2)
library(stats)
countData <- as.matrix(subset(counts, select = c(-width)))
colData <- read.table(coldata_file, header = T, sep = '\t', 
                      stringsAsFactors = TRUE)
designFormula <- "~ group"
dds <- DESeqDataSetFromMatrix(countData = countData, 
                              colData = colData, 
                              design = as.formula(designFormula))
dds <- dds[ rowSums(DESeq2::counts(dds)) > 1, ]
dds <- DESeq(dds)
DEresults = results(dds, contrast = c("group", 'CASE', 'CTRL'))
DEresults <- DEresults[order(DEresults$pvalue),]
write.csv(DEresults, "~/Downloads/DEresults.csv")

DEresults_trim <- 0
flag = 0
for(i in 1: dim(DEresults)[1]){
  a = sum(DEresults[1:i, ]$padj < 0.05)
  b = sum(DEresults[1:i, ]$pvalue < 0.05)

  if(flag == 0 && a / b < 0.99){
    flag = 1
    DEresults_trim <- DEresults[1:i-1, ]
    write.csv(DEresults[1:i-1, ], "~/Downloads/DEresults(FDR = 1%).csv")
  }
}


### step 4. t.test DEG 획득 ###
pbmc <- CreateSeuratObject(counts = tpm, project = "pbmc3k", min.cells = 3, min.features = 200)
pbmc <- FindVariableFeatures(object = pbmc, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(x = pbmc)
markers_t.test <- FindMarkers(object = pbmc, test.use ='t', only.pos = TRUE, min.pct = 0.1, logfc.threshold = 0.25, group.by = 'orig.ident', ident.1 = 'CASE')
write.csv(markers_t.test, "~/Downloads/t.test.csv")

markers_t.test_trim <- 0
flag = 0
for(i in 1: dim(markers_t.test)[1]){
  a = sum(markers_t.test[1:i, ]$p_val_adj < 0.05)
  b = sum(markers_t.test[1:i, ]$p_val < 0.05)

  if(flag == 0 && a / b < 0.99){
    flag = 1
    markers_t.test_trim <- markers_t.test[1:i-1, ]
    write.csv(markers_t.test[1:i-1, ], "~/Downloads/t.test(FDR = 1%).csv")
  }
}


### step 5. 유전자 교집합 ###
inter <- intersect(rownames(DEresults_trim), rownames(markers_t.test_trim))
write.csv(inter, "~/Downloads/inter.csv")

 

⑴ 결과 

DESeq2 결과 7,050개의 DEG가 발견되었고, t.test 결과 9개의 DEG가 발견되었습니다. 그리고 겹치는 유전자는 9개입니다. 

⑵ 해석 

차이가 나는 이유는 기본적으로 가정이 다르기 때문입니다. DESeq2는 음이항분포(negative binomial distribution)를 가정으로 하고, t.test는 student's T test를 가정으로 하기 때문입니다. 음이항분포가 RNA-seq 데이터에 더 적합한 가정으로 알려져 있고, 표본이 크면 두 분포가 유사하지만 표본이 작으면 t.test의 부정확성이 더 두드러진다고 알려져 있습니다. 다행히도 t.test 결과 전부가 DESeq2에 발견되어 통계분석 종류에 관계 없이 robustness가 관찰된다고 할 수 있으나, t.test에서 적은 DEG가 관찰된 이유는 이런 부정확성 때문으로 여겨집니다.

 

입력: 2022.06.08 13:50