對(duì)day49，day50兩次學(xué)習(xí)的總結(jié)赘理。
拿到reads表達(dá)矩陣后宦言，在自己電腦上用DESeq2這個(gè)R包進(jìn)行數(shù)據(jù)分析。
一商模、先要準(zhǔn)備表達(dá)矩陣和分組表格
1奠旺，表達(dá)矩陣
count之后的txt文件進(jìn)行整理，只留下geneid和四個(gè)樣本列施流。
其中有一些行不完整响疚，有NA的數(shù)據(jù)，在后面讀入的時(shí)候一定要用fill=TURE參數(shù)瞪醋。
2忿晕，分組表格
可以用txt設(shè)置分組表格，但是可能是用WPS建立表格趟章，轉(zhuǎn)成txt有格式問(wèn)題杏糙，一直報(bào)錯(cuò)：“Error in make.names(col.names, unique = TRUE) :
invalid multibyte string at '<ff><fe>s'
In addition: Warning messages:
1: In read.table(file = file, header = header, sep = sep, quote = quote, :
line 1 appears to contain embedded nulls”
最后解決辦法是，轉(zhuǎn)成scv格式蚓土，用read.csv 命令讀入宏侍。group_list一列的列名必須是“group_list”。
二蜀漆、在R中編輯谅河，依此運(yùn)行下面的腳本：
#step1: read counts
rm(list = ls()) #首先清空R中的變量
setwd("/Users/meraner/Documents/bioinfo")
a=read.table('/Users/meraner/Documents/bioinfo/counts2.txt', fill=T, header=T, row.names=1, stringsAsFactors = F)
a = na.omit(a) #含有NA值的行刪除
index1=c(rowSums(a)>0)
exprSet=a[index1,]
write.csv(exprSet,'filter_reads_matrix.csv' ) #過(guò)濾后的數(shù)據(jù)保存
說(shuō)明
由于count文件中可能會(huì)有很一些行里面有缺省NA，就是沒(méi)比對(duì)成功的。所以用fill=TURE绷耍，即可讀入吐限，這些數(shù)據(jù)顯示為NA。在Mac電腦的R里面需要用fill=T肮邮肌诸典！不然會(huì)報(bào)錯(cuò)： object 'TURE' not 。
接著可以刪除有NA的行崎苗。即用no.omit命令把NA的行去除狐粱，避免后面分析中的干擾。這些NA數(shù)據(jù)的基因名通常都很奇怪胆数，很長(zhǎng)一串后續(xù)分析有點(diǎn)兒煩肌蜻，去除更好。
index1為邏輯型必尼。就是計(jì)算每一行的和蒋搜，得到的數(shù)大于零就是ture，等于零就是false判莉，再把ture的那些行單獨(dú)拎出來(lái)建立一個(gè)expreSet表達(dá)框豆挽，就是過(guò)濾過(guò)后的數(shù)據(jù)，用來(lái)進(jìn)行后面的分析骂租。
#step2: set group & build dds
library(DESeq2)
coldata <- read.csv ('/Users/meraner/Documents/bioinfo/group.csv', row.names=1, stringsAsFactors =T)
group_list=coldata[,1] #取這一列的數(shù)據(jù)為變量group_list
colnames(exprSet)<-rownames(coldata) #把表達(dá)矩陣的列名重設(shè)為簡(jiǎn)單格式且和coldata里的樣本名保持一致
dds <- DESeqDataSetFromMatrix(countData = exprSet, colData = coldata, design = ~group_list) #這個(gè)命令創(chuàng)建DEseqDataSet（dds）祷杈。
說(shuō)明
創(chuàng)建一個(gè)數(shù)據(jù)框coldata，只有一列內(nèi)容渗饮，即分組信息但汞。行名為樣本名，順序必須和前面表達(dá)矩陣從左到右的順序一致互站。
建立一個(gè)變量group_list（名字必須是這個(gè)）私蕾，里面內(nèi)容就是每個(gè)樣本的分組信息，順序必須和前面表達(dá)矩陣從左到右的順序一致胡桃。這個(gè)變量還必須是一個(gè)因子型的踩叭。所以用了stringsAsFactors =T。
countData用于說(shuō)明數(shù)據(jù)來(lái)源翠胰，colData用于說(shuō)明不同組數(shù)據(jù)的實(shí)驗(yàn)操作類(lèi)型容贝，design用于聲明自變量，即誰(shuí)和誰(shuí)進(jìn)行對(duì)比
#step3: DEG
dds2 <- DESeq(dds)
res <- results(dds2,contrast = c("group_list","siSUZ12","con"))
resOrdered<- res[order(res$padj),]
DEG=as.data.frame(resOrdered)
write.csv(DEG,file = "treat_vs_con.csv")
DEGsig <-subset(res, padj < 0.001 & abs(log2FoldChange) > 1)
write.csv(DEGsig,file= "treat_vs_con_sig.csv")
說(shuō)明
用DESeq命令對(duì)dds數(shù)據(jù)進(jìn)行運(yùn)算及分析,共三步：文庫(kù)大小估計(jì)之景；離散程度估計(jì)斤富；統(tǒng)計(jì)檢驗(yàn)。生成結(jié)果文件是按照siSUZ12組vs. con組進(jìn)行比較的锻狗。#order函數(shù)是給出從小到大排序后的位置(默認(rèn)升序), 將結(jié)果文件按照res文件中的padj這一列進(jìn)行降序排列满力，其中$符號(hào)表示res中的padj這一列焕参。再按照自定義閾值提取差異基因。
#step4：normalization matrix
rld <- rlogTransformation(dds2) ## 得到經(jīng)過(guò)normlization的表達(dá)矩陣油额。
exprSet_new=assay(rld) #生成歸一化的矩陣文件
write.csv(exprSet_new,'Normalized_matrix.csv' ) #導(dǎo)出數(shù)據(jù)叠纷，保存為csv格式。
說(shuō)明
rlogTransformation是DESeq2包中的一個(gè)命令（函數(shù)）潦嘶，在help文檔中查看涩嚣，它能把原始含有reads的表達(dá)數(shù)值，根據(jù)整個(gè)樣本和每個(gè)基因的長(zhǎng)度等衬以，計(jì)算出相對(duì)表達(dá)缓艳。并表示為log2格式，因此最后導(dǎo)出的數(shù)據(jù)里會(huì)有負(fù)數(shù)看峻。比如-1.7，就是2^（-1.7）=0.3衙吩。
assay是SummarizedExperiment包中的函數(shù)互妓，可以把陣提取出來(lái)。
#step5: visulization
#fig1
hist(exprSet_new)
#fig2
tiff(filename = "exprBOX.tiff",width = 1500, height = 1500, units = "px", res = 150)
par(cex = 1.2) #指定符號(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_new, col = cols,main="expression value",las=2)
dev.off()
#fig3
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()
}
#fig4
tiff(filename = "PCA.tiff", width = 1500, height = 1500, units = "px", res = 150)
plotPCA(rld, intgroup="group_list")
dev.off()
說(shuō)明
hist() 用來(lái)做直方圖（柱形圖）坤塞，應(yīng)該用歸一化的數(shù)據(jù)來(lái)看冯勉。
boxplot是用來(lái)看每個(gè)樣本的整體表達(dá)水平。
plotCounts是用來(lái)看單個(gè)感興趣基因的表達(dá)差異
plotPCA用來(lái)做PCA聚類(lèi)分析
生成的圖片都保存在工作目錄中摹芙，就是前面setwd的那個(gè)路徑灼狰。
