RNA Dizilimi için Üç Diferansiyel İfade Analizi Yöntemi: limma, EdgeR, DESeq2

Özet

RNA dizilimi için diferansiyel ekspresyon analiz yöntemlerinin ayrıntılı bir protokolü sağlanmıştır: limma, EdgeR, DESeq2.

Özet

RNA dizilimi (RNA-seq), genetik değişim ve karmaşık biyolojik süreçler arasındaki ilişkiyi ortaya çıkarabildiği ve tümörlerin tanı, prognostik ve terapötiklerinde büyük değere sahip olduğu için transkriptomikte en yaygın kullanılan teknolojilerden biridir. RNA-seq verilerinin diferansiyel analizi anormal transkripsiyonları tanımlamak için çok önemlidir ve limma, EdgeR ve DESeq2 diferansiyel analiz için etkili araçlardır. Bununla birlikte, RNA-seq diferansiyel analizi, tıp eğitimi müfredatında eksik olan R dili ve uygun bir yöntem seçme yeteneği ile belirli beceriler gerektirir.

Burada, sırasıyla limma, DESeq2 ve EdgeR aracılığıyla kolanjiokarsinom (CHOL) ve normal dokular arasında farklı olarak ifade edilen genleri (DEG' ler) tanımlamak için ayrıntılı protokol sunuyoruz ve sonuçlar volkan arazilerinde ve Venn diyagramlarında gösteriliyor. Limma, DESeq2 ve EdgeR'ın üç protokolü benzerdir, ancak analiz süreçleri arasında farklı adımlara sahiptir. Örneğin, doğrusal bir model limma istatistikleri için kullanılırken, negatif binom dağılımı edgeR ve DESeq2'de kullanılır. Ayrıca, normalleştirilmiş RNA-seq sayısı verileri EdgeR ve limma için gereklidir, ancak DESeq2 için gerekli değildir.

Burada, üç diferansiyel analiz yöntemi için ayrıntılı bir protokol sunuyoruz: limma, EdgeR ve DESeq2. Üç yöntemin sonuçları kısmen çakışıyor. Her üç yöntemin de kendi avantajları vardır ve yöntem seçimi yalnızca verilere bağlıdır.

Giriş

RNA dizilimi (RNA-seq), transkriptomikte birçok avantajla (örneğin, yüksek veri tekrarlanabilirliği) en yaygın kullanılan teknolojilerden biridir ve karmaşık biyolojik süreçlerin işlevleri ve dinamikleri hakkında anlayışımızı önemli ölçüde artırmıştır^1,2. Farklı olarak ifade edilen genler (DEG' ler) olarak da bilinen farklı biyolojik bağlam altında sapma transkriptlerinin tanımlanması, RNA-seq analizinde önemli bir adımdır. RNA-seq, patogenezle ilgili moleküler mekanizmaların ve biyolojik fonksiyonların derinlemesine anlaşılmasını mümkün kılar. Bu nedenle, ayırıcı analiz tümörlerin tanı, prognostik ve terapötikleri için değerli olarak kabul edilmiştir^3,4,5. Şu anda, RNA-seq diferansiyel ekspresyon analizi için daha açık kaynaklı R / Biyoiletken paketleri geliştirilmiştir, özellikle limma, DESeq2 ve EdgeR^1,6,7. Bununla birlikte, diferansiyel analizi, R dili ile belirli beceriler ve tıp eğitimi müfredatında eksik olan uygun yöntemi seçme yeteneğini gerektirir.

Bu protokolde, Kanser Genom Atlası'ndan (TCGA) çıkarılan kolanjiokarsinom (CHOL) RNA-seq sayım verilerine dayanarak, CHOL ve normal dokular arasındaki DEG'leri tanımlamak için R programı¹¹tarafından en bilinen yöntemlerden üçü (sırasıyla limma⁸, EdgeR⁹ ve DESeq2¹⁰⁾ gerçekleştirilmiştir. Limma, EdgeR ve DESeq2'nin üç protokolü benzerdir, ancak analiz süreçleri arasında farklı adımlara sahiptir. Örneğin, normalleştirilmiş RNA-seq sayısı verileri EdgeR ve limma^8,9için gereklidir, DESeq2 ise normalleştirme¹⁰yerine verileri düzeltmek için kendi kitaplık tutarsızlıklarını kullanır. Ayrıca, edgeR RNA-seq verileri için özel olarak uygundur, limma ise mikroarraylar ve RNA-seq için kullanılır. Doğrusal bir model limma tarafından DEG'leri değerlendirmek için benimsenmiştir¹², edgeR'daki istatistikler ampirik Bayes tahmini, kesin testler, genelleştirilmiş doğrusal modeller ve yarı olasılık testleri⁹dahil olmak üzere negatif binom dağılımlarına dayanmaktadır.

Özetle, sırasıyla limma, DESeq2 ve EdgeR kullanarak RNA-seq diferansiyel ekspresyon analizinin ayrıntılı protokollerini sağlıyoruz. Bu makaleye atıfta bulunarak, kullanıcılar RNA-seq diferansiyel analizini kolayca gerçekleştirebilir ve verileri için uygun diferansiyel analiz yöntemlerini seçebilirler.

Protokol

NOT: R-studio programını açın ve R dosyası "DEGs.R" yükleyin, dosya Ek dosyalardan / Komut Dosyalarından edinilebilir.

1. Verilerin indirilmesi ve önceden işlenmesi

1. Kanser Genom Atlası'ndan (TCGA) kolanjiokarsinom (CHOL) yüksek verimli dizileme (HTSeq) sayısı verilerini indirin. Bu adım aşağıdaki R kodu ile kolayca elde edilebilir.
   1. R paketlerini yüklemek için Çalıştır'ı tıklatın.
   2. R paketlerini yüklemek için Çalıştır'ı tıklatın.
      if(!requireNamespace("BiocManager", sessizce=DOĞRU))
      + install.packages("BiocManager")
      BiocManager::install(c("TCGAbiolinks", "SummarizedExperiment"))
   3. Çalışma dizinini ayarlayın.
      kütüphane (TCGAbiolinks)
      library(SummarizedExperiment)
      setwd("C:/Kullanıcılar/LIUSHIYI/Masaüstü")
   4. Kanser türünü seçin.
      kanser <- "TCGA-CHOL"
   5. Verileri indirmek için "GDCquery.R" dosyasındaki R kodunu çalıştırın. "GDCquery.R" dosyası Tamamlayıcı dosyalardan/Komut Dosyalarından edinilebilir:
      kaynak("Tamamlayıcı dosyalar/Komut Dosyaları/GDCquery.R")
      kafa(cnt)
      ##TCGA-3X-AAVA-01A-11R-A41I-07
      ##ENSG00000000003 4262
      ##ENSG00000000005 1
      ##ENSG00000000419 1254
      ##ENSG00000000457 699
      ##ENSG00000000460 239
      ##ENSG00000000938 334
      NOT: Yürütmeden sonra, CHOLHTSeq sayısı verileri indirilir ve satırların topluluk gen kimliklerini ve sütunların örnek kimlikleri temsil ettiği "cnt" olarak adlandırılır. Örnek kimliklerdeki 14-15 numaralarına dikkat edin; 01 ile 09 arasında değişen sayılar tümörleri, 10 ile 19 arasında değişen sayılar ise normal dokuları gösterir.
2. Topluluk gen kimliklerini gen sembollerine dönüştürün.
   1. Ek açıklama dosyasını depolama yoluna göre R'ye alın. Ek açıklama dosyası (gencode.v22.annotation.gtf) Tamamlayıcı dosyalardan edinilebilir.
      gtf_v22 <- rtracklayer::import('Tamamlayıcı dosyalar/gencode.v22.annotation.gtf')
   2. R kodunu "gtf_v22 çalıştır. Ek dosyalardan/Komut Dosyalarından alınabilen R" dosyası:
      source("Tamamlayıcı dosyalar/Komut Dosyaları/gtf_v22. R")
   3. Topluluk gen kimliklerini gen sembollerine dönüştürmek için "ann" işlevini uygulayın.
      cnt=ann(cnt,gtf_v22)
3. Düşük ekspresyonlu genleri filtreleme
   1. R paketi "edgeR"ı yüklemek için Çalıştır'ı tıklatın.
      BiocManager::install("edgeR")
   2. R paketi "edgeR"ı yüklemek için Çalıştır'ı tıklatın.
      kitaplık(edgeR)
   3. En az iki örnekten birinden büyük milyon (BGBM) değerine sahip genleri tutmak için aşağıdaki R kodunu çalıştırın.
      < satır tutSums(cpm(cnt)>1)>=2
      cnt <- as.matrix(cnt[keep,])
      NOT: Farklı sıralama derinliklerinin neden olduğu sapmayı ortadan kaldırmak için okuma sayısı yerine milyon başına sayım (BGBM) değeri kullanılır.

2. "Limma" ile diferansiyel ekspresyon analizi

1. R paketi "limma"yı yüklemek için Çalıştır'ı tıklatın.
   BiocManager::yükle("limma")
2. R paketlerini "limma", "edgeR" yüklemek için Çalıştır'ı tıklatın.
   kütüphane(limma)
   kitaplık(edgeR)
3. Tasarım matrisini oluşturmak için aşağıdaki R kodunu çalıştırın.
   grup <- substring(colnames(cnt),14,15) # Extract group information
   grup [grup %in% "01"] <- "Cancer" # set '01' as tumor tissue
   grup [grup %in% "11"] <- "Normal" # set '11' as normal tissue
   grup <- factor (group, levels = c("Normal","Cancer"))
   1. Tasarım matrisini oluşturun.
      tasarım <- model.matrix (~grup)
      satır adları(tasarım) <- sütun adları(cnt)
   2. DGEList nesnesini oluşturun.
      dge <- DGEList(sayı = cnt, grup = grup)
   3. Verileri normalleştirin.
      dge <- calcNormFactors(dge, yöntem = "TMM")
   4. Limma eğilim yöntemi tabanlı diferansiyel ifade çözümlemesi gerçekleştirmek için aşağıdaki R kodunu çalıştırın.
      dge
      ##An "DGEList" sınıfının nesnesi
      ##$counts
      ##TCGA-3X-AAVA-01A-11R-A41I-07
      ##TSPAN6 4262
      ##DPM1 1254
      ##SCYL3 699
      ##C1orf112 239
      ##FGR 334
   5. BGBM değerini hesaplayın.
      logdge <- cpm(dge, log=TRUE, prior.count=3)
   6. Verileri tahmin etmek veya değişkenler arasındaki ilişkiyi çıkarmak için doğrusal bir model sığdırmak için Çalıştır'ı tıklatın.
      fit <- lmFit (logdge, tasarım)
   7. Bayesian'a göre T değerini, F değerini ve günlük oranlarını hesaplayın.
      fit <- eBayes (fit, trend=DOĞRU)
   8. Sonuç tablosunu ayıklayın.
      res_limma<- as.data.frame(topTable(sığdırma,n=Inf))
      
      kafa(res_limma)
      ## logFC AveExpr t P.Value adj. P.Val B
      ##RP11-252E2.2 -4.899493 -2.488589 -20.88052 2.386656e-25 4.931786e-21 47.28823
      ##BX842568,1 -4,347930 -2,595205 -20,14532 1,082759e-24 1,118706e-20 45,83656
      ##CTC-537E7.3 -5.154894 -2.143292 -19.59571 3.452354e-24 2.216114e-20 44.72001
      ##RP11-468N14.3 -6.532259 -2.029714 -19.49409 4.289807e-24 2.216114e-20 44.51056
      ##AP006216,5 -4,507051 -2,670915 -19,25649 7,153356e-24 2,956339e-20 44,01704
      ##RP11-669E14.4 -4.107204 -2.828311 -18.93246 1.448209e-23 4.987633e-20 43.33543
      #The ifade analizinin sonucu, gen kimliği, log2 kat değişim değeri (logFC), deneydeki genin ortalama log2 ifade düzeyi (AveExpr), değiştirilmiş t istatistiği (t), relavent p değeri (P.Value), yanlış keşif hızı (FDR) düzeltilmiş p değerini (adj) içeren "res_limma" kaydedilir. P.Val) ve farklı ifade edilen genlerin günlük oranları (B)
      NOT: "edgeR"ın "calcNormFactors()" işlevi, örnek hazırlama veya kitaplık yapımı ve diziliminin neden olduğu etkiyi ortadan kaldırmak için verileri normalleştirmek için kullanılmıştır. Tasarım matrisinin yapımında, deneysel tasarımın (örneğin, doku tipi: normal veya tümör dokuları) matrisin örnek kimlikleriyle eşleşmesi gerekir. limma-trend, sıralama derinliği aynı olan veriler için uygundur, limma-voom ise uygundur: (i) örnek kitaplık boyutu farklı olduğunda; (ii) TMM tarafından normalleştirilmeyen veriler; (iii) verilerde çok fazla "gürültü" vardır. Pozitif logFC, genin deneyde yukarı doğru düzenlendiği anlamına gelirken, negatif sayı genin aşağı regüle edildiği anlamına gelir.
   9. DEG'leri tanımlayın.
      res_limma$sig <- as.factor(
      ifelse(res_limma$adj. P.Val < 0.05 & abs(res_limma$logFC) > 2,
      ifelse(res_limma$logFC > 2 ,'up','down'),'not')) # 0,05 < adj.p Değeri ve |log2FC| >= 2, DEG'leri tanımlamak için eşiklerdir
      özet(res_limma$sig)
      ##down yukarı değil
      ##1880 ​17341 1443
   10. Sonuç tablosunu bir dosyaya çıktısını verin.
       yazma.csv(res_limma, dosya = 'result_limma.csv')
   11. R paketi "ggplot2"yi yüklemek için Çalıştır'ı tıklatın.
       install.packages("ggplot2")
   12. R paketi "ggplot2"yi yüklemek için Çalıştır'ı tıklatın.
       kitaplık(ggplot2)
   13. "Volkan"ın R kodunu çalıştır. R" volkan arsa oluşturmak için. "Volkan" dosyası. R" Ek dosyalardan edinilebilir.
       kaynak("Tamamlayıcı dosyalar/Komut Dosyaları/volkan. R")
       volkan(res_limma,"logFC","adj. P.Val",2,0.05)
       NOT: Genler log2FC ve adj-p değerlerine göre farklı pozisyonlara eşlenebilir, yukarı düzenlenmiş DEG'ler kırmızı renklendirilir ve aşağı düzenlenmiş DEG'ler yeşil renklendirilir.
   14. Volkan arsasını kaydetmek için Dışarı Aktar'ı tıklatın.
       NOT: Volkan çizimleri farklı formatlarda oluşturulabilir ve indirilebilir (örneğin, pdf, TIFF, PNG, JPEG formatı). Genler log2FC ve adj p değerlerine göre farklı pozisyonlara eşlenebilir, yukarı düzenlenmiş DEG'ler (log2FC > 2, adj p < 0.05) kırmızı renktedir ve aşağı düzenlenmiş DEG'ler (log2FC < -2, adj p < 0.05) yeşil renkte, DEG olmayanlar gri renkte renklendirilir.

3. "edgeR" ile diferansiyel ifade analizi

1. R paketi "edgeR"ı yüklemek için Çalıştır'ı tıklatın.
   kitaplık(edgeR)
2. Tasarım matrisi oluşturmak için aşağıdaki R kodunu çalıştırın.
   grup <-alt dize(sütun adları(cnt),14,15)
   grup [grup %in% "01"] <- "Kanser"
   grup [grup %in% "11"] <- "Normal"
   group=faktör(grup, düzey = c("Normal","Kanser"))
   tasarım <-model.matrix(~grup)
   satır adları(tasarım) = sütun adları(cnt)
3. DGEList nesnesini oluşturmak için Çalıştır'ı tıklatın.
   dge <- DGEList(counts=cnt)
4. Verileri normalleştirin.
   dge <- calcNormFactors(dge, yöntem = "TMM")
5. Gen ifade değerlerinin dağılımını tahmin etmek için Çalıştır'ı tıklatın.
   dge <- estimateDisp(dge, tasarım, sağlam = T)
6. Verileri saymak için modeli sığdırmak için Çalıştır'ı tıklatın.
   fit <- glmQLFit(dge, tasarım)
7. İstatistiksel bir test yap.
   fit <- glmQLFTest(fit)
8. Sonuç tablosunu ayıklayın. Sonuç, günlük katlama değişim değeri, günlük BGBM, F, p değeri ve FDR düzeltilmiş p değerini içeren "res_edgeR" kaydedilir.
   res_edgeR=as.data.frame(topTags(sığdırma, n=Inf))
   kafa(res_edgeR)
   ## logFC logCPM F PValue FDR
   ##GCDH -3,299633 5,802700 458.5991 1,441773e-25 2,979280e-21
   ##MSMO1 -3,761400 7,521111 407.0416 1,730539e-24 1,787993e-20R
   ##CL1 -3,829504 5,319641 376.5043 8,652474e-24 5,516791e-20
   ##ADI1 -3,533664 8,211281 372.6671 1,067904e-23 5,516791e-20
   ##KCNN2 -5,583794 3,504017 358.6525 2,342106e-23 9,679455e-20
   ##GLUD1 -3,287447 8,738080 350.0344 3,848408e-23 1,194406e-19
   #The sonuç, günlük katlama değişim değeri (logFC), günlük BGBM, F, p değeri ve FDR düzeltilmiş p değerini içeren "res_edgeR" kaydedilir
9. DEG'leri tanımlayın.
   res_edgeR$sig = as.factor(
   ifelse(res_edgeR$FDR < 0,05 & abs(res_edgeR$logFC) > 2,
   ifelse(res_edgeR$logFC > 2 ,'yukarı','aşağı'),'değil'))
   özet(res_edgeR$sig)
   ##down yukarı değil
   ##1578 15965 3121
10. Sonuç tablosunu bir dosyaya çıktısını verin.
    yazma.csv(res_edgeR, dosya = 'res_edgeR.csv')
11. Volkan arsa oluşturun.
    volkan(res_edgeR,"logFC","FDR",2,0.05)
12. Volkan arsasını kaydetmek için Dışarı Aktar'ı tıklatın.

4. "DESeq2" ile diferansiyel ifade analizi

1. R paketleri "DESeq2"yi yüklemek için Çalıştır'ı tıklatın.
   BiocManager::yükle("DESeq2")
2. R paketlerini "DESeq2" yüklemek için Çalıştır'ı tıklatın.
   kitaplık(DESeq2)
3. Gruplandırma faktörünü belirlemek için aşağıdaki R kodunu çalıştırın.
   grup <-alt dize(sütun adları(cnt),14,15)
   grup [grup %in% "01"] <- "Kanser"
   grup [grup %in% "11"] <- "Normal"
   group=faktör(grup, düzey = c("Normal","Kanser"))
4. DESeqDataSet nesnesini oluşturun.
   dds <-DESeqDataSetFromMatrix (cnt, DataFrame(grup), tasarım = ~grup)
   Dds
   ##class: DESeqDataSet
   ##dim: 20664 45
   ##metadata(1): sürüm
   ##assays(1): sayılar
   ##rownames(20664): TSPAN6 DPM1 ... RP11-274B21.13 LINC01144
   ##rowData adları(0):
   ##colnames(45): TCGA-3X-AAVA-01A-11R-A41I-07 ...
   ##colData adları(1): grup
5. Analizi gerçekleştirin.
   dds <- DESeq(dds)
6. Sonuç tablosunu oluşturun.
   res_DESeq2 <- data.frame(sonuçlar(dds))
   
   kafa(res_DESeq2)
   ## baseMean log2FoldChange lfcSE stat pvalue padj
   ##TSPAN6 4704.9243 -0,8204515 0,3371667 -2,433370 1,495899e-02 2,760180e-02
   ##DPM1 1205.9087 -0,3692497 0,1202418 -3,070894 2,134191e-03 4,838281e-03
   ##SCYL3 954.9772 0,2652530 0,2476441 1,071106 2,841218e-01 3,629059e-01
   ##C1orf112 277.7756 0,7536911 0,2518929 2,992109 2,770575e-03 6,101584e-03
   ##FGR 345.8789 -0,6423198 0,3712729 -1,730047 8,362180e-02 1,266833e-01
   ##CFH 27982.3546 -3.8761382 0,5473363 -7,081823 1,422708e-12 1,673241e-11
   NOT: Sonuç, normalleştirilmiş okuma sayısı (baseMean), log fold Change değeri (log2FoldChange), log fold değişiklik standart hatası (lfcSE), Wald istatistiği (stat), orijinal p değeri (pvalue) ve düzeltilmiş p değeri (padj) ortalamasını içeren "res_DESeq2" olarak kaydedilir.
7. DEG'leri tanımlayın.
   res_DESeq2$sig = as.factor(
   ifelse(res_DESeq2$padj < 0,05 & abs(res_DESeq2$log2FoldChange) > 2,
   ifelse(res_DESeq2$log2FoldChange > 2 ,'yukarı','aşağı'),'değil'))
   özet(res_DESeq2$sig)
   ##down yukarı değil
   ##1616 16110 2938
8. Sonuç tablosunu bir dosyaya çıktısını verin.
   yazma.csv(res_DESeq2, dosya = 'res_DESeq2.csv')
9. Volkan arsa oluşturun.
   volkan(res_DESeq2,"log2FoldChange","padj",2,0.05)
10. Volkan arsasını kaydetmek için Dışarı Aktar'ı tıklatın.

5. Venn diyagramı

1. R paketini "VennDiagram" yüklemek için Çalıştır'ı tıklatın.
   install.packages("VennDiagram")
2. R paketini "VennDiagram" yüklemek için Çalıştır'ı tıklatın.
   kütüphane (VennDiagram)
3. Düzenlenmiş DEG'lerin Venn diyagramını yapın.
   grid.newpage()
   grid.draw(venn.diagram(list(Limma=satır adları(res_
   limma[res_limma$sig=="yukarı",]),
   edgeR=satır adları(res_edgeR[res_edgeR$sig=="yukarı",]),
   DESeq2=satır adları(res_DESeq2[res_DESeq2$sig==
   "yukarı",])),
   NULL,yükseklik = 3,genişlik = 3,birim = "in",
   col="siyah",lwd=0,3,fill=c("#FF6666","#FFFF00",
   "#993366"),
   alpha=c(0,5, 0,5, 0,5),main = "Yukarı düzenlenmiş DEG'ler"))
4. Venn diyagramını kaydetmek için Dışarı Aktar'ı tıklatın.
5. Aşağı düzenlenmiş DEG'lerin Venn diyagramını yapın.
   grid.newpage()
   grid.draw(venn.diagram(list(Limma=satır adları(res_
   limma[res_limma$sig=="aşağı",]),
   edgeR=satır adları(res_edgeR[res_edgeR$sig==
   "aşağı",]),
   DESeq2=satır adları(res_DESeq2[res_DESeq2$sig=="aşağı",])),
   NULL,yükseklik = 3,genişlik = 3,birim = "in",
   col="siyah",lwd=0,3,fill=c("#FF6666","#FFFF00",
   "#993366"),
   alpha=c(0,5, 0,5, 0,5),main = "Aşağı düzenlenmiş DEG'ler"))
6. Venn diyagramını kaydetmek için Dışarı Aktar'ı tıklatın.

Sonuçlar

Volkan arsası ve Venn diyagramının özellikle kullanıldığı diferansiyel ifade analizinin sonucunu görselleştirmek için çeşitli yaklaşımlar vardır. limma, CHOL ve normal dokular arasında |logFC|≥2 ve adj ile 3323 DEG tespit etti. P.Val <0.05 eşik olarak, aralarında 1880'inIN CHOL dokularında aşağı regüle edildiği ve 1443'ü yukarı regüle edildiği (Şekil 1a). Bu arada edgeR, 1578 aşağı düzenlenmiş DEG'leri ve 3121 yukarı düzenlenmiş DEG'leri tanımladı (<...

Tartışmalar

Kanserlerde bol miktarda sapma transkriptleri RNA-seq diferansiyel analizi ile kolayca tanımlanabilir⁵. Bununla birlikte, R dili ile belirli beceriler ve uygun yöntemleri seçme kapasitesi gerektirdiğinden, RNA-seq diferansiyel ifade analizinin uygulanması genellikle kısıtlanır. Bu sorunu gidermek için, en bilinen üç yönteme (limma, EdgeR ve DESeq2) ve RNA-seq diferansiyel ifade çözümlemesi uygulamak için öğreticilere ayrıntılı bir giriş sağlıyoruz. Bu, her üç yöntemdeki ...

Malzemeler

Name	Company	Catalog Number	Comments
R		version 3.6.2	free software
Rstudio			free software