100個(gè)GEO基因表達(dá)芯片或轉(zhuǎn)錄組數(shù)據(jù)處理

寫在前邊

雖然現(xiàn)在是高通量測(cè)序的時(shí)代，但是GEO、ArrayExpress等數(shù)據(jù)庫(kù)儲(chǔ)存并公開大量的基因表達(dá)芯片數(shù)據(jù)，還是會(huì)有大量的需求去處理芯片數(shù)據(jù)，并且建模或驗(yàn)證自己所研究基因的表達(dá)情況，芯片數(shù)據(jù)的處理也可能是大部分剛學(xué)生信的道友入門R語(yǔ)言數(shù)據(jù)處理的第一次實(shí)戰(zhàn)，因此準(zhǔn)備更新100個(gè)基因表達(dá)芯片或轉(zhuǎn)錄組高通量數(shù)據(jù)的處理。

數(shù)據(jù)信息檢索

可以看到GSE28623是 芯片數(shù)據(jù)，因此可以使用GEOquery包下臨床信息，然后從網(wǎng)頁(yè)下載 原始的基因表達(dá)數(shù)據(jù)用 R 標(biāo)準(zhǔn)化處理

image.png

image.png

使用GEOquery包下載 臨床 數(shù)據(jù)

BiocManager::install('ScienceAdvances/Canton')
Canton::using(using, tidyverse, GEOquery, magrittr, data.table, AnnoProbe, clusterProfiler, org.Hs.eg.db, org.Mm.eg.db,ggdendro,ComplexHeatmap)

注：using作用是一次性加載多個(gè)R包，不用寫雙引號(hào)，并且不在屏幕上打印包的加載信息

處理表型數(shù)據(jù)

這部分是很關(guān)鍵的，可以篩選一下分組表型信息，只保留自己需要的樣本，作為后續(xù)分析的樣本（根據(jù)自己的研究目的篩選符合要求的樣本）

geo_accession <- "GSE28623"
eSet <- GEOquery::getGEO(geo_accession, AnnotGPL = F, getGPL = F)

pdata=Biobase::pData(eSet[[1]]) %>% 
    dplyr::mutate(
        Sample = geo_accession,
        Group =  case_when(str_detect(`group:ch1`,"latently")~"LTBI",str_detect(`group:ch1`,"non-infected")~"HC",str_detect(`group:ch1`,"tuberculosis")~"TB",TRUE~NA),
        Sex = str_to_title(`gender:ch1`)
    ) %>%
    dplyr::select(Sample, Group,Sex, everything())

處理表達(dá)譜數(shù)據(jù)

讀取原始數(shù)據(jù)

raw_dir <- "GSE28623_RAW"
x <- limma::read.maimages(
    files = list.files(raw_dir, pattern = "^GSM", full.names = T),
    source = "agilent",
    green.only = TRUE,
    other.columns = "gIsWellAboveBG"
)

使用 limma 標(biāo)準(zhǔn)化處理芯片數(shù)據(jù)

y <- limma::backgroundCorrect(x, method = "normexp") %>% limma::normalizeBetweenArrays(method = "quantile")

刪除 Control 探針，沒有基因名字的探針，重復(fù)的探針

Control <- y$genes$ControlType == 1L
NoSymbol <- is.na(y$genes$GeneName)
Isdup <- duplicated(y$genes$GeneName)
yfilt <- y[!Control & !NoSymbol & !Isdup, ]

獲取表達(dá)量數(shù)據(jù)

fdata <- yfilt@.Data[[1]]

刪除樣本名中的路徑字符串 GSE28623/GSM709520_251485039549_1_4.txt -> GSM709520

colnames(fdata) %<>% str_extract("GSM\\d*")

給數(shù)據(jù)框添加基因名作為行名

rownames(fdata) <- yfilt$genes$GeneName

數(shù)據(jù)質(zhì)控

qcplot為自定義函數(shù)，作用是繪制用于質(zhì)控判斷的圖，如PCA、top20基因熱圖、<font style="color:rgba(0, 0, 0, 0.85);">樹狀圖</font>，PCA圖，可以三組區(qū)別區(qū)別不明顯，可以考慮直接用作者處理好的數(shù)據(jù)雖然后負(fù)數(shù)值

qcplot(
  data=fdata, 
  group=pdata$Group, 
  w=18,h=18
) 

image.png

保存數(shù)據(jù)

common_samples <- base::intersect(colnames(fdata),pdata$Sample)
fdata=fdata[,common_samples]

fwrite(as.data.table(fdata,keep.rownames="Feature"), file = stringr::str_glue("{geo_accession}_fdata.csv.gz"))
pdata %<>% dplyr::filter(Sample %in% common_samples)
fwrite(pdata, file = stringr::str_glue("{geo_accession}_pdata.csv"))

Reference

http://www.itdecent.cn/p/19ff109b6409