寫在前邊
雖然現(xiàn)在是高通量測序的時(shí)代，但是GEO、ArrayExpress等數(shù)據(jù)庫儲存并公開大量的基因表達(dá)芯片數(shù)據(jù)，還是會(huì)有大量的需求去處理芯片數(shù)據(jù)，并且建?；蝌?yàn)證自己所研究基因的表達(dá)情況，芯片數(shù)據(jù)的處理也可能是大部分剛學(xué)生信的道友入門R語言數(shù)據(jù)處理的第一次實(shí)戰(zhàn)，因此準(zhǔn)備更新100個(gè)基因表達(dá)芯片或轉(zhuǎn)錄組高通量數(shù)據(jù)的處理。
數(shù)據(jù)信息檢索
可以看到GSE36376是基因表達(dá)芯片數(shù)據(jù)，因此可以使用GEOquery包下
使用GEOquery包下載數(shù)據(jù)
remotes::install_github('ScienceAdvances/using')
using::using(using, tidyverse, GEOquery, magrittr, data.table, AnnoProbe, clusterProfiler, org.Hs.eg.db, org.Mm.eg.db,ggdendro,ComplexHeatmap)
注：using作用是一次性加載多個(gè)R包，不用寫雙引號，并且不在屏幕上打印包的加載信息
因?yàn)槲募螅赗內(nèi)下載失敗，可通過圖片中的方法下載文件，GEOquery::getGEO直接讀取本地的文件。
geo_accession <- "GSE36376"
eSet <- getGEO(filename=stringr::str_glue('{geo_accession}_series_matrix.txt.gz'), AnnotGPL = F, getGPL = F)
gpl <- eSet@annotation
處理表型數(shù)據(jù)
下載表型信息
這部分是很關(guān)鍵的，可以篩選一下分組表型信息，只保留自己需要的樣本，作為后續(xù)分析的樣本（根據(jù)自己的研究目的篩選符合要求的樣本）
data <- readxl::read_xls('GSE5281_sample_characteristics.xls') %>%
dplyr::rename(Sample=`GEO Accession:`,Tissue=`Organ Region:`) %>%
dplyr::mutate(
Group=ifelse(`Sample Description:`=='AD control','Control','AD'),
Sex=str_to_title(`Sex:`),
Age = str_remove_all(`Age:`,'(years)|(days)|( )')%>% as.numeric(),
) %>%
dplyr::filter(!is.na(Group)) %>%
dplyr::select(Sample,Group,Age,Sex,Tissue,everything())
處理表達(dá)譜數(shù)據(jù)
數(shù)據(jù)大小大于50需要取log
exprs_mtx <- exprs(eSet)
if(max(exprs_mtx, na.rm = TRUE)<50 | min(exprs_mtx, na.rm = TRUE)<0){
message("基因表達(dá)最大值小于50或者最小值小于0不需要log轉(zhuǎn)化")
}else {
message("基因表達(dá)最大值大于50需要log轉(zhuǎn)化")
exprs_mtx <- log2(exprs_mtx+1)
}
probe_exprs <- as.data.table(exprs_mtx, keep.rownames = "ProbeID")
探針與基因Symbol對應(yīng)關(guān)系
從AnnoProbe包中獲取探針與GeneID對應(yīng)關(guān)系
AnnoProbe::checkGPL(gpl)
probe2symbol <- AnnoProbe::idmap(gpl = gpl, type = "bioc", mirror = "tencent", destdir = outdir) %>%
dplyr::rename(ProbeID = probe_id, Feature = symbol) # pipe", "bioc", "soft"
ID轉(zhuǎn)換
把表達(dá)矩陣中的探針名轉(zhuǎn)換為基因名；transid是我寫的一個(gè)R函數(shù)，有需要可以聯(lián)系我，加入交流群
fdata <- transid(probe2symbol, probe_exprs)
保存數(shù)據(jù)
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"))
數(shù)據(jù)質(zhì)控
qcplot為自定義函數(shù)，作用是繪制用于質(zhì)控判斷的圖，如PCA、top20基因熱圖、樹狀圖，PCA圖可以看到兩組樣本有差異但是也有重疊，熱圖看到兩組樣本一些基因表達(dá)趨勢有差異，因?yàn)闃颖据^多觀察樹狀圖圖可能不大容易
qcplot(
data=fdata %>% column_to_rownames('Feature'),
pdata=column_to_rownames(pdata,'Sample') %>% dplyr::select(Group),
w=18,h=9
)
