>library(DESeq2)
>raw_count_filt<-read.table("HTseq.QC.sort.by.n.count",header = T,row.names = 1)
> head(raw_count_filt)
KKO.mLtb1 KO.muta KO.mutc KO.WTa KO.WTb KO.WTc
0610006L08Rik 1 0 0 0 0 0
0610007P14Rik 999 1234 1293 1234 1663 1270
0610009B22Rik 359 393 274 381 385 288
0610009E02Rik 4 6 1 6 3 4
0610009L18Rik 8 3 7 11 6 10
0610009O20Rik 635 550 561 564 692 493
#構建dds對象
condition <- factor(c('treated','treated','treated','untreated','untreated','untreated'))
dds <- DESeqDataSetFromMatrix(countData, DataFrame(condition), design= ~ condition )
> dds
class: DESeqDataSet
dim: 48587 6
metadata(1): version
assays(3): counts mu cooks
rownames(48587): 0610006L08Rik
0610007P14Rik ... n-TSaga9 n-TStga1
rowData names(27): baseMean baseVar ...
deviance maxCooks
colnames(6): KKO.mLtb1 KO.muta ... KO.WTb
KO.WTc
colData names(2): condition sizeFactor
#數(shù)據(jù)標準化處理
rld <- rlogTransformation(dds) #DEseq2自己的方法標準化數(shù)據(jù)`
exprSet_new=assay(rld) #提取DEseq2標準化后的數(shù)據(jù)
write.table(exprSet_new, file="FZH.DESeq2.normalization.txt", sep="\t",quote=F)
> head(exprSet_new)
KKO.mLtb1 KO.muta KO.mutc KO.WTa KO.WTb KO.WTc
0610006L08Rik -1.967768 -1.970802 -1.970600 -1.970944 -1.970991 -1.970597
0610007P14Rik 10.076339 10.268512 10.434751 10.183503 10.430670 10.419946
0610009B22Rik 8.440090 8.509859 8.344399 8.417481 8.402646 8.380791
0610009E02Rik 1.971115 1.985315 1.951507 1.981837 1.961226 1.975130
0610009L18Rik 2.902465 2.867296 2.902476 2.917000 2.883140 2.925268
0610009O20Rik 9.237139 9.120903 9.246608 9.065483 9.208009 9.142025
plotPCA(rld, intgroup=c('condition')) #DEseq2自帶函數(shù)
dev.copy(png,'Deseq2_pca.png')
dev.off()
Deseq2_pca.png
通過上圖驴一，可以看見組層次樣本相似信息篙程，無法精確到具體是哪一個樣本浅役？如果想要精確到具體是哪一個樣本（圖中點各自代表哪一個樣本）绘梦，該怎么做呢嗤军？
>pcaData <- plotPCA(rld, intgroup=c("condition"), returnData=TRUE)
> pcaData
PC1 PC2 group condition name
KKO.mLtb1 -9.3479355 -0.4078916 treated treated KKO.mLtb1
KO.muta 2.3993577 4.0581652 treated treated KO.muta
KO.mutc 0.8091405 2.3055847 treated treated KO.mutc
KO.WTa 1.0457699 0.2362714 untreated untreated KO.WTa
KO.WTb 2.1439032 -3.2297884 untreated untreated KO.WTb
KO.WTc 2.9497641 -2.9623413 untreated untreated KO.WTc
plot(pcaData[,1:2],pch=19,col=c("red","red","red","blue","blue","blue"))
text(pcaData[,1],pcaData[,2]+0.2,row.names(pcaData),cex=0.5)
FZH.deseq2_pca.name.png
