一区二区三区日韩精品-日韩经典一区二区三区-五月激情综合丁香婷婷-欧美精品中文字幕专区

分享

基因表達分析(上) 基因表達

 wangprince2017 2018-07-17

什么是基因表達,如下是來自于維基百科的解釋:

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.

http://upload-images./upload_images/2013053-1af3f38efbc4ff83.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/580https://upload-images./upload_images/2013053-1af3f38efbc4ff83.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/580

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ù)雜交信號的強弱來判斷基因表達的程序。

https://upload-images./upload_images/2013053-e7b6a62c3e0d3cae.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/700

microarray

但是microarray檢測的基因數(shù)量完全取決于你的探針設計的數(shù)量,而且難以研究mRNA的可變剪切。

RNA-Seq

RNA-Seq是目前基因表達分析最常用的技術。分為以下幾步

1.分離所有mRNA

2.逆轉(zhuǎn)錄mRNA成cDNA

3.對cDNA測序

4.比對參考基因組

RNA-Seq實驗設計中的“重復”包括:技術重復生物學重復
重復是為了檢測組間和組內(nèi)的變異,對于假設檢驗至關重要。

(1)技術重復為了估計測量技術(RNA-Seq)的變異。

(2)生物學重復是為了發(fā)現(xiàn)生物組內(nèi)的變異。
簡單的說,兩組的基因表達的變化只有比組內(nèi)變異還大時才能認為時顯著的。

RNA-Seq的概率分布
相同基因在不同細胞的表達水平服從log-normal(對數(shù)正態(tài))分布,由定量PCR驗證。

(注:這與相同細胞不同基因表達的分布不同)

但是大多數(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ù)庫找找看:

1.http://www./arrayexpress/

2.https://www.ncbi.nlm./geo/

差異表達(differential expression,DE)基因分析

通過研究基因的差異表達,我們可以發(fā)現(xiàn)

①細胞特異性的基因;

②發(fā)育階段特異性的基因;

疾病狀態(tài)相關的基因;

④環(huán)境相關的基因;

...

基本方法就是以生物學意義的方式計算基因表達量,然后通過統(tǒng)計學分析表達量尋找具有統(tǒng)計學顯著性差異的基因,從而

①選擇合適的基因

②衡量結(jié)果的可靠性

分析方法

尋找差異表達基因有三種方式:
第一種是計算Fold change(倍數(shù)變化),十分簡單粗暴的方法,計算方法如下:

E = mean(group1)

B=mean(group2)

FC = (E-B)/min(E,B)

說人話就是,基因A和基因B的平均值之差與兩者中較小的比值。選擇2-3倍的基因作為結(jié)果

(為什么是2-3倍,就是大家約定俗成)。

http://upload-images./upload_images/2013053-5e7e56bf5084f33c.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/673https://upload-images./upload_images/2013053-5e7e56bf5084f33c.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/673

fold change

但是簡單粗暴的用2到3倍作為閾值,對于低表達的基因,3倍也是噪音,那些高表達的基因,1.1倍都是生物學顯著了。更重要的沒有考慮到組內(nèi)變異,沒有統(tǒng)計學意義。所以發(fā)文章肯定這個圖只能作為附錄了。

第二種就是統(tǒng)計檢驗,寫文章的時候總需要給出一個p值告訴主編這個結(jié)果可信的(雖然p值也存在爭論)。
復習一下:
p值指的碰巧是拒絕零假設機會。P值越大假陽性越低,同時真實結(jié)果也可能會剔除。

注: 基因表達分析的零假設是: 基因在不同處理下的表達量相同。

對于基因芯片的數(shù)據(jù)而言,由于樣本服從正態(tài)分布,所以可以用t-test(雙處理)或anova分析(多處理以上)。

T檢驗適用于只有兩個處理的實驗設計,如植物葉片在相同處理第一天和第二天的基因表達差異。

https://upload-images./upload_images/2013053-778332ebdd08d66f.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/473

T-Test 實驗設計

進行T-test檢驗時要注意:是雙尾檢驗(存在差異)還是單尾檢驗(顯著性上調(diào)或下降),兩個樣本的總體是不是等方差(標準T檢驗還是Welch’s test)

如果存在多于兩個處理(條件),就需要用到ANOVA分析了。ANOVA分析能主要是研究結(jié)果之間的差異是如何引起的,具體請移步到我寫方差分析教程。
   對于基因表達而言,研究目標是,對于同一個基因而言,他們之間的差異是處理不同造成,還是因為系統(tǒng)誤差造成。

https://upload-images./upload_images/2013053-77f16e3c9dafcf56.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/700

單因素

https://upload-images./upload_images/2013053-b60f84f70e63b6a7.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/700

雙因素

當然你可以研究,不同基因的表達差異是由因為處理不同,還是基因不同,還是系統(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分析。
   我們可以通過火山圖來看看如何確定區(qū)間

https://upload-images./upload_images/2013053-2afae264bafebcec.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/700

火山圖

實戰(zhàn)演練

為了方便理解,我們選用同一組數(shù)據(jù)進行實際操作。默認你懂基本的R語言操作,如安裝R包,查看幫助文件等。
正式學習之前,先感謝一下
Bioconductor。 根據(jù)他們提供流程,我學習到如何進行RNA-Seq分析。

·        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ù),不需要提供基因長度信息。
    邏輯就是分別讀取每一個文件的count列,然后賦予文件名。

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)

https://upload-images./upload_images/2013053-29a942b7df79024d.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/360

平滑前

log.cts.one <- log2(cts + 1)

meanSdPlot(log.cts.one, ranks = FALSE)

https://upload-images./upload_images/2013053-b61beaa7d2fc8f01.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/360

平滑后

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)

https://upload-images./upload_images/2013053-2872940fbd05c2de.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/432

不同轉(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)

https://upload-images./upload_images/2013053-e105e80e704afd6e.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/581

heatmap

同樣的可以用Poisson Distance (Witten 2011)計算距離,計算方式如下:

library("PoiClaClu")

poisd <- PoissonDistance(t(counts(dds)))

PCA分析

還有一種可視化樣本-樣本距離的方法就是主成分分析。PCA分析我打算找點資料好好理解之后再寫,這個說下有這個方法。
DESeq2提供了專門的方法用于作圖,還很好看呢!

plotPCA(rld,intgroup=c("condition"))

https://upload-images./upload_images/2013053-0cdc8f67a4eeb435.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/581

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的標準誤。其他項都比較好理解。
最后還可以看一些總結(jié)性的內(nèi)容

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"))

https://upload-images./upload_images/2013053-4e01b6e17d7a6965.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/581

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)

https://upload-images./upload_images/2013053-f18267c2a1abcff9.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/581

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é)果中把這些基因都放在后面?

④嘗試著把所有基因相關的文獻都都讀讀看(勸你放棄這個念頭)?

歡迎補充

這些想法都是非常順理成章的,但是不要著急。
首先,差異表達找到的基因往往很多,你簡單的粗暴去找每一個基因的詳細資料,顯然不太現(xiàn)實;
其次,如果我們單純覺得某一個基因和你研究的課題相關,或者說你其實已經(jīng)找到了一個有可能的基因(或者你只是希望用一些高大上的實驗驗證一下)那么這個行為是不是有太多主觀性,存在一些偏見。
當然,你覺得基因就是你要找的,可是萬一它只是碰巧來打醬油的呢,這不就是很尷尬了。
所以為了讓審稿人相信你的結(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(推薦大家看一遍)作者將研究方法歸為三種:

https://upload-images./upload_images/2013053-9e6c72d024bc0bee?imageMogr2/auto-orient/strip%7CimageView2/2/w/700

其中第三種方法想的很好就是難度很大。而且貼心的把每一種方法有哪些工具都總結(jié)出來了:

https://upload-images./upload_images/2013053-9da6e8cd419c223f?imageMogr2/auto-orient/strip%7CimageView2/2/w/700

Over-Repressentation Analysis(ORA)

ORA是目前商業(yè)化最多的方法。為了說明他的基本思想,我要舉一個喜聞樂見的例子:讀書無用論。

            https://upload-images./upload_images/2013053-3594b9b9d36ce381.jpeg?imageMogr2/auto-orient/strip%7CimageView2/2/w/436

這是我百度找到搜狐財經(jīng)一篇文章《大數(shù)據(jù)告訴你真正的有錢人是什么樣》的有錢人的學歷分布情況,高學歷人群(本科以上,因為本科生太多了)所占的比例是9.4%,其他都是一般學歷占90.6%。這時候,有些公眾號就可以開始不帶腦子的說了,讀書沒什么用呀,有錢人中都是一般學歷的呀,以后讀書讀到大學就行了,甚至也可以不上本科呀(34.2%本科和本科以下)。
    你每年回家總能回去看到有人炫耀說,雖然我有錢,可是讀書太少了,都不能和你們讀書人比的。你總感覺哪里不對勁,但是卻又不太方便說出來。

實際上,這就是因為沒有考慮到背景。因為高學歷本身人數(shù)就不多,當然在有錢人里面的人數(shù)也就相應不多了。我們要證明有錢人更多是富集高學歷這一部分。

類別

有錢人

整體

高學歷

10

50

一般學歷

90

950


100

1000

H0: 是否有錢和學歷高無關
Ha: 學歷高還是有點用的
然后做一個Fisher精確檢驗,看看p值。

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()。

類別

gene list

Genome

in anno group

10

50

not in anno group

290

19950


300

20000

上述的基本思想就是統(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"

https://upload-images./upload_images/2013053-aba7e105f72eb93f.gif?imageMogr2/auto-orient/strip%7CimageView2/2/w/440

the install screen for GSEA

有如下特點:

①計算所有輸入基因集合的分數(shù),而不是單個基因

②不需要設置cutoff

③找到一組相關的基因

提供了更加穩(wěn)健的統(tǒng)計框架

GSEA是一款圖形化的軟件,根據(jù)他們提供的教程,然后點呀點,就會得到如下結(jié)果。

下圖就是需要好好理解的部分。

https://upload-images./upload_images/2013053-f5aa582b2677bbcd.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/700

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é)果會有如下
ID:          基因本體論的ID
Description: 描述
GeneRatio:   在GO詞條所占的比例
BgRatio :    在背景所占的比例
pvalue:      假設是正確但是被拒絕的概率
p.adjust:    采用BH方法進行多重試驗p值校正
qvalue:      Q值=被拒絕但卻是正確的概率
geneID:     列出在這個GO詞條下的我們提供的基因
count:      數(shù)量,如果沒有約分,就是GeneRatio的分子了。

所以從GeneRatio:24/199,BgRatio:231/11711可以反推出下表:

類別

gene list

Genome

in anno group

24

199

not in anno group

231

11711


255

11910

最后再畫一張美美的點圖,根據(jù)GeneRatio所做。

dotplot(ego, showCategory=30)

https://upload-images./upload_images/2013053-dd7aa9d6db37b272.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/700

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"

https://upload-images./upload_images/2013053-9eb4719e76c7501c.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/700

    本站是提供個人知識管理的網(wǎng)絡存儲空間,所有內(nèi)容均由用戶發(fā)布,不代表本站觀點。請注意甄別內(nèi)容中的聯(lián)系方式、誘導購買等信息,謹防詐騙。如發(fā)現(xiàn)有害或侵權(quán)內(nèi)容,請點擊一鍵舉報。
    轉(zhuǎn)藏 分享 獻花(0

    0條評論

    發(fā)表

    請遵守用戶 評論公約

    類似文章 更多

    色一欲一性一乱—区二区三区| 男女激情视频在线免费观看| 高清不卡视频在线观看| av在线免费观看一区二区三区| 亚洲中文字幕日韩在线| 亚洲妇女黄色三级视频| 国产欧美日韩精品自拍| 午夜精品黄片在线播放| 欧美成人久久久免费播放| 黄片在线免费观看全集| 亚洲少妇人妻一区二区| 亚洲熟女乱色一区二区三区| 99久久国产亚洲综合精品| 国产精品偷拍一区二区| 日韩中文字幕在线不卡一区| 日本99精品在线观看| 国产成人精品资源在线观看| 欧美日韩国产欧美日韩| 都市激情小说在线一区二区三区| 五月激情五月天综合网| 国产午夜福利在线免费观看| 精品少妇人妻av一区二区蜜桃| 国产原创激情一区二区三区| 久久久精品区二区三区| 亚洲中文字幕有码在线观看| 男人和女人黄 色大片| 精品一区二区三区乱码中文| 麻豆视频传媒入口在线看| 中国美女草逼一级黄片视频| 国产色一区二区三区精品视频| 综合久综合久综合久久| 91亚洲人人在字幕国产| 五月情婷婷综合激情综合狠狠| 日韩国产精品激情一区| 亚洲专区中文字幕视频| 亚洲精品国产第一区二区多人| 99一级特黄色性生活片| 午夜激情视频一区二区| 国产色第一区不卡高清| 九九热视频免费在线视频| 东北老熟妇全程露脸被内射|