>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
#数据标准化处理
rld <- rlogTransformation(dds) #DEseq2自己的方法标准化数据`
exprSet_new=assay(rld) #提取DEseq2标准化后的数据
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自带函数
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
