目的：
對illumina數據進行處理，利用 RNA-Seq 發現新的 RNA 變體和剪接位點，或量化 mRNA 以進行基因表達分析等。對兩組或多組樣本的轉錄組數據，通過差異表達分析和對所發現的差異表達基因集合進行功能富集分析以推斷生物學功能。
數據準備：
數據下載：
- Human genome(GRCh38/hg3)：Index of /goldenPath/hg38/chromosomes (ucsc.edu)
- Mouse genome (GRCm39/mm39)：Index of /goldenPath/mm39/bigZips (ucsc.edu)
- D. melanogaster genome(BDGP Release 6 + ISO1 MT/dm6)：Index of /goldenPath/dm6/bigZips (ucsc.edu)
2）註釋文件（.gtf/gff）下載
- Human genome Annotation: Index of /pub/release-104/gtf/homo_sapiens/ (ensembl.org)
數據質控：
可採用Trim Galore或fastqc進行質控：
Fastqc軟件參數：
#示例：fastqc -o qc/ *.fastq.gz
# 基本格式
# fastqc [-o output dir] [--(no)extract] [-f fastq|bam|sam] [-c contaminant file] seqfile1 .. seqfileN
# 主要是包括前面的各種選項和最後面的可以加入N個文件
# -o --outdir FastQC生成的報告文件的儲存路徑，生成的報告的文件名是根據輸入來定的
# --extract 生成的報告默認會打包成1個壓縮文件，使用這個參數是讓程序不打包
# -t --threads 選擇程序運行的線程數，每個線程會佔用250MB內存，越多越快咯
# -c --contaminants 污染物選項，輸入的是一個文件，格式是Name [Tab] Sequence，裏面是可能的污染序列，如果有這個選項，FastQC會在計算時候評估污染的情況，並在統計的時候進行分析，一般用不到
# -a --adapters 也是輸入一個文件，文件的格式Name [Tab] Sequence，儲存的是測序的adpater序列信息，如果不輸入，目前版本的FastQC就按照通用引物來評估序列時候有adapter的殘留
# -q --quiet 安靜運行模式，一般不選這個選項的時候，程序會實時報告運行的狀況。
Trim Galore軟件參數：
trim_galore -q 20 --phred33 --stringency 3 --length 20 -e 0.1 \
--paired $dir/cmp/01raw_data/$fq1 $dir/cmp/01raw_data/$fq2 \
--gzip -o $input_data
--hardtrim5：用於從序列的3’端切除鹼基，示意如下
before: 5’-CCTAAGGAAACAAGTACACTCCACACATGCATA-3’
--hardtrim5 20: CCTAAGGAAACAAGTACACT
通過hardtrim5參數可以將序列截取成固定長度。
--hardtrim3參數：從序列的5’端切除鹼基。
--quality：設定Phred quality score閾值，默認為20。
--phred33：選擇-phred33或者-phred64，表示測序平台使用的Phred quality score。
--adapter：輸入adapter序列。也可以不輸入，Trim Galore!會自動尋找可能性最高的平台對應的adapter。自動搜選的平台三個，也直接顯式輸入這三種平台，即--illumina、--nextera和--small_rna。
--stringency：設定可以忍受的前後adapter重疊的鹼基數，默認為1。可以適度放寬，因為後一個adapter幾乎不可能被測序儀讀到。
--length：設定輸出reads長度閾值，小於設定值會被拋棄。
--paired：對於雙端測序結果，一對reads中，如果有一個被剔除，那麼另一個會被同樣拋棄，而不管是否達到標準。
--retain_unpaired：對於雙端測序結果，一對reads中，如果一個read達到標準，但是對應的另一個要被拋棄，達到標準的read會被單獨保存為一個文件。
--gzip和--dont_gzip：清洗後的數據zip打包或者不打包。
--output_dir：輸入目錄。需要提前建立目錄，否則運行會報錯。
-- trim-n : 移除read一端的reads。
一般步驟：
- 去除read 3‘端低質量鹼基
- 去除adapter(接頭)序列
有三種常見接頭：
Illumina: AGATCGGAAGAGC
Small RNA: TGGAATTCTCGG
Nextera: CTGTCTCTTATA
- 去除很短的序列
- 其他
基本步驟：
Reads align
$ hisat2 -p 8 --dta -x chrX_data/indexes/chrX_tran -1
chrX_data/samples/ERR188044_chrX_1.fastq.gz -2
chrX_data/samples/ERR188044_chrX_2.fastq.gz -S
ERR188044_chrX.sam
Sort sam to bam
samtools sort -@ 8 -o ERR188044_chrX.bam ERR188044_chrX.sam
Assemble transcripts
$ stringtie -p 8 -G chrX_data/genes/chrX.gtf -o
ERR188044_chrX.gtf –l ERR188044 ERR188044_chrX.bam
Merge transcripts
stringtie --merge -p 8 -G chrX_data/genes/chrX.gtf -o
stringtie_merged.gtf chrX_data/mergelist.txt
Examine the transcripts compare to reference annotation
$ gffcompare –r chrX_data/genes/chrX.gtf –G –o merged
stringtie_merged.gtf
Estimate transcript abundances ,create table counts for Ballgown:
$ stringtie –e –B -p 8 -G stringtie_merged.gtf -o
ballgown/ERR188044/ERR188044_chrX.gtf ERR188044_chrX.bam
differential expression analysis
R version 4.2.0 (2022-04-22 ucrt) -- "Vigorous Calisthenics"
Copyright (C) 2022 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)#以GO富集分析為例
rm(list = ls())
#load data
DEGdir = "D:/Learn_Summary/Bioinformat/RNA_seq/R_plot/DEG"
setwd(DEGdir)
info <- read.csv("Deseq2_HC_vs_Treat_diff.addExpression.xls",
header=TRUE,sep="\t",fill=TRUE)[,1:4,1]
info$sig = as.factor(ifelse(info$pvalue < 0.05 & abs(info$log2FoldChange) >= 1.6,
ifelse(info$log2FoldChange> 1.6 ,'Up','Down'),'NoSignifi'))
up <- subset(info, sig == 'Up')
up <- up[order(up$pvalue), ][1:70, ] #取前70個
down <- subset(info, sig == 'Down')
down <- down[order(down$pvalue), ][1:70, ] #取前70個
info <- rbind(up,down)
output.gene_id = data.frame(gene_id = info$GeneName)
write.table(output.gene_id,file="gene_name.txt",col.names = F,row.names = F,sep = "\t",quote = F)
#cuffdiff_result = read.table(file="./hela_gene_exp.diff",header = T,sep = "\t")
#cuffdiff_result$sample_1 = "treat"
#cuffdiff_result$sample_2 = "ctrl"
#dim(cuffdiff_result)
dim(info)
# select_vector = (info$value_1 > 1 | cuffdiff_result$value_2 > 1) & (abs(infot$log2.fold_change.) >= 1.5) & (info$p_value < 0.05)
# cuffdiff_result.sign = cuffdiff_result[select_vector,]
# dim(cuffdiff_result.sign)
# 安裝包
# source("https://bioconductor.org/biocLite.R")
# BiocManager::install("clusterProfiler") #用來做富集分析
# BiocManager::install("topGO") #畫GO圖用的
# BiocManager::install("Rgraphviz")
# BiocManager::install("pathview") #看KEGG pathway的
# BiocManager::install("org.Hs.eg.db") #這個包裏存有人的註釋文件
# install.packages("R.utils")
# 載入包
library(clusterProfiler)
library(topGO)
library(Rgraphviz)
library(pathview)
library(org.Hs.eg.db)
# library()
DEG$sig = as.factor(ifelse(DEG$pvalue < 0.05 & abs(DEG$log2FoldChange) >= 1.6,
ifelse(DEG$log2FoldChange> 1.6 ,'Up','Down'),'NoSignifi'))
up <- subset(DEG, sig == 'Up')
up <- up[order(up$pvalue), ][1:70, ]
down <- subset(DEG, sig == 'Down')
down <- down[order(down$pvalue), ][1:70, ]
#提取symbol ID --> 轉換為ENTREZID
# info <- read.csv("Deseq2_HC_vs_Treat_diff.addExpression.xls",header=TRUE,sep="\t",fill=TRUE)
DEG.gene_symbol = as.character(output.gene_id$gene_id) #獲得基因 symbol ID
DEG.entrez_id = mapIds(x = org.Hs.eg.db,
keys = DEG.gene_symbol,
keytype = "SYMBOL",
column = "ENTREZID")
DEG.entrez_id <- na.omit(DEG.entrez_id) #除去轉換失敗的NA
#BP（Biological process）層面上的富集分析
erich.go.BP = enrichGO(gene = DEG.entrez_id,
OrgDb = org.Hs.eg.db,
keyType = "ENTREZID",#設定讀取的gene_ID類型
ont = "BP",
pvalueCutoff = 0.5,#設定p值閾值
qvalueCutoff = 0.5,#設定q值閾值
readable = T)
write.table(as.data.frame(erich.go.BP@result), file="GO_BP.xls", sep="\t", row.names=F)
##分析完成後，作圖
setwd(paste(DEGdir,"/plot",sep = ""))
pdf(file = "GO_BP_dotplot.pdf",width = 10,height = 15)
dotplot(erich.go.BP)
dev.off()
pdf(file = "GO_BP_barplot.pdf",width = 15,height = 20)
barplot(erich.go.BP)
dev.off()
#CC（Cellular component）層面上的富集分析
erich.go.CC = enrichGO(gene = DEG.entrez_id,
OrgDb = org.Hs.eg.db,
keyType = "ENTREZID",
ont = "CC",
pvalueCutoff = 0.5,
qvalueCutoff = 0.5)
write.table(as.data.frame(erich.go.CC@result), file="GO_CC.xls", sep="\t", row.names=F)
## 畫圖
pdf(file = "GO_CC_barplot.pdf",width = 15,height = 20)
barplot(erich.go.CC)
dev.off()
pdf(file = "GO_CC_dotplot.pdf",width = 15,height = 20)
dotplot(erich.go.CC)
dev.off()
#MF（Mollecular Function）層面上的富集分析
erich.go.MF = enrichGO(gene = DEG.entrez_id,
OrgDb = org.Hs.eg.db,
keyType = "ENTREZID",
ont = "MF",
pvalueCutoff = 0.5,
qvalueCutoff = 0.5)
write.table(as.data.frame(erich.go.MF@result), file="GO_MF.xls", sep="\t", row.names=F)
## 畫圖
pdf(file = "GO_MF_barplot.pdf",width = 15,height = 20)
barplot(erich.go.MF)
dev.off()
pdf(file = "GO_MF_dotplot.pdf",width = 15,height = 20)
dotplot(erich.go.MF)
dev.off()
#all_GO分析
GO = enrichGO(gene = DEG.entrez_id,
OrgDb = org.Hs.eg.db,
keyType = "ENTREZID",
ont = "ALL",
pvalueCutoff = 0.5,
qvalueCutoff = 0.5,
readable =T)
## 畫圖
# pdf(file = "GO_barplot.pdf",width = 15,height = 20)
# barplot(GO)
# dev.off()
# pdf(file = "GO_dotplot.pdf",width = 15,height = 20)
# dotplot(GO)
# dev.off()
pdf(file = "GO_all_barplot.pdf",height = 10,width = 9)
barplot(GO, split="ONTOLOGY",title = "GO Pathway",font.size = 9)#柱狀圖
dev.off()
# pdf(file = "KEGG_all_barplot.pdf",height = 35,width = 35)
# barplot(erich.kegg,showCategory = 40,title = 'KEGG Pathway',font.size = 10)
# dev.off()
pdf(file = "GO_all_dotplot.pdf",height = 10,width = 9)
dotplot(GO, split="ONTOLOGY",title = "GO Pathway",
font.size=9,size = NULL)#點狀圖
dev.off()
# pdf(file = "KEGG_all_dotplot.pdf",height = 35,width = 35)
# dotplot(erich.kegg)
# dev.off()
#轉換為樹圖
pdf(file="GO_BP.tree.pdf",width = 10,height = 12)
# par(cex=0.1)
plotGOgraph(erich.go.BP)+
theme(text=element_text(size=20))
dev.off()
# cex_label_category=2
pdf(file = "GO_CC.tree.pdf",width = 15,height = 20)
plotGOgraph(erich.go.CC)
dev.off()
pdf(file = "GO_MF.tree.pdf",width = 15,height = 20)
plotGOgraph(erich.go.MF)
dev.off()#保存為pdf
基本分析：
有參考基因組，需要組裝新的轉錄本：
HISAT2 -> StringTie -> Ballgown(或)htseq-count + DESeq2
注：樣本量如果大於10，推薦使用TACO，代替cuffmerge或stringtie
HISAT2的使用：
示例：
$ hisat2 -p 8 --dta -x chrX_data/indexes/chrX_tran -1
chrX_data/samples/ERR188044_chrX_1.fastq.gz -2
chrX_data/samples/ERR188044_chrX_2.fastq.gz -S
ERR188044_chrX.sam
主要參數：
-x <hisat2-idx>
參考基因組索引文件的前綴（目錄及文件前綴）。
-1 <m1>
雙端測序結果的第一個文件。若有多組數據，使用逗號將文件分隔。Reads的長度可以不一致。
-2 <m2>
雙端測序結果的第二個文件。若有多組數據，使用逗號將文件分隔，並且文件順序要和-1參數對應。Reads的長度可以不一致。
-U <r>
單端數據文件。若有多組數據，使用逗號將文件分隔。可以和-1、-2參數同時使用。Reads的長度可以不一致。
–sra-acc <SRA accession number>
輸入SRA登錄號，比如SRR353653，SRR353654。多組數據之間使用逗號分隔。HISAT將自動下載並識別數據類型，進行比對。
-S <hit>
指定輸出的SAM文件（目錄及名稱）。
-t/--time print wall-clock time taken by search phases
-p/--threads <int> number of alignment threads to launch 指定進程數目（8/16）
Samtools的使用
示例：
#version 1.3 or newer
samtools sort -@ 8 -o ERR188044_chrX.bam ERR188044_chrX.sam
#version 1.2 or older
samtools view -bS ERR188044_chrX.sam > ERR188044_chrX_unsorted.sam
samtools sort -@ 8 ERR188044_chrX_unsorted.bam ERR188044_chrX
主要參數：
參考：Samtools使用大全
Stringtie的使用
主要參數：
|
參數 |
描述 |
|
-G <ref_ann.gff> |
使用註釋好的gtf文件輔助組裝，在-e未設置的條件下，輸出中包括註釋文件中的轉錄本和預測的新型轉錄本 |
|
-o [<path/>]<out.gtf> |
輸出文件的名字，最好是全路徑，默認輸出為標準輸出 |
|
-l <label> |
為輸出的轉錄本設置前綴名，默認為STRG |
|
-p <int> |
線程數，默認為1 |
|
-A <gene_abund.tab> |
對輸出的gtf統計基因表達量，並以一個tab分割的文件輸出，這裏需要提交輸出的文件名 |
|
-C <cov_refs.gtf> |
對輸出的gtf中屬於-G提交的參考gtf的轉錄本統一輸出到該文件，這裏需要提交一個文件名 |
|
-B |
是否需要輸出Ballgown可以識別的文件,在-b設置的情況下，使用-o的路徑輸出 |
|
-b <path> |
對Ballgown輸出的文件指定一個路徑保存 |
|
-e |
我認為是最需要注意的參數！！只統計可以匹配-G提交的參考gtf中的轉錄本，不再對新的轉錄本做預測，這可以加快程序的運行速度 |
|
-m <int> |
對預測的轉錄本設置最小長度，默認為200 |
參考：StringTie的使用説明
Ballgown的使用
參考：Ballgown.pdf
有參考基因組，不需組裝新的轉錄本：
參考：https://github.com/hnnd/NGS/blob/master/pip_b.md
無參考基因組：
trinity > rsem > DEseq2
參考：https://github.com/hnnd/NGS/blob/master/pip_c.md
其他分析：
- 可變剪切分析
- 融合基因分析
- GO富集分析
- Pathway富集分析
基本分析流程[1]：
Concept
Abundance estimate：表現為FPKM值，即每百萬映射讀取的轉錄片段每千鹼基。
Total：clean data雙端總reads數目Duplicate：重複的reads數目。
Mapped：比對到參考基因組上的總reads數目（比例）。
Properly mapped：比對到參考基因組且方向正確的reads數目（比例）。
PE mapped：雙端reads比對到參考基因組上的reads數目（比例）。
SE mapped：僅單端read比對到參考基因組上的reads數目（比例）。
with mate mapped to a different chr：比對到不同染色體的reads數目。
with mate mapped to a different chr (mapQ>=5)：比對到不同染色體且比對質量不低於5的reads數目。
Average_sequencing_depth:比對到參考基因組的平均測序深度（測序數據量/基因組大小）。
Coverage：比對數據對全基因組區域的覆蓋度（鹼基覆蓋長度佔全基因組鹼基總長的比例）。
Coverage_at_least_4X:全基因組區域中鹼基覆蓋深度不低於4X的比例。
Coverage_at_least_10X：全基因組區域中鹼基覆蓋深度不低於10X的比例。
Coverage_at_least_20X：全基因組區域中鹼基覆蓋深度不低於20X的比例。
RNA-seq可行的分析流程與方向
全轉錄組是指特定細胞在特定狀態下所能轉錄出來的所有RNA的總和，包括mRNA和非編碼RNA（non-coding RNA），針對非編碼RNA的研究主要集中在具有調控作用的miRNA，lncRNA和circRNA。
基於高通量測序的RNA測序（RNA-Seq）是一種高度靈敏和準確的檢測全轉錄組表達的工具，研究人員可以在單次分析中同時分析同一樣本中的mRNA，miRNA，lncRNA，circRNA，檢測已知的和發現新的特徵，使轉錄本亞型、基因融合、單核苷酸位點變異和其他特徵不受先驗知識的限制。研究人員能夠檢測疾病狀態、治療反應、各種環境條件下以及一系列其他研究設計中以前無法檢測的變化。全面快速地獲得RNA序列和表達丰度信息，通過關聯分析，還可深入挖掘生命現象背後的轉錄調控問題。
靶向RNA測序則是對感興趣的轉錄本表達進行高度準確且特異測定的方法，同時提供定量和定性信息，能實現表達分析、融合基因鑑定。通過對特定基因的表達譜進行分析，可以評估疾病相關的變異和表觀遺傳學改變。基因融合和基因表達變化分析為癌症功能相關的改變提供了特定的視角。
Small RNA主要包括miRNA、piRNA、tsRNA（tRF&tiRNA）、snRNA和snoRNA等，是一類不具有蛋白編碼能力的RNA分子，能調控基因表達，在細胞生長、發育和代謝等基礎生物學過程中都扮演着重要的角色，甚至在癌症等相關疾病形成過程中也起着關鍵的作用。高通量測序技術省去了煩瑣的Small RNA克隆文庫構建過程，可以一次性生成上百萬條Small RNA序列，無需預先假設，能夠快速鑑定特定條件下表達的已知Small RNA並發現新的Small RNA，以單鹼基分辨率鑑定變異，同時還可以研究不同條件下Small RNA的表達差異，並可以聯合轉錄組測序表達譜數據進行關聯分析，在轉錄和轉錄後水平研究基因調控。
基因結構水平分析：與參考基因組比對、新轉錄本預測、可變剪切分析、SNP/InDel註釋分析、融合基因分析、癌基因註釋等；基因表達水平分析：RNA-seq整體質量評估、差異基因分析、GO/KEGG富集分析、蛋白互作網絡分析、轉錄因子註釋等。
1、原始數據質控
以原始數據為研究對象，對於低質量序列，未檢測序列，接頭序列進行過濾，並對於過濾前後數據的鹼基質量、GC含量、長度分佈、接頭留存和Duplication比率等指標進行分析。
2、RNA基因組比對（RNA Mapping）
採用Hisat2/Mapsplice/Star/Tophat2等算法進行基因組比對，得到基因組比對的bam文件，並基於bam文件進行信息統計，得到基因組比對率、reads在基因結構和染色體上的分佈結果。
3、表達量統計（Expression）
採用HTSeq以及基因組註釋的gff3文件，根據單端或雙端測序類型，選擇RPKM或FPKM的標化方式對基因表達量進行統計。基於統計結果，分析得到樣本間相關性、 RPKM/FPKM密度和豐度等分析結果，反映單個樣本基因表達水平分佈和離散程度，以及不同樣本整體基因表達水平的差異。
4、差異基因篩選（Dif Gene Analysis）
採用DESeq2/DESeq/EBSeq/EdgeR/Limma等算法進行差異篩選，得到滿足差異倍數（Fold Change）以及FDR閾值的差異基因，並基於差異篩選結果以及樣本的FPKM或RPKM，進行火山圖分析（Volcano Plot）以及聚類圖分析（Heatmap）。
5、功能分析（GO Analysis）
採用NCBI/UNIPROT/SWISSPROT/AMIGO等GO數據庫，對於差異基因進行功能分析，從而得到差異基因所顯著性富集的功能條目（GO Term）。
6、信號通路分析（Pathway Analysis）
整合KEGG、NCBI、EMBL等數據庫，深入優化所需算法，對差異基因進行信號通路分析，得到差異基因所顯著性富集的信號通路條目。
7、GO-Tree分析
採用GO數據庫中GO-term的上下級層級從屬關係，進行GO-Tree繪製，得到顯著性差異功能的功能簇以及層級從屬關係。
8、Path-Act-Network分析
採用KEGG數據庫記載的信號通路上下游調控關係，進行Path-Act-Network繪製，得到宏觀上的顯著性信號通路的上下游調控關係。
9、共表達網絡分析（Co-Exp-Network Analysis）高級分析
基於GO Analysis和Pathway Analysis得到的顯著性條目，以及研究者感興趣條目，以這些條目中基因的表達值為研究目標，進行共表達網絡和K-Core分析，從而得到基因間的相關性和基因的核心度，再以Co-Expression.txt和K-Core為研究對象，採用Cytoscape進行圖形化展示，得到Co-Expression-Network。
10、基因間相互作用關係網絡(Gene-Act-Network Analysis) 高級分析
基於GO Analysis和Pathway Analysis得到的顯著性條目，以研究者感興趣的相關表型基因為研究對象，採用KEGG數據庫基因間關係註釋，幫助研究者繪製Gene-Act-Network，快速定位“核心”基因。
11、韋恩分析
以各分組間的基因為研究對象，採用韋恩作圖分析的方法，可找出各分組間共有或者特有的差異表達基因並進行深入分析。
12、趨勢分析
以各差異分組間的韋恩基因的FPKM值為研究對象，採用STEM算法，進行趨勢分析，得到按照樣本邏輯順序所在趨勢。
13、加權基因共表達網絡分析（WGCNA）分析
基於加權的表達相關性，進行層級聚類分析，並根據設定標準切分聚類結果，獲得不同的基因模塊，採用聚類樹的分枝和不同顏色來鑑定高度協同變化的基因集。如果有表型信息，還可以計算基因模塊與表型相關性，鑑定性狀相關的模塊，並根據基因集的內連性和基因集與表型之間的關聯鑑定候補生物標記基因或治療靶點。
關於GO富集時的pathway選擇
不要僅僅基於Pathway富集分析的結果解讀數據，人為的解讀和挑選是必不可少的。因為生物數據的解讀，在現階段更多是生物學問題，而不是數學問題。原因大體如下：
（1）基因調控是個系統，不要僅僅看成1個1個孤立的pathway。
“在植物體內其實根本就不存在pathway，什麼脱落酸通路，水楊酸通路，其實這些調控因子都是相互聯通，相互影響的，是個整體。只是我們人類為了研究方便，人為將這些系統拆分各個子集。 ” 所以，如果你真的將pathway看成1個個破碎的途徑，以為某種處理只會影響某個pathway，富集分析必須在數學上或統計學上得到1個指向性很強的結論，那是不大可能的。
具體説了，説基因調控是個系統，可以從兩個層面進行解讀：
a）1個基因的改變可以造成整個系統的改變；
舉幾個例子：
把1個生命活動必須的蛋白敲除後，整個細胞會發生紊亂。而植物抗病應激，也往往是1個受體蛋白識別了病原的外源蛋白，然後導致整個細胞系統的變化。
b）1個基因往往有多個功能，但執行具體的功能往往是不同蛋白複合物共同作用的結果。
例如。基因X理論上在不同情況下，有可能參與A、B、C通路。在某個生物處理下，或許基因X 只在A通路里起作用。但如果進行基因註釋的話，X同樣也會被註釋到B、C。所以，富集分析的結果總是會涉及特別多的通路。例如，研究人的項目，無論什麼研究背景，常常會富集到 帕金森綜合症通路。不是你的材料真的得了帕金森綜合症，只是那些與你實驗處理相關的基因，在一定條件下也可以參與到帕金森綜合症的過程，所以被註釋到了這個通路里。
小結：所以，我們也看到了。無論什麼實驗處理，總有可能導致整個系統的變化。同時，基因的通路註釋也有欺騙性。那麼，從這一堆冗餘信息中，想得到與我們研究相關的結論，離不開人為的篩選也解讀。從那個複雜的整體中，篩選出核心的局部片段，這是個技術活。“這樣的話是否存在一個問題就是在結果的解釋上比較主觀，也會因自身背景知識的不足而漏掉一些新穎的結果”。那當然，同樣的結果1個外行可能什麼都沒有看見。但1個資深的學者可能會把握到很精彩的內容。好像任何領域都是如此，除了提高內功好像沒有其他捷徑。
（2）pathway富集分析的統計假設，並非在任何情況下都適用
pathway富集分析，在生物學上的假設是：1個pathway 上游基因的改變，會導致下游相關基因改變，從而改變通路中大量基因的表達，達到統計學上富集的效果。但很多pathway中，基因A、B、C並不是相互調控的關係，而是共同參與某個過程的不同部分。
例如，代謝物X的合成修飾。基因A、B、C分步驟參與合成的3個步驟。基因A 給X前體加了 羥基，然後傳遞到下游；基因B又給X前體加了苯環，再傳遞到下游；基因C又給X的前體 加了個乙酰基 ，完成X的合成。那麼，基因A、B、C是參與了的相同的通路。如果基因A發生表達量變化，會直接調控影響B、C的表達量變化嗎？ 看來很有可能不會，所以從RNA-seq差異分析的富集分析結果中，這個通路是不顯著的。那麼基因A的表達變化是否有生物學意義？ 當然有，因為代謝物X的合成的確受影響了。
類似的例子，理論上DNA差異甲基化的結果，就不能看pathway富集分析的結果。1個pathway 1個基因的DNA甲基化變化，就足以改變這個通路的基因表達，而不需要整個通路的甲基化都發生變化。DNA甲基化、組蛋白CHIP-seq的結果，其實只看功能註釋、或通路註釋就足夠了，不需要考慮富集。
所以，我們還是要觀察、理解某個核心pathway中基因的相互作用，才能判斷其中的基因變化是否有生物學意義， 而不僅僅看富集分析的p value或Q value。
（3）目前的pathway是不完整的。
目前KEGG等數據庫，收錄的是已有的研究結果。但這些pathway的信息，遠沒有到達完善的水準。大部分通路只是瞭解1個大概的調控途徑，而中間有什麼轉錄因子參與、是否還有其他代謝物的生成，都是不知道的。這些通路的完整性，也會影響pathway富集分析結果。例如，基因A發生變化了，看起來下游基因沒有變化。也許是還有其他的調控在起作用，只是這些調控作用現在還不知道而已。
總結：pathway 和 GO富集分析結果的解讀，應該從生物學意義的角度出發，P value 和 Q value只是個參考而已，那些不顯著的通路也值得解讀（從功能註釋的角度解讀，而不是從富集分析的角度解讀）。只要結果可以解釋，有意義，不用太迷信P value。
If you want these code to practice it ，please contact me：zengcj2035@163.com
Here are some relevant information：
- https://github.com/xiangsui5/RNA-seq
- https://files.cnblogs.com/files/blogs/698237/expriment_report.rar?t=1697439634&download=true
