. 1 安裝R包數據包： ALL, CLL , pasilla, airway ################################# # 1.安裝R包 ################################# source ("https:///biocLite.R" )"http://mirrors.ustc.edu.cn/bioc/" )"ALL" ,"CLL" , "pasilla" , "airway" )) #數據包 "limma" ,"DESeq2" , "clusterProfiler" )) #軟件包 "reshape2" ) #工具包 "ggplot2" ) #繪圖包 #另外可以檢測某個包是否存在，只有不存在時才會安裝 if (! require ('CLL' )){"http://mirrors.ustc.edu.cn/bioc/" )'CLL' ,ask = F , suppressUpdates = T )#意思就是設定安裝鏡像，然后隱藏安裝信息，中間過程不詢問 . 2 了解ExpressionSet對象CLL包中有data(sCLLex)，是一個表達芯片數據對象，其中包含許多信息
第一行的ExpressionSet就是表達矩陣，查看它使用exprs ，用dim 查看矩陣大小
################################# # 2.了解ExpressionSet對象 ################################# library (CLL)) #隱藏具體加載信息 #?CLL 可以看到要取出CLL數據包中的數據，需要用data函數 # CLL數據集是慢性淋巴細胞白血病(Chroniclymphocytic leukemia)，它采用了Affymetrix公司的HG_U95Av2表達譜芯片(含有12625個探針組)，共測量了22個樣品("CLL2.CEL" "CLL3.CEL" "CLL4.CEL" "CLL5.CEL" "CLL6.CEL"等)，每個樣品來自一個癌癥病人，所有病人根據健康狀態分為兩組：穩定期(Stable)組和進展期(Progressive)也稱為惡化期組 ## 2.1 取出數據 "sCLLex" )#sCLLex ## 2.2 獲得表達矩陣 #提取出來表達矩陣，賦給e ## 2.3 結構性探索 # 查看結構 # 查看前6行 #查看表達矩陣大小 #結果可以看到，包含了12625個探針，22個樣本 ## 2.4 具體內容探索【做到完完全全了解加載的對象】 #查看樣本編號 #查看所有表型變量 1 :100 ] #查看基因芯片編碼 #看看有沒有重復，去重結果還是還是12625個 #將樣本表型信息從數據框中提取出來【取出來的是表型、樣本的數據框】 2 ]) #從數據框中只要表型信息 #統計表型信息 #看看數據質量如何 0.7 )if (n.sample>40 ) par(cex = 0.5 )1.2 )"expression value" ,las=2 )使用str查看對象的結構，使用head查看對象的前6行（默認）
. 3 安裝并了解hgu95av2.db包官網： http://www./packages/release/data/annotation/html/hgu95av2.db.html
安裝 biocLite("hgu95av2.db")
這個數據庫中共有36個包，每個包都可以當成一個列表操作，可以用as.list 函數展示數據
3.1 了解hgu95av2.db 【這是一個關于hgu95av2芯片的注釋包】 # Affymetrix Human Genome U95 Set annotation data # 官方鏈接：https:///packages/release/data/annotation/manuals/hgu95av2.db/man/hgu95av2.db.pdf "hgu95av2.db" )library (hgu95av2.db)"package:hgu95av2.db" )#產生36個映射數據（探針id轉為36種主流id） #將hgu95av2包含的具體信息輸出為字符串【只用hgu95av2()可能輸出，但不易整理】 #結果可以看到，每一個子集都是有keys-value構成的，就相當于list的結構。因此對于hgu95av2的操作都要使用as.list() 1 ]) #變成列表進行查看第一組元素 #結果得到的就是 #$`1000_at` #[1] "MAPK3" #不懂hgu95av2SYMBOL是什么意思，就查找 #結果得到：hgu95av2SYMBOL is an R object that provides mappings between manufacturer identifiers and gene abbreviations【就是基因id與探針號的關系，這里顯示的就是MAPK3這個基因對應的探針id是1000_at】 3.2 探索hgu95av2CHR：探針id與染色體編號對應的關系 C = hgu95av2CHR #首先明白CHR它的數據結構，左邊一列是探針id，右邊一列是chr id #計算左側的數量，同理Rlength() #查看右側的名稱，同理Lkeys() "6" ,"8" ) #只保留chr6、chr8的數據 2 ]) #統計chr6上有668個探針，chr8上有402個探針 3.3 探索hgu95av2SYMBOL：探針id與基因名（縮寫）的關系 3.3.1 先對探針id進行簡單的探索：【輸出探針】 s = hgu95av2SYMBOL#一個函數搞定左右兩邊的各種統計【包括總數、過濾后的】 #比如：左邊探針總數是12625，能匹配上的是11430個； #右邊：基因名的總數是30071，實際上只有8585種 #【【那么，為什么基因名過濾前后差別如此大呢？】】 #[下面]目的：對參考注釋進行過濾，去掉沒有map的探針 #手動過濾，然后進行統計 #統計：與summary結果相符，有11430個探針有對應基因名稱的【剩余1165沒有名稱】 3.3.2 再對基因名進行探索【匹配、列表、轉數據框、查特定基因、基因數】 ss = as.list(s[mapped_probes]) #將匹配上的基因輸出為列表 #toTable把list轉換成數據框 【后面會經常用到sst】 #看一下sst的列名【相當于summary(s)中的Lkeyname=probe_id、Rkeyname=symbol】 "^TP53$" ,sst$symbol),] #查找TP53基因對應的探針id #查看總共基因數（注意unique的使用，為了避免重復基因出現） #table將基因名字符量化為數字，sort從小到大排序，tail取最大6個 #table函數統計不同探針數量對應的基因數 #【【為什么有5個基因會分別有8個探針，而大部分6555個基因只對應一個探針？】】 為什么有5個基因會分別有8個探針，而大部分6555個基因只對應一個探針？
A：不管是Agilent芯片，還是Affymetrix芯片，上面設計的探針都非常短。最長的如Agilent芯片上的探針，往往都是60bp，但是往往一個基因的長度都好幾Kb。因此一般多個探針對應一個基因，取最大表達值探針來作為基因的表達量
. 4 過濾、整合表達矩陣4.1 過濾 #[下面]目的：對表達矩陣進行過濾，去掉沒有map的探針 # table(rownames(e)%in%sst$probe_id) #找到sCLLex表達矩陣(e)在hgu95av2.db包中沒有交叉的探針 # %>%是管道符號，相當于linux的“|”; %in%表示兩者求交集 in %mapped_probes,] #對原始表達矩陣過濾 # e2 = e[match(rownames(e), mapped_probes, nomatch = 0),] #使用match過濾 4.2 整合 #[下面]目的：一個基因對應一個探針 #現狀分析：多個探針對應一個基因的情況存在 #解決途徑：只保留在所有樣本里面平均表達量最大的那個探針【一般采用均值即可，當然也可用最大值】 #（之前得到的sst矩陣來的時候就是根據mapped_probes得到的，也就是過濾好的，因此可以直接用作index） function (x) rownames(x)[which.max(rowMeans(x))]) #矩陣用by，向量用tapply #獲得每個基因獨特的探針 in %uniprobes,] #完成表達矩陣過濾 2 ]4.3 重塑【reshape2:矩陣=》數據框】 head(efilt) #可以看到現在還是一個矩陣的樣子，行為基因名，列為樣本名 #我們想把這個矩陣變成標準的tidy data，三列：第一列基因名，第二列樣本名，第三列表達量 library (reshape2)#先將原來矩陣“融化 'symbol' ,'sample' ,'value' ) #重新命名三列 #最后再加一列表型信息（就是剛得到表達矩陣時提取的group_list）,每一個樣本有8585行【nrow(efilt)統計得到】，我們這里重塑的數據框是把所有樣本按次序堆疊起來，因此每個樣本的表型應該將group_list中對應的表型重復nrow(efilt)這些次 #例如，第一個樣本是CLL11.CEL，它的表型是progres。現在我們把CLL11.CEL放到重塑的sample列，應該就有8585行都是CLL11.CEL，然后再有8585個樣本CLL12.CEL。因此CLL11.CEL對應的表型信息也就應該是8585個，也就是融化前矩陣的行數 #融化前后做個對比就看出來了 . 5 畫圖探索作出樣本和基因表達量之間的關系圖，主要基于ggplot
5.1 boxplot 5.2 violin plot 5.3 Histogram 5.4 Density . 6 做一些統計mean,median,max,min,sd,var,mad，t檢驗
6.1 利用apply函數進行統計 他需要矩陣，也就是之前得到的efilt，按行進行統計 即可
e_mean = tail(sort(apply(efilt,1 ,mean)),30 )1 ,median)), 30 )1 ,max)),30 )1 ,min)),30 )1 ,sd)),30 )1 ,var)),30 )1 ,mad)),30 ) #絕對中位差來估計方差,先計算出數據與它們的中位數之間的偏差，然后這些偏差的絕對值的中位數就是mad 6.2 看表達量top30基因之間的重合部分 用UpSetR包結合之前做的top30基因各種統計，適用于樣本數量大于5的情況
其實我們平常做的韋恩圖也是這個意思，找交集，但是韋恩圖樣本數量一般都會控制在
install.packages("UpSetR" )library ("UpSetR" )#先找出7個統計量的共同基因名 #分別將7個統計量關于共同基因名的統計值 in % names(e_mean) ,1 ,0 ),in % names(e_median) ,1 ,0 ),in % names(e_max) ,1 ,0 ),in % names(e_min) ,1 ,0 ),in % names(e_sd) ,1 ,0 ),in % names(e_var) ,1 ,0 ),in % names(e_mad) ,1 ,0 )7 ,sets.bar.color = "#56B4E9" )6.3 批量T檢驗——為了得到pvalue【后續分析重點】 有了pvale就能有padj值；另外還需要對照、處理兩組均值，這樣就能有log2FC
gl=as.factor(group_list) #將最初得到的的表型數據因子化 1 ]) #levels(group_list)[1]返回第一個因子progres,從group_list中選出progres的元素，用which來獲取他們所在的位置【目的是為下面分別得到兩種表型的樣本作準備】 2 ]) #返回第二個因子stable #將表型為progres的樣本選出來，因為這里是要求t值，可以命名為e矩陣的t值，即et #將表型為stable的樣本選出來 #按列合并 1 , function (x){# 多組樣本的t檢驗 "BH" ) #多重比較時校正p值 #progres是對照組 #stable是使用藥物處理后的——處理組 4 ]),] #從小到大排序 一般來講，下游分析使用的差異表達矩陣應該是log2后的結果 ，它的計算公式是log2FC = log2 (mean(處理組/對照組))
這里為什么可以之間相減？ 芯片數據的特點 ：小樣本和大變量，因此數據分布呈偏態、標準差大。而對數轉換能使上調、下調的基因連續分布在0的周圍，更加符合正態分布，同時對數轉換可以使熒光信號強度的標準差減少，方便下游分析
. 7 做一些分析 表達量、聚類、PCA、火山圖7.1 按mad指標選表達量前30（top30）的基因,做熱圖可視化 做熱圖前需要將矩陣數據中心化+標準化【目的為了向數據中心靠攏，減小數據之間的差別】
中心化：數據減去均值后得到的; 標準化則是在中心化后的數據基礎上再除以數據的標準差
top30_gene=names(e_mad)#得到top30的表達量矩陣 7.2 聚類分析圖 過濾后的樣本聚類 過濾后的基因聚類 使用factoextra包 如果看到聚類分析的結果枝干太長，那么就要換種聚類方法
7.3 PCA分析 【目的：簡化變量的個數】本質是降維，本來應該有22個變量，現在我們變成了22個主成分,一般前面的幾個主成分就能解釋所有的數據。解釋：https://wenku.baidu.com/view/c22d1539a31614791711cc7931b765ce05087a6f.html
7.3.1 使用一鍵式ggfortify + prcomp 關于ggfortify的使用：https://wenku.baidu.com/view/e5dc63fb763231126fdb1100.html
結果中：不同顏色代筆不同分組;
7.3.2 使用fviz_pca_ind包 7.4 差異分析火山圖【也就是logFC加上-log10(Pvals)的散點圖】 首先構建差異分析矩陣 suppressMessages(library (limma)) #加載包不顯示啰嗦的信息 關于limma：
limma是基于R和Bioconductor平臺的分析芯片數據的綜合軟件包，該包功能齊全、教程完善、使用率極高，幾乎成為了芯片數據處理流程的代名詞；
關于比較矩陣: 參考https://github.com/bioconductor-china/basic/blob/master/makeContrasts.md
7.4.1 構建一個非差異比較矩陣【只有一組處理和對照，所以可以不用比較矩陣】 design1=model.matrix(~factor(group_list)) 3 ) #設置結果的小數位數為3 2 ,adjust='BH' ,n=Inf ) #coef要么必須等于2， 要么是個字符串；關于adjust的設置，說明書中13.1章有描述，BH是最流行的設置 # topTable 默認顯示前10個基因的統計數據；使用選項n可以設置,n=Inf就是不設上限，全部輸出 #去除缺失值 7.4.2 構建差異比較矩陣 #【當有多組處理、對照時就發揮作用了，例如兩個細胞系A、B，每個細胞系都有處理、對照組】 ##先分組## 0 +factor(group_list))##后比較## "-" ),levels = design2) #注意寫法：處理和對照之間用-隔開 #假如是A、B兩組處理、對照的情況，makeContrasts(A_trt-A_con,B_trt-B_con),levels = design2) #另外，paste0相當于paste(..., sep = "", collapse = ) 將...對象拼起來，并且sep指定對象們無縫拼接，得到的結果用collapse分隔 # 如果paste(unique(group_list),sep = "-") =》"progres." "stable" # 如果paste0(unique(group_list),collapse = "-") =》"progres.-stable" 【一個整體，只是中間有分隔】 ##得矩陣## 1 , n=Inf ) #如果是多組，還是上面的例子：如果需要第一組A的差異基因，就用coef = 1；如果需要第二組B的差異基因，就用coef = 2；但是如果是多組比較，結果容易什么差異基因都沒有，要注意批次效應！使用h-cluster來辨識批次效應，用combat來校正 #write.csv(nrDEG2,"limma_notrend.results.csv",quote = F) #就結果來看，DEG_mtx1和DEG_mtx2是一樣的 ######################################################### 7.4.3 準備火山圖 畫火山圖第一步，設定閾值，選出UP、DOWN、NOT表達基因
DEG=DEG_mtx22 *sd(abs(logFC)) ) # with()函數適用于當同名變量出現多次，避免程序定位錯誤的情況，with()的括號內外，信息是完全隔絕的。因為需要使用DEG中的logFC，但logFC只在DEG中出現，當然可以單獨用DEG$logFC選出來設為全局變量，但是感覺有點“興師動眾”，此時用with就是最好的情況 關于閾值的設置：一般情況都是設為2，但具體背景不同還是應該設置不同的閾值。根據網站https://www./about-bmj/resources-readers/publications/statistics-square-one/2-mean-and-standard-deviation 的介紹。mean+2SD可以反映95%以上的觀測值，這是比較可信的，如果再嚴格一點，設為mean+3SD，就可以反映97%以上的觀測
DEG$result = as.factor(ifelse(DEG$P.Value < 0.05 & abs(DEG$logFC) > logFC_cutoff,'UP' ,'DOWN' ),'NOT' )#這是兩個ifelse判斷嵌套。先了解ifelse的結構，ifelse（條件，yes，no），如果滿足條件，那么返回yes/或者執行yes所處的下一個命令；反之返回no #這里外層的ifelse中DEG$P.Value < 0.05 & abs(DEG$logFC) > logFC_cutoff是判斷條件，這個就是看p值和logFC是不是達到了他們設定的閾值【p是0.05，logFC是logFC_cutoff】，如果達到了就進行下一個ifelse，達不到就返回NOT； #第二層ifelse也是上來一個條件：DEG$logFC > logFC_cutoff，如果達到了，就返回UP即上調基因，達不到就是下調DOWN # 最后將判斷結果轉位因子型，得到DOWN、UP、NOT的三種因子 畫火山圖第二步，設定圖形的標題
this_tile <- paste0('Cutoff for logFC is ' ,round(logFC_cutoff,3 ), #round保留小數位數 '\nThe number of up gene is ' ,nrow(DEG[DEG$result =='UP' ,]) ,'\nThe number of down gene is ' ,nrow(DEG[DEG$result =='DOWN' ,])最后才開始畫火山圖
library (ggplot2)#注意DEG的列名不同的分析軟件可能命名不同，比如p值，有的是P.value，有的是P_value #選出p值并且進行對數轉換 #其實火山圖就是不同顏色的散點圖，只不過低于閾值的居多，并且設為深顏色，所以看著像火山噴發 0.4 , size=1.75 ) +20 )))+"log2 fold change" ) + ylab("-log10 p-value" ) +15 ,hjust = 0.5 ))+'blue' ,'black' ,'red' )) ## 這里要注意和之前設置的result三個因子相對應，DOWN就設為blue，NOT就設為black
