WGCNA分析，簡(jiǎn)單全面的最新教程詳細(xì)介紹了WGCNA構(gòu)建共表達(dá)網(wǎng)絡(luò)的原理和過(guò)程，但有時(shí)由于樣本太少，效果不理想，找不到合適的Power值。實(shí)驗(yàn)室?guī)煹軒熋媒?jīng)常在討論樣本較少的時(shí)候如何篩選共表達(dá)基因（轉(zhuǎn)錄因子的突變體/轉(zhuǎn)基因材料和野生型材料在2-3個(gè)時(shí)間點(diǎn)進(jìn)行轉(zhuǎn)錄組測(cè)序），加上3個(gè)生物學(xué)重復(fù)勉強(qiáng)構(gòu)成18個(gè)樣本。這樣的實(shí)驗(yàn)設(shè)計(jì)，已經(jīng)足夠算出各個(gè)時(shí)間點(diǎn)或者處理間的差異基因和任意基因之間的相關(guān)系數(shù)。本文介紹少樣本如何構(gòu)建共表達(dá)無(wú)向網(wǎng)絡(luò)，篩選共表達(dá)基因的流程。

授權(quán)轉(zhuǎn)載自SimonCat，有增刪。

0 實(shí)驗(yàn)設(shè)計(jì)

思路：（1）得到脅迫處理誘導(dǎo)的基因和突變體野生型間的差異基因作為候選基因集（2）兩兩計(jì)算候選基因的皮爾遜相關(guān)系數(shù)（3）使用igrah計(jì)算各個(gè)基因的中心度。

從RNA-seq開始構(gòu)建共表達(dá)網(wǎng)絡(luò)。

1 使用軟件FastQC進(jìn)行測(cè)序質(zhì)量評(píng)估 (NGS基礎(chǔ) - FASTQ格式解釋和質(zhì)量評(píng)估)

$ ls *.fastq.gz | xargs -n1 -P2 fastqc $1

備注： -n1 表示一個(gè)參數(shù)傳遞給fastqc，-P表示線程數(shù)

2 使用trimmomatic質(zhì)量清理

script.trimmo.sh，對(duì)測(cè)序數(shù)據(jù)進(jìn)行質(zhì)量控制

$ echo 'nohup java -jar trimmomatic-0.36.jar SE $1 $1.trim.fil.gz LEADING:20 TRAILING:20 AVGQUAL:25 SLIDINGWINDOW:10:30 MINLEN:36 > $1.trim.nohup' > script.trimmo.sh

3 使用hisat2比對(duì) (轉(zhuǎn)錄組分析工具哪家強(qiáng)？)

使用比對(duì)軟件hisat2建立索引和比對(duì)

$ hisat2-build TAIR10_Chr.all.fasta TAIR10_Chr.all
$ echo ’nohup $WD/tools/hisat2-2.0.5/hisat2 -x TAIR10_Chr.all -U $1 -S $1.sam > $1.align.stat.txt’ > script.align.sh
$ ls *.trim.fil.gz | xargs -n1 -P2 sh script.align.sh

4使用htseq-count，獲得每個(gè)基因的counts數(shù)目

$ echo 'htseq-count -s no $1 TAIR10.gtf > $1.count ' > script_count.sh
$ ls *sam|xargs -n1 sh script_count.sh
# 注意用grep -v 去除掉htseq-count生成文件中的無(wú)用信息
$ csvtk join -t treat_1h_rep1.count treat_1h_rep2.count treat_1h_rep3.count control_1h_rep1.count  control_1h_rep2.count  control_1h_rep3.count > 0h_count.txt
$ csvtk join -t treat_6h_rep1.count treat_6h_rep2.count treat_6h_rep3.count control_6h_rep1.count  control_6h_rep2.count  control_6h_rep3.count > 6h_count.txt
$ csvtk join -t treat_12h_rep1.count treat_12h_rep2.count treat_12h_rep3.count control_12h_rep1.count  control_12h_rep2.count  control_12h_rep3.count > 12h_count.txt
$ csvtk join -t 0h_count.txt 6h_count.txt 12h_count.txt > total_count.txt

備注：PE數(shù)據(jù)必須先根據(jù)序列的名字或者位置進(jìn)行排序，本次項(xiàng)目是單端，所以隨意。

5 使用DESeq2 計(jì)算差異基因 （DESeq2差異基因分析和批次效應(yīng)移除）

將0、6、12 h的count的table依次導(dǎo)入,分別計(jì)算這3個(gè)時(shí)間點(diǎn)的差異基因。

過(guò)濾掉低表達(dá)的基因（所有樣本的和 >1），log2FC >1 with adjusted p-values<0.01 篩選差異基因

>library(DESeq2)

> colname(table)=c('0h_c_rep1','0h_c_rep2','0h_c_rep3','0h_t_rep1','0h_t_rep2','0h_t_rep3')

>condition <- c('c','c','c','t','t','t')

>dds <- DESeqDataSetFromMatrix(countData = countData, colData = group, design = ~condition )

>res <- na.exclude(as.data.frame(results))

> write.table(filter,'DGE_0h.txt', quote=F,sep='\t', col. names = NA)

6 使用EBSeq 對(duì)count數(shù)目進(jìn)行標(biāo)準(zhǔn)化

（生信寶典注：標(biāo)準(zhǔn)化后的數(shù)據(jù)可以也直接從DESeq2獲取）

獲得每個(gè)基因的counts數(shù)目之后，就可以計(jì)算兩兩基因的皮爾遜相關(guān)系數(shù)。使用EBSeq的函數(shù)median normalization 對(duì)count數(shù)目進(jìn)行標(biāo)準(zhǔn)化，再進(jìn)行對(duì)數(shù)化。

> if (!requireNamespace('BiocManager', quietly = TRUE))

BiocManager::install('EBSeq', version = '3.8')

> NormData <- GetNormalizedMat(counts, MedianNorm(counts))

> Norm.interest <- NormData.log[rownames(filter),] # filter是步驟5獲得的基因集

7 使用psych 計(jì)算兩兩相關(guān)性

psych在CRAN官網(wǎng)進(jìn)行下載。使用corr.test (也可以使用WGCNA的bicor函數(shù)，如果用Python更快)函數(shù)對(duì)差異基因進(jìn)行兩兩的相關(guān)性計(jì)算。如果用所有基因計(jì)算數(shù)據(jù)量將大得可怕。因此我們只挑選差異基因進(jìn)行分析。如我們對(duì)這個(gè)2 x 3的實(shí)驗(yàn)組進(jìn)行兩兩差異計(jì)算，得到一個(gè)差異基因的總表，挑選這些基因進(jìn)行計(jì)算。以相關(guān)性 > 0.9 p小于0.01篩選出共表達(dá)基因

>library('psych')

> Norm.interest.corr$p[lower.tri( Norm.interest.corr$p,diag=TRUE)]=NA

> Correlation <- as.data.frame(as.table(Norm.interest.corr $r))

> colnames(Cor.table) <- c('gene1','gene2','cor','p.adj')

8 使用igraph和Cytoscape可視化網(wǎng)絡(luò)

使用igrah對(duì)前面篩出的互作基因Cor.table.filt進(jìn)行網(wǎng)絡(luò)分析，degree 是指節(jié)點(diǎn) (這里指基因)的連接度，即一個(gè)點(diǎn)有多少條邊相連，degree centrality是某個(gè)節(jié)點(diǎn)的度除以網(wǎng)絡(luò)中所有點(diǎn)能構(gòu)成的連接數(shù)目，能反應(yīng)一個(gè)基因的中心度。介度中心性 (betweenness centrality) 是指一個(gè)節(jié)點(diǎn)充當(dāng)其它兩個(gè)節(jié)點(diǎn)中間人的次數(shù)“被經(jīng)過(guò)”的感覺，用“被經(jīng)過(guò)次數(shù)”除以總的ties，即n（n-1）/2，計(jì)算方法見下圖 （來(lái)源于 易生信陳亮博士的擴(kuò)增子和宏基因組授課，本期推文同期有陳亮博士的“一文學(xué)會(huì)網(wǎng)絡(luò)分析 igraph更詳細(xì)假設(shè)”）。輸出的文件Node_nw_st.txt和之前獲得的Cor.table.filt兩個(gè)表格共同用于Cytoscape可視化了。

# By Chen Liang
node_pro <- function(igraph){  
  # 節(jié)點(diǎn)度
  igraph.degree <- igraph::degree(igraph)
  # 節(jié)點(diǎn)度中心性
  igraph.cen.degree <- round(centralization.degree(igraph)$res,3)
  # 節(jié)點(diǎn)介數(shù)中心性
  igraph.betweenness <- round(centralization.betweenness(igraph)$res,3)
  # 節(jié)點(diǎn)中心性
  igraph.closeness <- round(centralization.closeness(igraph)$res,3)
  #igraph.node.pro <- cbind(igraph.degree,igraph.closeness,igraph.betweenness,igraph.cen.degree)
  #colnames(igraph.node.pro)<-c('igraph.degree','igraph.closeness','igraph.betweenness','igraph.cen.degree')
  igraph.node.pro <- data.frame(Node=names(igraph.degree), igraph.degree,igraph.closeness,igraph.betweenness,igraph.cen.degree)
}
net_node_attr <- node_pro(net)[,1:4]
write.table(net_node_attr, file='Node_nw_st.txt', row.names=F, col.names=T,quote=F,sep='\t')

 具體可視化見：

• network3D 交互式網(wǎng)絡(luò)生成

• Cytoscape教程1

• Cytoscape之操作界面介紹

• 新出爐的Cytoscape視頻教程

• Cytoscape制作帶bar圖和pie圖節(jié)點(diǎn)的網(wǎng)絡(luò)圖

• Cytoscape: MCODE增強(qiáng)包的網(wǎng)絡(luò)模塊化分析

• 一文學(xué)會(huì)網(wǎng)絡(luò)分析

R統(tǒng)計(jì)和作圖

• Graphpad，經(jīng)典繪圖工具初學(xué)初探

• 在R中贊揚(yáng)下努力工作的你，獎(jiǎng)勵(lì)一份CheatShet

• 別人的電子書，你的電子書，都在bookdown

• R語(yǔ)言 - 入門環(huán)境Rstudio

• R語(yǔ)言 - 熱圖繪制 (heatmap)

• R語(yǔ)言 - 基礎(chǔ)概念和矩陣操作

• R語(yǔ)言 - 熱圖簡(jiǎn)化

• R語(yǔ)言 - 熱圖美化

• R語(yǔ)言 - 線圖繪制

• R語(yǔ)言 - 線圖一步法

• R語(yǔ)言 - 箱線圖（小提琴圖、抖動(dòng)圖、區(qū)域散點(diǎn)圖）

• R語(yǔ)言 - 箱線圖一步法

• R語(yǔ)言 - 火山圖

• R語(yǔ)言 - 富集分析泡泡圖

• R語(yǔ)言 - 散點(diǎn)圖繪制

• R語(yǔ)言 - 韋恩圖

• R語(yǔ)言 - 柱狀圖

• R語(yǔ)言 - 圖形設(shè)置中英字體

• R語(yǔ)言 - 非參數(shù)法生存分析

• R語(yǔ)言 - 繪制seq logo圖

• WGCNA分析，簡(jiǎn)單全面的最新教程

• 一文看懂PCA主成分分析

• 富集分析DotPlot，可以服

• 基因共表達(dá)聚類分析和可視化

• R中1010個(gè)熱圖繪制方法

• 還在用PCA降維？快學(xué)學(xué)大牛最愛的t-SNE算法吧, 附Python/R代碼

• 一個(gè)函數(shù)抓取代謝組學(xué)權(quán)威數(shù)據(jù)庫(kù)HMDB的所有表格數(shù)據(jù)

• 文章用圖的修改和排版

• network3D: 交互式桑基圖

• network3D 交互式網(wǎng)絡(luò)生成

• Seq logo 在線繪制工具——Weblogo

• 生物AI插圖素材獲取和拼裝指導(dǎo)

• ggplot2高效實(shí)用指南 (可視化腳本、工具、套路、配色)

• 圖像處理R包magick學(xué)習(xí)筆記

• SOM基因表達(dá)聚類分析初探

• 利用gganimate可視化全球范圍R-Ladies（R社區(qū)性別多樣性組織）發(fā)展情況

• 一分鐘繪制磷脂雙分子層：AI零基礎(chǔ)入門和基本圖形繪制

• AI科研繪圖（二）：模式圖的基本畫法

• 你知道R中的賦值符號(hào)箭頭(<-)和等號(hào)（=）的區(qū)別嗎？

• R語(yǔ)言可視化學(xué)習(xí)筆記之ggridges包

• 利用ComplexHeatmap繪制熱圖(一)

• ggplot2學(xué)習(xí)筆記之圖形排列

• R包reshape2，輕松實(shí)現(xiàn)長(zhǎng)、寬數(shù)據(jù)表格轉(zhuǎn)換

• 用R在地圖上繪制網(wǎng)絡(luò)圖的三種方法

• PCA主成分分析實(shí)戰(zhàn)和可視化 附R代碼和測(cè)試數(shù)據(jù)

• iTOL快速繪制顏值最高的進(jìn)化樹！

• 12個(gè)ggplot2擴(kuò)展包幫你實(shí)現(xiàn)更強(qiáng)大的可視化

• 編程模板-R語(yǔ)言腳本寫作：最簡(jiǎn)單的統(tǒng)計(jì)與繪圖，包安裝、命令行參數(shù)解析、文件讀取、表格和矢量圖輸出

• R語(yǔ)言統(tǒng)計(jì)入門課程推薦——生物科學(xué)中的數(shù)據(jù)分析Data Analysis for the Life Sciences

• 數(shù)據(jù)可視化基本套路總結(jié)

• 你知道R中的賦值符號(hào)箭頭<-和等號(hào)=的區(qū)別嗎？

• 使用dplyr進(jìn)行數(shù)據(jù)操作30例

• 交集intersect、并集union、找不同setdiff

• R包reshape2，輕松實(shí)現(xiàn)長(zhǎng)、寬數(shù)據(jù)表格轉(zhuǎn)換

• 1數(shù)據(jù)類型（向量、數(shù)組、矩陣、 列表和數(shù)據(jù)框）

• 2讀寫數(shù)據(jù)所需的主要函數(shù)、與外部環(huán)境交互

• 3數(shù)據(jù)篩選——提取對(duì)象的子集

• 4向量、矩陣的數(shù)學(xué)運(yùn)算

• 5控制結(jié)構(gòu)

• 6函數(shù)及作用域

• 7認(rèn)識(shí)循環(huán)函數(shù)lapply和sapply

• 8分解數(shù)據(jù)框split和查看對(duì)象str

• 9模擬—隨機(jī)數(shù)、抽樣、線性模型

• 1初識(shí)ggplot2繪制幾何對(duì)象

• 2圖層的使用—基礎(chǔ)、加標(biāo)簽、注釋

• 3工具箱—誤差線、加權(quán)數(shù)、展示數(shù)據(jù)分布

• 4語(yǔ)法基礎(chǔ)

• 5通過(guò)圖層構(gòu)建圖像

• 6標(biāo)度、軸和圖例

• 7定位-分面和坐標(biāo)系

• 8主題設(shè)置、存儲(chǔ)導(dǎo)出

• 9繪圖需要的數(shù)據(jù)整理技術(shù)

• 創(chuàng)建屬于自己的調(diào)色板

• 28個(gè)實(shí)用繪圖包，總有幾個(gè)適合你

• 熱圖繪制

• R做線性回歸

• 繪圖相關(guān)系數(shù)矩陣corrplot

• 相關(guān)矩陣可視化ggcorrplot

• 繪制交互式圖形recharts

• 交互式可視化CanvasXpress

• 聚類分析factoextra

• LDA分析、作圖及添加置信-ggord

• 解決散點(diǎn)圖樣品標(biāo)簽重疊ggrepel

• 添加P值或顯著性標(biāo)記ggpubr

• Alpha多樣性稀釋曲線rarefraction curve

• 堆疊柱狀圖各成分連線畫法：突出組間變化

• 沖擊圖展示組間時(shí)間序列變化ggalluvial

• ?；鶊Driverplot

• 微生物環(huán)境因子分析ggvegan

• 五彩進(jìn)化樹與熱圖更配ggtree

• 多元回歸樹分析mvpart

• 隨機(jī)森林randomForest 分類Classification 回歸Regression

• 加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析WGCNA

• circlize包繪制circos-plot

• R語(yǔ)言搭建炫酷的線上博客系統(tǒng)

易生信系列培訓(xùn)課程，掃碼獲取免費(fèi)資料

更多閱讀

畫圖三字經(jīng) 生信視頻 生信系列教程 

心得體會(huì) TCGA數(shù)據(jù)庫(kù) Linux Python 

高通量分析 免費(fèi)在線畫圖 測(cè)序歷史 超級(jí)增強(qiáng)子

生信學(xué)習(xí)視頻 PPT EXCEL 文章寫作 ggplot2

海哥組學(xué) 可視化套路 基因組瀏覽器

色彩搭配 圖形排版 互作網(wǎng)絡(luò)

自學(xué)生信