什么是基因表達,如下是來自于維基百科的解釋: Gene expression is the process by which information from a gene is used in the synthesis of a functional gene product. These products are often proteins, but in non-protein coding genes such as transfer RNA (tRNA)or small nuclear RNA (snRNA)genes, the product is a functional RNA. The process of gene expression is used by all known life—eukaryotes (including multicellular organisms),prokaryotes (bacteria and archaea), and utilized by viruses—to generate the macromolecular machinery for life. Flow of genetic information 研究方法 定量PCR 這部分我不太懂,所以就放幾段百度百科和維基百科的定義。 百度百科 定量PCR(即時聚合酶鏈鎖反應,Real-time Polymerase Chain Reaction,簡稱Real-time PCR、即時PCR), 又稱定量即時聚合酶鏈鎖反應(Quantitative real time polymerase chain reaction,簡稱 Q-PCR/qPCR/rt-qPCR、定量即時PCR、即時定量PCR),是一種在DNA擴增反應中,以螢光染劑偵測每次聚合酶鏈鎖反應(PCR)循環(huán)后產(chǎn)物總量的方法技術,有廣義概念和狹義概念。 廣義概念的定量PCR技術是指以外參或內(nèi)參為標準,通過對PCR終產(chǎn)物的分析或PCR過程的監(jiān)測,進行PCR起始模板量的定量。 狹義概念的定量PCR技術(嚴格意義的定量PCR技術)是指用外標法(熒光雜交探針保證特異性)通過監(jiān)測PCR過程(監(jiān)測擴增效率)達到精確定量起始模板數(shù)的目的,同時以內(nèi)對照有效排除假陰性結(jié)果(擴增效率為零)。 維基百科 A real-time polymerase chain reaction (Real-Time PCR), also known as quantitative polymerase chain reaction (qPCR), is a laboratory technique of molecular biology based on the polymerase chain reaction (PCR). It monitors the amplification of a targeted DNA molecule during the PCR, i.e. in real-time, and not at its end, as in conventional PCR. Real-time PCR can be used quantitatively (quantitative real-time PCR), and semi-quantitatively, i.e. above/below a certain amount of DNA molecules (semi quantitative real-time PCR). 優(yōu)點:靈敏性高,準確性高,通量也還行。 一般而言,RNA-Seq和microassay分析得到的差異表達基因最終也需要通過這種實驗方法進行驗證。 但是一般適用于驗證實驗,而不是用于探索性實驗。 microarray<small>基因矩陣</small> 基因芯片的概念在上個世紀80年代就已經(jīng)提出來了, 被評為1998年度自然科學領域十大進展之一。他的基本原理通過設計專門的短核苷酸作為探針,把這些探針固定在專門的基片表面,然后用樣本的cDNA進行雜交,根據(jù)雜交信號的強弱來判斷基因表達的程序。 microarray 但是microarray檢測的基因數(shù)量完全取決于你的探針設計的數(shù)量,而且難以研究mRNA的可變剪切。 RNA-Seq RNA-Seq是目前基因表達分析最常用的技術。分為以下幾步 1.分離所有mRNA 2.逆轉(zhuǎn)錄mRNA成cDNA 3.對cDNA測序 4.比對參考基因組 RNA-Seq實驗設計中的“重復”包括:技術重復和生物學重復 (1)技術重復為了估計測量技術(RNA-Seq)的變異。 (2)生物學重復是為了發(fā)現(xiàn)生物組內(nèi)的變異。 RNA-Seq的概率分布 (注:這與相同細胞不同基因表達的分布不同) 但是大多數(shù)基因表達實驗都是用一群細胞,幾乎沒有相應分布提出。 RNA-Seq試驗中,抽樣得到的raw read counts服從泊松分布。并且同一樣本在兩次試驗中的結(jié)果不同,這稱為shot noise。這種變異在RNA-Seq技術重復間成為Possion noise。 生物學上不同的樣本間的差異服從負二項(negative binomial)分布,有時稱gamma-Poisson分布。 由于RNA-Seq count數(shù)據(jù)也表現(xiàn)出zero inflation(大量值為0)的特征,所以很難擬合到負二項分布,所以有文章認為要用Poisson-Tweedie family建模。 RNA-seq數(shù)據(jù)和microassay在差異表達分析上的區(qū)別: 1.RNA-Seq觀察到的數(shù)據(jù)是抽樣過程中產(chǎn)生的離散(discrete)count形式。 也就是說總體是恒定的,表達量越高的基因在抽樣結(jié)果中所占的比例越大。 表達量低的基因可能即便有也無法被檢測出來。 當然,重新對相同文庫進行測序,還是有可能找到更多表達的轉(zhuǎn)錄本 2.microassay檢測的是熒光信號的連續(xù)度量。 由于使用固定的核酸序列去雜交,所以不是一種“零和游戲”,只要能雜交,就能被檢測。 (但如果沒有設計相應的引物,就不能檢測到可能的基因) 研究意義 1.在不同背景下比較mRNA水平 ①同一物種,不同組織: 研究基因在不同部分的表達情況 ②同一物種,同一組織: 研究基因在不同處理下,不同條件下的表達變化 ③同一組織,不同物種: 研究基因的進化關系 ④時間序列實驗: 基因在不同時期的表達情況與發(fā)育的關系 2.基因分類: 找到細胞特異,疾病相關,處理相關的基因表達模式,用于診斷疾病和預測等 3.基因網(wǎng)絡和通路: 基因在細胞活動中的功能,基因間的相互作用。 公共數(shù)據(jù)庫 當然,如果你要研究一個基因的功能時,不要先急著去花錢找公司測序, 先去一些基因表達公共數(shù)據(jù)庫找找看: 差異表達(differential expression,DE)基因分析 通過研究基因的差異表達,我們可以發(fā)現(xiàn) ①細胞特異性的基因; ②發(fā)育階段特異性的基因; ③疾病狀態(tài)相關的基因; ④環(huán)境相關的基因; ... 基本方法就是以生物學意義的方式計算基因表達量,然后通過統(tǒng)計學分析表達量尋找具有統(tǒng)計學顯著性差異的基因,從而 ①選擇合適的基因 ②衡量結(jié)果的可靠性 分析方法 尋找差異表達基因有三種方式: E = mean(group1) B=mean(group2) FC = (E-B)/min(E,B) 說人話就是,基因A和基因B的平均值之差與兩者中較小的比值。選擇2-3倍的基因作為結(jié)果 (為什么是2-3倍,就是大家約定俗成)。 fold change 但是簡單粗暴的用2到3倍作為閾值,對于低表達的基因,3倍也是噪音,那些高表達的基因,1.1倍都是生物學顯著了。更重要的沒有考慮到組內(nèi)變異,沒有統(tǒng)計學意義。所以發(fā)文章肯定這個圖只能作為附錄了。 第二種就是統(tǒng)計檢驗,寫文章的時候總需要給出一個p值告訴主編這個結(jié)果可信的(雖然p值也存在爭論)。 注: 基因表達分析的零假設是: 基因在不同處理下的表達量相同。 對于基因芯片的數(shù)據(jù)而言,由于樣本服從正態(tài)分布,所以可以用t-test(雙處理)或anova分析(多處理以上)。 T檢驗適用于只有兩個處理的實驗設計,如植物葉片在相同處理第一天和第二天的基因表達差異。 T-Test 實驗設計 進行T-test檢驗時要注意:是雙尾檢驗(存在差異)還是單尾檢驗(顯著性上調(diào)或下降),兩個樣本的總體是不是等方差(標準T檢驗還是Welch’s test) 如果存在多于兩個處理(條件),就需要用到ANOVA分析了。ANOVA分析能主要是研究結(jié)果之間的差異是如何引起的,具體請移步到我寫方差分析教程。 單因素 雙因素 當然你可以研究,不同基因的表達差異是由因為處理不同,還是基因不同,還是系統(tǒng)誤差,還是其中一些的交互作用。 上面都是針對基因芯片的樣本服從正態(tài)分布進行的統(tǒng)計檢驗?,F(xiàn)在的RNA-Seq,它的抽樣過程是離散的,結(jié)果是count,服從泊松分布,樣本間的差異是服從負二向分布,顯然不能按照上述方法分析。 方差分析(ANOVA)和線性回歸分析(regression)都是同一時期發(fā)展的兩套緊密相連的理論。方差分析考量的是離散型自變量(因子)對連續(xù)型應變量(響應變量)的模型分析,而線性回歸分析只要求響應變量是連續(xù)的,對于自變量無要求。如果響應變量不是連續(xù)型分布,就要使用更加一般化的廣義線性模型(generalized linear model),通過一個連接函數(shù)變換響應變量期望,將響應變量的期望與自變量建立線性關系。 因此,我們可以用廣義線性模型去分析RNA-seq前期分析得到的離散型結(jié)果(count) 方差分析一般用于分析有計劃的實驗結(jié)果,比如說不同處理下的水稻產(chǎn)量?;貧w分析一般分析沒有計劃的數(shù)據(jù),比如說你可以找到大量體檢的數(shù)據(jù),只分析其中性別和身高對體重的影響。所以兩者各有側(cè)重,不要拿大炮轟蚊子。 統(tǒng)計檢驗相對于fold change具有統(tǒng)計意義,不需要參考樣本,需要處理隨機取樣。但是需要重復(ANOVA推薦4-10重復),但由于資金和材料等原因,不一定能夠滿足。此外,對于1000個基因,就要做1000次ANOVA或t-test,最后的p值會有一定的假陽性,因此要做p值矯正(FDR)篩選。 注:推薦在統(tǒng)計檢驗前過濾表達量低,也就如果一個基因在所有樣本中count均低于某一閥值,請在分析前剔除。這個閥值也是約定俗成,一般設置為3. 第三種:Fold Change + 統(tǒng)計檢驗。說一比較尷尬的事情,在統(tǒng)計檢驗中你找到越多的差異表達基因,在p值矯正之后,你反而找不到差異表達基因。也就是說,如果在結(jié)果中存在大量濫竽充數(shù)的所謂的DE基因,那么在嚴格的p值矯正篩選后,反而會誤刪真實的DE基因。 因此在p值矯正之前,你先要手動剔除一部分明顯就是假陽性的DE基因。這個步驟就需要用到前面的fold-change分析。 火山圖 實戰(zhàn)演練 為了方便理解,我們選用同一組數(shù)據(jù)進行實際操作。默認你懂基本的R語言操作,如安裝R包,查看幫助文件等。 · http://www./help/workflows/RnaSeqGeneEdgeRQL/ · http://www./help/workflows/rnaseqGene/ · http://www./help/workflows/RNAseq123/ · http://www./help/workflows/ExpressionNormalizationWorkflow/ 數(shù)據(jù)介紹和下載 用到的數(shù)據(jù)來自于一篇文獻“A pooled shRNA screen for regulators of primary mammary stem and progenitor cells identifies roles for Asap1 and Prox1”。 實驗目的: The aim of this study is to determine the absolute and relative expression levels of mRNA transcripts across multiple flow cytometrically sorted epithelial cell types including freshly isolated CD24+CD29hi mammary stem cell-enriched basal cells (MaSC/basal), CD24+CD29loCD61+ luminal progenitor-enriched (LP) and the CD24+CD29loCD61- mature luminal-enriched (ML) cell populations. Additionally, a comparison between these primary cell types and cultured MaSC/Basal-derived mammosphere cells (mammosphere) and the CommaD-βGeo (CommaDβ) cell line was performed. 實驗設計: Total RNA was extracted and purified from sorted luminal or basal populations from the mammary glands of female virgin 8- to 10-week-old FVB/N mice (3 independent samples for population), MaSC/Basal cells cultured for 1 week under mammsophere conditions and CommaDβ cells grown under maintenance conditions (Deugnier et al. 2006) and their transcriptomes analysed by RNA-Seq. | 根據(jù)文章的介紹,數(shù)據(jù)是100 bp的單端read,用Rsubread比對到mouse reference genome(mm10), 然后使用featureCounts統(tǒng)計每個基因的count數(shù)。然后用TMM進行標準化,轉(zhuǎn)換成log2 counts per million.最后用limma包對每個樣本每個基因的平均表達值以觀察水平權(quán)重的線性模型進行擬合,并用T檢驗找到不同群體的差異表達基因。以FDR + log2-fold-change對基因排序。 我們的任務是下載他們featureCounts的結(jié)果,不用limma,改用DESeq2分析。 # install.packages("R.utils") library("R.utils") url <- "https://www.ncbi.nlm./geo/download/?acc=GSE63310&format=file" utils::download.file(url, destfile="GSE63310_RAW.tar", mode="wb") utils::untar("GSE63310_RAW.tar", exdir = ".") files <- c("GSM1545535_10_6_5_11.txt", "GSM1545536_9_6_5_11.txt", "GSM1545538_purep53.txt","GSM1545539_JMS8-2.txt", "GSM1545540_JMS8-3.txt", "GSM1545541_JMS8-4.txt", "GSM1545542_JMS8-5.txt", "GSM1545544_JMS9-P7c.txt", "GSM1545545_JMS9-P8c.txt") for(i in paste(files, ".gz", sep="")) R.utils::gunzip(i, overwrite=TRUE) 數(shù)據(jù)讀取 先定一個包含所有文件的向量,方便后續(xù)調(diào)用。 files <- c("GSM1545535_10_6_5_11.txt", "GSM1545536_9_6_5_11.txt", "GSM1545538_purep53.txt", "GSM1545539_JMS8-2.txt", "GSM1545540_JMS8-3.txt", "GSM1545541_JMS8-4.txt", "GSM1545542_JMS8-5.txt", "GSM1545544_JMS9-P7c.txt", "GSM1545545_JMS9-P8c.txt") 查看下文件的數(shù)據(jù)格式,分別記錄了每個基因的EntrezID,基因長度和數(shù)量 read.delim(files[1], nrow=5) # EntrezID GeneLength Count #1 497097 3634 1 #2 100503874 3259 0 #3 100038431 1634 0 #4 19888 9747 0 #5 20671 3130 1 之后我們需要用DESeqDataSetFromMatrix為DESeq2提供數(shù)據(jù),命令形式如下: library("DEseq2") DESeqDataSetFromMatrix(countData, colData, design, tidy = FALSE,ignoreRank = FALSE, ...) 其中countData存放counts數(shù)據(jù),colData存放樣本信息的數(shù)據(jù),design就是實驗設計。 首先從各個count matrix文件中讀取count,基因長度部分可以舍棄,因為DESeq2只需要為標準化的count數(shù)據(jù),不需要提供基因長度信息。 sample1 <- read.delim(files[1],header = T)[,c(1,3)] count.table <- data.frame(sample1) for ( f in files[2:length(files)]){ column <- read.delim(f,header=T,row.names = 1)[,2] count.table <- cbind(count.table, column) } head(count.table) rownames(count.table) <- count.table[,1] count.table <- count.table[-1] samplenames <- substring(files, 12, nchar(files)-4) colnames(count.table) <- samplenames 然后要提供colData,其中colData存放列名要和countData的行名相同。 # colData condition <- as.factor(c("LP", "ML", "Basal", "Basal", "ML", "LP", "Basal", "ML", "LP")) lane <- as.factor(rep(c("L004","L006","L008"), c(3,4,2))) colData <- data.frame(condition = condition, lane=lane,row.names = samplenames) all(rownames(colData) %in% colnames(count.table)) 最后導入數(shù)據(jù): dds <- DESeqDataSetFromMatrix(countData = count.table, colData = colData, design = ~ condition ) 數(shù)據(jù)預處理 預過濾 雖然DESeq2會自動屏蔽那些低count的基因,但是剔除那些幾乎不存在基因的部分能夠提高運行速度。 dds <- dds [ rowSums(counts(dds)) > 1, ] rlog和方差齊性轉(zhuǎn)換 許多常見的多維數(shù)據(jù)探索性分析的統(tǒng)計分析方法,例如聚類和主成分分析要求,在那些同方差性的數(shù)據(jù)表現(xiàn)良好。所謂的同方差性就是雖然平均值不同,但是方差相同。 但是對于RNA-Seq count數(shù)據(jù)而言,當均值增加時,方差期望也會提高。也就說直接對count matrix或標準化count(根據(jù)測序深度調(diào)整)做PCA分析,由于高count在不同樣本間的絕對差值大,也就會對結(jié)果有很大影響。簡單粗暴的方法就是對count matrix取log后加1。這個1也是約定俗成,看經(jīng)驗了。 隨便舉個栗子看下效果: lambda <- 10^seq(from = -1, to = 2, length = 1000) cts <- matrix(rpois(1000*100, lambda), ncol = 100) library("vsn") meanSdPlot(cts, ranks = FALSE) 平滑前 log.cts.one <- log2(cts + 1) meanSdPlot(log.cts.one, ranks = FALSE) 平滑后 DESeq2為count數(shù)據(jù)提供了兩類變換方法,使得不同均值的方差趨于穩(wěn)定: regularized-logarithm transformation or rlog(Love, Huber, and Anders 2014)和 variance stabilizing transformation(VST) (Anders and Huber 2010) 用于處理含有色散平均趨勢負二項數(shù)據(jù)。 到底用啥:數(shù)據(jù)集小于30 -> rlog,大數(shù)據(jù)集 -> VST。還有這個處理過程不是用于差異檢驗的,在DESeq分析中會自動選擇最合適的所以你更不需要糾結(jié)了,記得用raw count。 rld <- rlog(dds, blind = FALSE) head(assay(rld), 3) vsd <- vst(dds, blind = FALSE) head(assay(vsd), 3) 處理的結(jié)果如下: library("dplyr") library("ggplot2") dds <- estimateSizeFactors(dds) df <- bind_rows( as_data_frame(log2(counts(dds, normalized=TRUE)[, 1:2]+1)) %>% mutate(transformation = "log2(x + 1)"), as_data_frame(assay(rld)[, 1:2]) %>% mutate(transformation = "rlog"), as_data_frame(assay(vsd)[, 1:2]) %>% mutate(transformation = "vst")) colnames(df)[1:2] <- c("x", "y") ggplot(df, aes(x = x, y = y)) + geom_hex(bins = 80) + coord_fixed() + facet_grid( . ~ transformation) 不同轉(zhuǎn)換的分布 結(jié)果就是轉(zhuǎn)換后更加集中了。 樣本距離 RNA-Seq分析第一步通常是評估樣本間的總體相似度。 ①那些樣本最接近 ②那些樣本差異較大 ③這與實驗設計預期符合么 這里使用R內(nèi)置的 dist 計算 rlog變換數(shù)據(jù)的Euclidean distance,然后用pheatmap根據(jù)結(jié)果畫熱圖 sampleDists <- dist(t(assay(rld))) library("pheatmap") library("RColorBrewer") sampleDistMatrix <- as.matrix( sampleDists ) colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255) pheatmap(sampleDistMatrix, clustering_distance_rows = sampleDists, clustering_distance_cols = sampleDists, col = colors) heatmap 同樣的可以用Poisson Distance (Witten 2011)計算距離,計算方式如下: library("PoiClaClu") poisd <- PoissonDistance(t(counts(dds))) PCA分析 還有一種可視化樣本-樣本距離的方法就是主成分分析。PCA分析我打算找點資料好好理解之后再寫,這個說下有這個方法。 plotPCA(rld,intgroup=c("condition")) PCA 能夠明顯的發(fā)現(xiàn)不同處理的距離離得很遠。 差異表達基因分析(DEA) 在我們定義好DESeqDataSe 就可以方便的調(diào)用DESeq分析. dds <- DESeq(dds) 在幾行輸出后信息后,分析就完成了,更多具體參數(shù)可以用?DESeq查看手冊 構(gòu)建結(jié)果表格 如果直接調(diào)用results,只會提取出design matrix最后兩個變量(這里是ML vs Basal)的 log2 fold changes and p values等結(jié)果。 > results(dds) log2 fold change (MAP): condition ML vs Basal Wald test p-value: condition ML vs Basal DataFrame with 27179 rows and 6 columns baseMean log2FoldChange lfcSE stat pvalue padj <numeric> <numeric> <numeric> <numeric> <numeric> <numeric> 497097 118.7522256 -7.0662875 0.5802473 -12.1780608 4.068401e-34 8.556128e-33 100503874 1.3308906 -2.6371037 1.3370584 -1.9723174 4.857338e-02 8.946685e-02 100038431 0.0843771 -0.1541107 1.2443948 -0.1238439 9.014388e-01 NA 19888 1.4833878 2.9503142 1.3351543 2.2097179 2.712475e-02 5.385925e-02 20671 26.4317752 -1.5256766 0.7754853 -1.9673829 4.913908e-02 9.034136e-02 ... ... ... ... ... ... ... 100861837 368.84237 0.33334299 0.3226665 1.0330883 0.3015626 0.4139846 100861924 22.79272 -1.25169702 0.8307249 -1.5067528 0.1318740 0.2107990 170942 852.61538 -0.07612745 0.2399318 -0.3172879 0.7510252 0.8235139 100861691 0.00000 NA NA NA NA NA 100504472 0.00000 NA NA NA NA NA 實際上results有如下眾多參數(shù) results(object, contrast, name, lfcThreshold = 0, altHypothesis = c("greaterAbs", "lessAbs", "greater", "less"), listValues = c(1, -1), cooksCutoff, independentFiltering = TRUE, alpha = 0.1, filter, theta, pAdjustMethod = "BH", filterFun, format = c("DataFrame", "GRanges", "GRangesList"), test, addMLE = FALSE, tidy = FALSE, parallel = FALSE, BPPARAM = bpparam(), ...) 比如可以指定比較對象Basal和PL,可以用mcols查看結(jié)果存儲的元數(shù)據(jù),了解列名的含義。 res <- results(dds,contrast = c("condition","Basal","LP")) mcols(res, use.names = TRUE) DataFrame with 6 rows and 2 columns type description <character> <character> baseMean intermediate mean of normalized counts for all samples log2FoldChange results log2 fold change (MAP): condition Basal vs LP lfcSE results standard error: condition Basal vs LP stat results Wald statistic: condition Basal vs LP pvalue results Wald test p-value: condition Basal vs LP padj results BH adjusted p-valuess 其中l(wèi)fcSE是Log2 fold change的標準誤。其他項都比較好理解。 summary(res) out of 22026 with nonzero total read count adjusted p-value < 0.1 LFC > 0 (up) : 5528, 25% LFC < 0 (down) : 5274, 24% outliers [1] : 43, 0.2% low counts [2] : 2953, 13% (mean count < 1) [1] see 'cooksCutoff' argument of ?results [2] see 'independentFiltering' argument of ?results 在limma分析結(jié)果分別是4127,4298,稍微多了那么2000個基因。其他結(jié)果也基本上都了2000個。 看起來對于22000多個基因而言,差異不算太大。 當然我們還可以通過 ①降低假陽性率(FDR, false discovery rate)閾值 ②提高log2 fold change的閾值 > res0.5 <- results(dds, contrast = c("condition","Basal","LP"),alpha=0.05) > table(res0.5$padj < 0.05) FALSE TRUE 8860 9750 > resLFC1 <- results(dds, lcontrast = c("condition","Basal","LP"),fcThreshold=1) > table(resLFC1$padj < 0.1) FALSE TRUE 8916 10958 于是乎結(jié)果就和limma差不多了。 多重試驗 在高通量數(shù)據(jù)分析中,我們通常不是用p值來拒絕原假設,更多是用來進行多重試驗矯正。 為什么要做多重試驗矯正? 如果對一個基因,我有99%的把握認為是差異表達的,也就是說1%的可能是錯的。那么假設有10000個基因,按照數(shù)學期望,有100個是假的。因此為了保證多重試驗結(jié)果的可靠性要對結(jié)果的p-value做矯正。 DESeq2用的Benjamini-Hochberg (BH) adjustment (Benjamini and Hochberg 1995),計算矯正后的P-value,用于回答如下問題: 如果我們選擇了所有小于或等于矯正p-value閾值的所有顯著性基因,假陽性比例( false discovery rate, FDR)是多少? FDR在高通量試驗中是比較有用的統(tǒng)計值,由于我們通常關注一類自己感興趣的基因,所以我們需要設置一個假陽性上限。 我們從結(jié)果中以1%作為顯著性,分別找出顯著性上調(diào)和下調(diào)的基因。 regSig <- subset(res, padj < 0.01) # Up head(resSig[ order(resSig$log2FoldChange), ]) # Down head(resSig[ order(resSig$log2FoldChange, decreasing = TRUE), ]) 可視化 人類對圖像比較敏感,對文字比較遲鈍。所以我們需要一些比較好看的圖來放到文章中解釋說明。 最簡單的就是Counts plot,看看特定基因的count數(shù)量。 topGene <- rownames(res)[which.min(res$padj)] plotCounts(dds, gene = topGene, intgroup=c("condition")) count plot 可以用MA-plot(也叫mean-difference plot,Bland-Altman plot)了解模型(如所有基因在不同處理比較的結(jié)果)的估計系數(shù)的分布 res <- lfcShrink(dds, contrast=c("contion","ML","PL"), res=res) plotMA(res, ylim = c(-5, 5)) 更加喜聞樂見的是基因聚類所提供的熱圖展示。我們可以找前20個樣本件差異比較大,然后看他們在不同樣本間的表達情況。 library("genefilter") topVarGenes <- head(order(rowVars(assay(rld)), decreasing = TRUE), 20) mat <- assay(rld)[ topVarGenes, ] mat <- mat - rowMeans(mat) anno <- as.data.frame(colData(rld)[, c("condition","lane")]) pheatmap(mat,annotation_col = anno) heatmap 大致可以發(fā)現(xiàn)同一組的基因顏色是相同的,也就是說表達量相近。 結(jié)語 找到差異表達的基因只是第一步,后續(xù)還需要對這些基因進行進一步的分析,如 ①基因富集分析 ②主成分分析 ③聚類分析 而這些內(nèi)容就是我接下來的學習重點,也是下次更新文章的主題。
回顧 首先對基因表達分析(上)做一個簡單的回顧 1.研究基因表達的有如下工具:RNA-Seq,microarray, qRT-PCR等(歡迎補充) 2.RNA-Seq,microarray一般用在探索性階段,qRT-PCR用于驗證 3.RNA-Seq和microarray由于他們的實驗方式不同,導致尋找差異表達基因的統(tǒng)計學方法也不同。其中microarray使用寡核苷酸作為探針進行雜交,基因表達量與亮度正相關,而亮度是一個連續(xù)型變量,因此大多認為結(jié)果是服從正態(tài)分布。而RNA-Seq的測序結(jié)果是一條條read,是一種離散抽樣過程,因此認為是服從泊松分布。 4.ANOVA和簡單線性模型都是廣義線性模型的特殊情況。ANOVA是研究名義型解釋變量和連續(xù)型解釋變量的關系,簡單線性模式是研究連續(xù)型解釋變量和連續(xù)型解釋變量的關系。而廣義線性模式?jīng)]特殊要求。 5.在3,4的背景下,microarray一般用t檢驗(兩個條件),ANOVA分析(多個條件),最常用limma(線性模型)進行檢驗。RNA-Seq有許多基于count的R包,如DESeq,DESeq2,(基于負二向分布廣義線性模型) 6.以上要求你每個條件都要有3個重復(目前投稿要求),你要是老板窮,一個重復都不給,那你去Google解決方案吧。 7.用R作差異表達分析大致分為以下幾步: 1)根據(jù)軟件包要求導入數(shù)據(jù); 2)數(shù)據(jù)預處理,把那些只有0或1計數(shù)結(jié)果的基因去掉,提高效率。 這一步還可以進行探索性數(shù)據(jù)分析; 3)跑程序,得到結(jié)果; 4)對結(jié)果進行可視化,看看基因聚類等結(jié)果,這一步不是必須的,但卻是展示數(shù)據(jù)最好的手段了。 如果這些內(nèi)容還有印象,那么我們繼續(xù)上次沒有寫完的內(nèi)容的其中一個部分-基因富集分析 為何要基因富集分析 在基因差異表達分析之后,你得到了好多p值特別?。ㄒ簿褪秋@著性很高)的基因,那么下一步你想做什么? ①選擇一些基因用于驗證? ②對其中基因進行后續(xù)研究? ③在結(jié)果中把這些基因都放在后面? ④嘗試著把所有基因相關的文獻都都讀讀看(勸你放棄這個念頭)? 歡迎補充 這些想法都是非常順理成章的,但是不要著急。 什么是基因富集分析 基因富集分析(gene set enrichment analysis)是在一組基因或蛋白中找到一類過表達的基因或蛋白。一般是高通量實驗,如基因芯片,RNA-Seq,蛋白質(zhì)組學(質(zhì)譜結(jié)果)的后續(xù)步驟。 基因富集分析需要我們提供某一類功能基因的集合用于背景,常用的注釋數(shù)據(jù)庫如: ①The Gene Ontology Consortium: 描述基因的層級關系 ②Kyoto Encyclopedia of Genes and Genomes: 提供了pathway的數(shù)據(jù)庫。 分析方法 在文獻Ten Years of Pathway Analysis: Current Approaches and Outstanding Challenges(推薦大家看一遍)作者將研究方法歸為三種: 其中第三種方法想的很好就是難度很大。而且貼心的把每一種方法有哪些工具都總結(jié)出來了: Over-Repressentation Analysis(ORA) ORA是目前商業(yè)化最多的方法。為了說明他的基本思想,我要舉一個喜聞樂見的例子:讀書無用論。 這是我百度找到搜狐財經(jīng)一篇文章《大數(shù)據(jù)告訴你真正的有錢人是什么樣》的有錢人的學歷分布情況,高學歷人群(本科以上,因為本科生太多了)所占的比例是9.4%,其他都是一般學歷占90.6%。這時候,有些公眾號就可以開始不帶腦子的說了,讀書沒什么用呀,有錢人中都是一般學歷的呀,以后讀書讀到大學就行了,甚至也可以不上本科呀(34.2%本科和本科以下)。 實際上,這就是因為沒有考慮到背景。因為高學歷本身人數(shù)就不多,當然在有錢人里面的人數(shù)也就相應不多了。我們要證明有錢人更多是富集高學歷這一部分。
H0: 是否有錢和學歷高無關 richer.pop <- matrix(data = c(10,90,50,950),nrow=2) fisher.test(richer.pop, alternative = "greater") Fisher's Exact Test for Count Data data: richer.pop p-value = 0.03857 alternative hypothesis: true odds ratio is greater than 1 95 percent confidence interval: 1.052584 Inf sample estimates: odds ratio 2.109244 p值小于0.05,看來我讀個博士讓我以后有錢概率變大了。 現(xiàn)在將我們上面的有錢人改成我們找到的基因,整體改成所有基因。高學歷表示屬于目標注釋基因集,一般學歷就是非注釋基因組.我們就是要判斷我們找到的基因更多是在目標注釋集中。 所以你需要列出下表,然后再做一個fisher.test()。
上述的基本思想就是統(tǒng)計學的白球黑球?qū)嶒灒?/span> 在一個黑箱里,有確定數(shù)量的黑白兩種球,你隨機抽?。ú环呕兀㎝個球中,其中兩種球的比例分別是多少? 除了用Fisher精確檢驗,還有其他統(tǒng)計方法: ①Hypergeometric (fisher精確檢驗用的就是超幾何檢驗)http://www./1225.html ②Binomial: 二項分布要求是有放回,無放回要求整體足夠大大到可以近似。 ③Chi-squared chisq.test(counts) ④Z ⑤Kolmogorov-Smirnov ⑥Permutation http://www./1237.html ORA的方法就是如此的簡單,但是有一個問題,就是你如何確定哪些基因是差異表達的,你還是需要設置一個人為的cutoff, 主觀能動性成分有點大。 Functional Class Scoring(FCS) FCS認為,“雖然個體基因表達改變之后會更多在通路中體現(xiàn),但是一些功能相關基因中較弱但協(xié)調(diào)的變化也有明顯的影響?!?/span> The hypothesis of functional class scoring (FCS) is that although large changes in individual genes can have significant effects on pathways, weaker but coordinated changes in sets of functionally related genes (i.e., pathways) can also have significant effects FCS分析方法稍微復雜了一點,他要求的輸入是一個排序的基因列表和一個基因集合。MIT, Broad Institute 2007年文獻就提供了這一方法的軟件"GSEA" the install screen for GSEA 有如下特點: ①計算所有輸入基因集合的分數(shù),而不是單個基因 ②不需要設置cutoff ③找到一組相關的基因 ④提供了更加穩(wěn)健的統(tǒng)計框架 GSEA是一款圖形化的軟件,根據(jù)他們提供的教程,然后點呀點,就會得到如下結(jié)果。 下圖就是需要好好理解的部分。 GESA 中間從藍色到紅色的過渡“帶”表示基因從上調(diào)到下調(diào)排列(排序可以按照fold change,也可以是p-value)。黑色像條形碼的豎線表示該位置的基因?qū)儆谀硞€指定通路。綠色有波動的曲線表示富集分數(shù),從0開始計算,屬于基因通路增加,不屬于則減少。最后看下黑色的條形碼是不是富集在一端。 那如何做統(tǒng)計檢驗呢? The final step in FCS is assessing the statistical significance of the pathway-level statistic. When computing statistical significance, the null hypothesis tested by current pathway analysis approaches can be broadly divided into two categories: i) competitive null hypothesis and ii) self-contained null hypothesis [3], [18], [22], [31]. A self-contained null hypothesis permutes class labels (i.e., phenotypes) for each sample and compares the set of genes in a given pathway with itself, while ignoring the genes that are not in the pathway. On the other hand, a competitive null hypothesis permutes gene labels for each pathway, and compares the set of genes in the pathway with a set of genes that are not in the pathway。 我們要檢驗的目標是基因富集在一端是因為于目標通路相關的基因都在一端富集。那么空假設就是,你把找到的基因隨便擺放也能看到富集現(xiàn)象。用比較專業(yè)的話說就是先生成一個零假設的數(shù)據(jù)分布,然后觀察實際數(shù)據(jù)在這個零假設分布下,是不是在尾端。 一些問題 統(tǒng)計檢驗的能力是有限的,所以還有很多問題存在解決。 ①我們希望找到生物學顯著的基因,但是生物學顯著和統(tǒng)計顯著兩者并不是完全相關 ②無論是ORA還是FCS都對背景(也就是這個物種一共有多少基因)有要求,但是隨著我們的研究深入,基因數(shù)量會改變。 有些軟件會直接設置一個很大的背景數(shù),從而讓p值很顯著,然后我們就開心地用他們的結(jié)果。 ③有些基因沒有注釋,也就是注釋缺失,處理方法就是扔(歡迎拍磚)。 ④有一些注釋項是其他項的子集。 實戰(zhàn):用clusterProfiler做富集分析 clusterProfiler是Y叔良心之作(當然他的良心之作還有很多),目前支持KEGG在線拉數(shù)據(jù),支持DAVID,支持Broad iNSTITUTE Molecular Signatures Databases, 支持GSEA。你問我支不支持,我當然是支持的啦。所以這里放上Y叔的公眾號了。 安裝: source("https:///biocLite.R") biocLite('clusterProfiler') 這里簡單舉例如何使用,更多內(nèi)容見https://guangchuangyu./clusterProfiler/, 這個部分你可以直接過掉,因為我只是跟著Y叔的代碼敲了一遍而已。 加載差異表達基因,我這里偷懶就隨機挑一些基因名(來自于DOSE包,Disease Ontology Semantic and Enrichment analysis )出來了。 library(org.Hs.eg.db) data(geneList) gene <- names(geneList)[abs(geneList) > 2] head(gene) [1] "4312" "8318" "10874" "55143" "55388" "991" Y叔為了展示他能夠處理不同命名方式的ID,用bitr(來自于clusterProfiler)進行生物學ID轉(zhuǎn)換 gene.df <- bitr(gene, fromType = "ENTREZID", toType = c("ENSEMBL", "SYMBOL"), OrgDb = org.Hs.eg.db) head(gene.df) ENTREZID ENSEMBL SYMBOL 1 4312 ENSG00000196611 MMP1 2 8318 ENSG00000093009 CDC45 3 10874 ENSG00000109255 NMU 4 55143 ENSG00000134690 CDCA8 5 55388 ENSG00000065328 MCM10 6 991 ENSG00000117399 CDC20 clusterProfiler: ORA 然后對不同命名的ID都做ORA富集分析。 ego <- enrichGO(gene = gene, universe = names(geneList), OrgDb = org.Hs.eg.db, ont = "CC", pAdjustMethod = "BH", pvalueCutoff = 0.01, qvalueCutoff = 0.05) ego2 <- enrichGO(gene = gene.df$ENSEMBL, OrgDb = org.Hs.eg.db, keytype = 'ENSEMBL', ont = "CC", pAdjustMethod = "BH", pvalueCutoff = 0.01, qvalueCutoff = 0.05) ego3 <- enrichGO(gene = gene.df$SYMBOL, OrgDb = org.Hs.eg.db, keytype = 'SYMBOL', ont = "CC", pAdjustMethod = "BH", pvalueCutoff = 0.01, qvalueCutoff = 0.05) 結(jié)果會有如下 所以從GeneRatio:24/199,BgRatio:231/11711可以反推出下表:
最后再畫一張美美的點圖,根據(jù)GeneRatio所做。 dotplot(ego, showCategory=30) clusterProfiler: GSEC clusterProfiler支持GSEC,而且很好用。 gsecc <- gseGO(geneList=geneList, ont="CC", OrgDb=org.Hs.eg.db, verbose=F) head(summary(gsecc) gseaplot(gsecc, geneSetID="GO:0000779" |
|