學(xué)習(xí)資料：B站視頻：生物技能樹 第八季轉(zhuǎn)錄組測(cè)序數(shù)據(jù)分析
得到表達(dá)矩陣之后就要考慮如何知道到差異表達(dá)基因（DEG）等脂，這個(gè)要在R中操作栽烂，具體是需要airway、DESeq2昙衅、edgeR這三個(gè)包度帮。
airway這個(gè)包是數(shù)據(jù)包歼捏，但我覺得沒啥用，它那些個(gè)命令只能用于它自己的數(shù)據(jù)笨篷，用自己從GEO上搞來的數(shù)據(jù)不香么瞳秽？雞肋。
DESeq2和EdgeR都是做基因差異表達(dá)分析的率翅，主要用于RNA-Seq數(shù)據(jù)寂诱，有些研究者用來處理類似的ChIP-Seq, shRNA、質(zhì)譜數(shù)據(jù)等安聘，原理都是在不同處理因素的組別之間找差異。
一瓢棒、安裝R包
1浴韭，執(zhí)行下面命令：
install.packages("BiocManager") #先安裝biocmanager
library(BiocManager)
BiocManager::install ('airway')#這個(gè)包安裝很費(fèi)勁，報(bào)錯(cuò)脯宿，但實(shí)際上也沒用念颈。后悔安裝之。
BiocManager::install ('DESeq2')
BiocManager::install ('edgeR')
Update all/some/none? [a/s/n]: 回答a连霉。
2榴芳，安裝airwa的時(shí)候出現(xiàn)報(bào)錯(cuò)：
解決辦法一——失敗：
remove.packages("cluster", lib="C:/Program Files/R/R-4.1.3/library")
remove.packages("MASS", lib="C:/Program Files/R/R-4.1.3/library")
remove.packages("Matrix", lib="C:/Program Files/R/R-4.1.3/library")
remove.packages("mgcv", lib="C:/Program Files/R/R-4.1.3/library")
remove.packages("nlme",lib="C:/Program Files/R/R-4.1.3/library")
remove.packages("survival", lib="C:/Program Files/R/R-4.1.3/library")
重新BiocManager::install ('airway')跺撼，依然報(bào)錯(cuò)窟感。
解決辦法二——成功：
下載該軟件到本地，再從右下角的install按鈕安裝歉井，如圖：
3柿祈，卸載R包
由于后來發(fā)現(xiàn)airway沒用，就練習(xí)了一下卸載。
remove.packages('airway')
二躏嚎、文件準(zhǔn)備
1蜜自，表達(dá)矩陣
通過前面計(jì)數(shù)后得到的txt文件，在excel里修整成為如下表卢佣，保留列名為樣本名重荠，行名為geneid。
2虚茶，分組準(zhǔn)備
建立txt文件
三戈鲁、讀入并對(duì)表達(dá)矩陣進(jìn)行過濾
rm(list = ls()) #首先清空R中的變量
setwd("D:/work/NGSanalysis/rnaseq")
a=read.table('D:/work/NGSanalysis/rnaseq/counts2.txt', header=T, row.names=1, stringsAsFactors = F) #行名，列名都設(shè)定了
index1=c(rowSums(a)>0) #刪除都是零的基因刪除
exprSet=a[index1,] #建立一個(gè)表達(dá)矩陣
write.csv(exprSet,'filter_reads_matrix.csv' )#過濾后的數(shù)據(jù)
四媳危、差異分析：接上面步驟
library(DESeq2)
group=read.table('D:/work/NGSanalysis/rnaseq/group.txt', row.names=1,stringsAsFactors =T )#讀入分組文件荞彼，為因子，有行名（樣本名）
group_list=group[,1] #取這一列的數(shù)據(jù)為變量group_list
coldata <- data.frame(row.names =colnames(exprSet), group_list) #創(chuàng)建一個(gè)數(shù)據(jù)框coldata待笑，只有一列內(nèi)容是分組信息鸣皂，行名為樣本名，
dds <- DESeqDataSetFromMatrix(countData = exprSet, colData = coldata, design = ~group_list) #這個(gè)命令創(chuàng)建DEseqDataSet（dds）暮蹂。countData用于說明數(shù)據(jù)來源寞缝，colData用于說明不同組數(shù)據(jù)的實(shí)驗(yàn)操作類型，design用于聲明自變量仰泻，即誰和誰進(jìn)行對(duì)比
dds2 <- DESeq(dds) #用于對(duì)dds數(shù)據(jù)進(jìn)行運(yùn)算及分析,共三步：文庫大小估計(jì)荆陆；離散程度估計(jì)；統(tǒng)計(jì)檢驗(yàn)
res <- results(dds2,contrast = c("group_list","siSUZ12","con")) #生成結(jié)果文件,siSUZ12組vs. con組集侯。
resOrdered<- res[order(res$padj),] #order函數(shù)是給出從小到大排序后的位置(默認(rèn)升序), 將結(jié)果文件按照res文件中的padj這一列進(jìn)行降序排列被啼，其中$符號(hào)表示res中的padj這一列
DEG=as.data.frame(resOrdered)
write.csv(DEG,file = "treat_vs_con.csv") #在工作目錄中生成結(jié)果文件
#按照自定義閾值提取差異基因
DEGsig <-subset(res, padj < 0.001 & abs(log2FoldChange) > 1)
write.csv(DEGsig,file= "treat_vs_con_sig.csv")
五、用DESeq2歸一化的表達(dá)矩陣
rld <- rlogTransformation(dds2) ## 得到經(jīng)過normlization的表達(dá)矩陣棠枉。
exprSet_new=assay(rld) #生成矩陣文件
write.csv(exprSet_new,'Normalized_matrix.csv' ) #導(dǎo)出數(shù)據(jù)浓体，保存為csv格式。
說明
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ù)扯躺，可以把陣提取出來冲茸。
