100個GEO基因表達芯片或轉錄組數據處理
寫在前邊
雖然現在是高通量測序的時代，但是GEO、ArrayExpress等數據庫儲存并公開大量的基因表達芯片數據，還是會有大量的需求去處理芯片數據，并且建模或驗證自己所研究基因的表達情況，芯片數據的處理也可能是大部分剛學生信的道友入門R語言數據處理的第一次實戰，因此準備更新100個基因表達芯片或轉錄組高通量數據的處理。
數據信息檢索
可以看到GSE24807是基因表達芯片數據，因此可以使用GEOquery包下
使用GEOquery包下載數據
remotes::install_github('ScienceAdvances/using')
using::using(tidyverse, GEOquery, magrittr, data.table, AnnoProbe, clusterProfiler, org.Hs.eg.db, org.Mm.eg.db)
注：using作用是一次性加載多個R包，不用寫雙引號，并且不在屏幕上打印包的加載信息
因為文件太大，在R內下載失敗，可通過圖片中的方法下載文件，GEOquery::getGEO直接讀取本地的文件。
geo_accession <- "GSE24807"
eSet <- getGEO(filename=stringr::str_glue('{geo_accession}_series_matrix.txt.gz'), AnnotGPL = F, getGPL = F)
gpl <- eSet@annotation
處理表型數據
這部分是很關鍵的，可以篩選一下分組表型信息，只保留自己需要的樣本，作為后續分析的樣本（根據自己的研究目的篩選符合要求的樣本）
pdata <- pData(eSet)
pdata %<>%dplyr::mutate(Sample = geo_accession,Group = case_when(`disease state:ch1`=='non-alcoholic steatohepatitis (NASH)'~'NASH',`disease state:ch1`=='normal (control)'~'Control',TRUE~NA)) %>%drop_na(Group) %>% dplyr::select(Sample,Group,everything())
處理表達譜數據
數據大小不大于50不需要取log
exprs_mtx <- exprs(eSet)
if(max(exprs_mtx, na.rm = TRUE)<50 | min(exprs_mtx, na.rm = TRUE)<0){message("基因表達最大值小于50或者最小值小于0不需要log轉化")
}else {message("基因表達最大值大于50需要log轉化")exprs_mtx <- log2(exprs_mtx+1)
}
probe_exprs <- as.data.table(exprs_mtx, keep.rownames = "ProbeID")
探針與基因Symbol對應關系
從https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL2895下載注釋文件GPL2895.annot.gz，從中獲取探針與GeneID對應關系
idmaps <- function(ann_file, ProbeID = "ID", Feature = "Symbol", skip = 229, pattern = "control") {temp <- fread(ann_file, skip = skip, nThread = 8)vars <- c(ProbeID, Feature)temp <- temp[, ..vars]data.table::setnames(temp, c("ProbeID", "Feature"))temp %<>% tidyr::separate_longer_delim(cols=Feature, delim=" /// ") %>% data.table::as.data.table()temp %<>% tidyr::separate_longer_delim(cols=Feature, delim="///") %>% data.table::as.data.table()temp <- temp[!is.null(Feature), ][!is.na(Feature), ][Feature != "", ][Feature != "---", ][!stringr::str_detect(string = Feature, pattern = pattern), ]return(as.data.frame(temp))
}
probe2symbol <- idmaps('GPL2895.annot.gz', Feature = "Gene symbol", skip = 27)
ID轉換
把表達矩陣中的探針名轉換為基因名；transid是我寫的一個R函數，有需要可以加我（18983376561），進入交流群
fdata <- transid(probe2symbol, probe_exprs)
保存數據
common_samples <- base::intersect(colnames(fdata),pdata$Sample)
fdata %<>% select(all_of(c("Feature",common_samples)))
fwrite(fdata, file = stringr::str_glue("{geo_accession}_{gpl}_fdata.csv.gz"))
pdata %<>% dplyr::filter(Sample %in% common_samples)
fwrite(pdata, file = stringr::str_glue("{geo_accession}_{gpl}_pdata.csv"))
