genes.results 파일이랑 isoform.results 파일 둘 다 사용할 수 있는 듯.
Gene단위 DEG 찾기.
# 파일 정보 가지는 데이터프레임 만들기.
files <- file.path(dir, paste0(samples, ".genes.results.gz"))
names(files) <- paste0("sample", 1:6)
# tximport로 불러오기.
# txIn = FALSE, txOut = FALSE 가 중요.
txi <- tximport(files, type = "rsem", txIn = FALSE, txOut = FALSE)
# row 이름이 샘플 이름, col 이름은 condition인 샘플 정보를 나타내는 데이터프레임 만들기.
condition <- factor(rep(c("A", "B"), each = 3)))
sampleTable <- data.frame(condition = condition)
rownames(sampleTable) <- colnames(txi$counts)
R
복사
DESeq2 DEG finding은 대부분 DESeqDataSet 오브젝트를 가지고 한다.
Error in DESeqDataSetFromTximport(txi, sample_table, ~condition): all(lengths > 0) is not TRUE 이런 에러가 뜨는데, txi에 estimated length가 0인 transcript가 있어서 그런 것이다.
그냥 txi$length[txi$length == 0] ← 1 로 해서 해결하라고 한다.
txi$length[txi$length == 0] <- 1
dds <- DESeqDataSetFromTximport(txi, sampleTable, ~condition)
# 보통 분석하기 전에 filtering을 거친다.
# 이건 샘플 수 등을 고려해서 자유롭게 바꿔도 될 듯.
count_threshold <- 10
keep <- rowSums(counts(dds)) >= count_threshold
dds <- dds[keep,]
# 핵심 파트.
dds <- DESeq(dds)
res <- results(dds, alpha=0.05)
summary(res) # res는 padj, pvalue, log2FoldChange 등의 column을 가진다.
# p-value histogram
hist(res$pvalue, breaks=100)
hist(res$padj, breaks=100)
R
복사
dir = '../pipelines/WTS-REAL-PROCESSING/result/03_rsem'
samples = c(
'MPN020', 'MPN212', 'MPN003', 'MPN041',
'MPN009', 'MPN122', 'MPN106',
'MPN005', 'MPN056', 'MPN208', 'MPN197', 'MPN068', 'MPN078'
)
condition = factor(c(
'LSDpos', 'LSDpos', 'LSDneg', 'LSDneg',
'LSDpos', 'LSDpos', 'LSDneg',
'LSDpos', 'LSDpos', 'LSDpos', 'LSDneg', 'LSDneg', 'LSDneg'
))
files = file.path(dir, paste0(samples, '.genes.results'))
names(files) = samples
txi = tximport(files, type='rsem', txIn=FALSE, txOut=FALSE)
# row 이름이 샘플 이름, col 이름은 condition인 샘플 정보를 나타내는 데이터프레임 만들기.
sampleTable <- data.frame(condition = condition)
rownames(sampleTable) <- colnames(txi$counts)
txi$length[txi$length == 0] = 1
dds = DESeqDataSetFromTximport(txi, sample_table, ~condition)
count_threshold = length(samples)
keep = rowSums(counts(dds)) >= count_threshold
dds = dds[keep,]
dds = DESeq(dds)
res = results(dds, alpha=0.05)
write.csv(res, '../note/0727_result/DESeq2_result.ALL.csv')
R
복사