接著day49
一改含、看看整體表達(dá)情況
原始矩陣為變量exprSet蝶涩,是從txt文件讀入得到的數(shù)據(jù)框
歸一化后矩陣為變量exprSet_new骤铃,用DESeq2進(jìn)行歸一化后計(jì)算而來沛善,是數(shù)值型矩陣蹲坷。從下面class()命令可以看出兩個(gè)變量屬性不同驶乾。
class(exprSet)
[1] "data.frame"
class(exprSet_new)
[1] "matrix" "array"
exprSet<-as.matrix(exprSet) #把data.frame格式轉(zhuǎn)唯matix格式
hist(exprSet)
hist(exprSet_new)
說明
hist() 用來做直方圖(柱形圖)邑飒,如果是一列數(shù)據(jù),就會(huì)畫成從小到大那樣的柱子级乐。是一個(gè)矩陣幸乒,就把全部數(shù)據(jù)按照大小和頻率列成柱形圖〈侥粒可以看出歸一化后的數(shù)據(jù)看著更均勻。
看一下這兩個(gè)矩陣的樣子:(這兩個(gè)矩陣都是把變量用write.csv給導(dǎo)出的聚唐。)
二丐重、看看不同樣本之間表達(dá)情況:
par(cex = 0.7) #指定符號(hào)的大小
par(mar=c(10, 5, 5, 5))
n.sample=ncol(exprSet) #樣本數(shù)
if(n.sample>40) par(cex = 0.5)
n.sample=ncol(exprSet)#樣本數(shù)
cols <- rainbow(n.sample*1.2)
boxplot(exprSet, col = cols,main="expression value",las=2)
boxplot(exprSet_new, col = cols,main="expression value",las=2)
說明
par函數(shù)cex參數(shù)是用來控制文字和點(diǎn)的大小。cex參數(shù)大了表明散點(diǎn)圖里的點(diǎn)和各種文字部分的字體都大杆查,比如上面命令中設(shè)定為0.7扮惦,在樣本很多的時(shí)候要縮小些字體,所以有了>40亲桦,改為0.5這個(gè)命令行崖蜜。
cex對(duì)所有的文字進(jìn)行統(tǒng)一設(shè)置外,針對(duì)不同的標(biāo)題客峭,還有對(duì)應(yīng)的cex系列參數(shù)的用法很多豫领,比如調(diào)節(jié)主標(biāo)題,副標(biāo)題字體大小舔琅,調(diào)節(jié)坐標(biāo)軸的字體大小等等恐,可以看下面鏈接。
參考教程:https://www.bbsmax.com/A/A7zg1Erl54/
而par函數(shù)里的mar參數(shù)則是設(shè)定圖形四周的留空大小备蚓。
參考教程:https://www.yht7.com/news/136868
兩個(gè)圖都是箱型圖课蔬,可歸一化以后的更好看些,在有些測(cè)序結(jié)果里看到過郊尝,主要是看各個(gè)樣本表達(dá)水平有沒有大的不同二跋。如果某個(gè)樣本比別人差太多,那就是有問題流昏,還有可能是這個(gè)細(xì)胞有什么特殊處理扎即。
原始數(shù)據(jù)那個(gè)圖太丑了。横缔。
三铺遂、查看單一基因在組間差異
plotCounts(dds2, "CD38", intgroup = "group_list") #把前面dds2中的nomalization之后的某個(gè)特定基因表達(dá)水平展示出來。
說明
plotCounts也是DESeq2包中的一個(gè)命令茎刚,dds是DESeqDataSet.襟锐, gene是一個(gè)特殊基因名稱,intgroup:在colData(x)中膛锭,用于進(jìn)行分組的名稱粮坞。
如果有很多想看的基因蚊荣,可以批量執(zhí)行。
sigGene= read.csv('treat_vs_con_sig.csv',row.names=1) #讀取文件
for(i in rownames(sigGene)){
genename=paste(i,sep="") #基因名
pdf(file=paste(i,'_counts.pdf',sep="")) #以pdf格式輸出
plotCounts(dds2, genename, intgroup = "group_list")
dev.off() #不在Rstudio里輸出莫杈,直接保存在當(dāng)前工作路徑下
}
很快就在工作目錄下生成了這么多文件互例。
挑一個(gè)打開看看。這個(gè)數(shù)值是歸一化后的counts值筝闹,找來表格里的數(shù)對(duì)照一下媳叨,還真沒錯(cuò)。如果樣本多一些关顷,畫的散點(diǎn)圖應(yīng)該挺好看的糊秆。
上面這個(gè)表里的數(shù)是log2之后的,所以原來的數(shù)據(jù)要用2(x)來計(jì)算一下议双。比如那個(gè)9.479736的痘番,轉(zhuǎn)成2(9.479736)=713.9781。
三平痰、PCA聚類
tiff(filename = "PCA.tiff", width = 1500, height = 1500, units = "px", res = 150)
plotPCA(rld, intgroup="group_list")
dev.off() #不在Rstudio里輸出汞舱,直接保存在當(dāng)前工作路徑下這個(gè)PCA.tiff文件了