我們的生信入門課程遇到一個(gè)學(xué)員提問(wèn)，他的提問(wèn)非常清晰，有前因后果：

GSE編號(hào)：GSE144169

實(shí)驗(yàn)方法：高通量測(cè)序法

物種：智人

問(wèn)題：

我想利用這個(gè)數(shù)據(jù)集尋找哪些基因是上調(diào)的，哪些基因是下調(diào)的，所以我做了如下嘗試

1.首先我用的是小潔老師用的轉(zhuǎn)錄組代碼，但是這里面只有兩個(gè)樣本，所以采用小潔老師上課的代碼行不通

2. 后來(lái)我看附加文件中，已經(jīng)將差異基因?qū)懗晌谋靖袷?，發(fā)在附加文件下，我就直接下載了，將文件導(dǎo)入，輸入列名，發(fā)現(xiàn)沒(méi)有p.value值，但是有l(wèi)og2FC值，所以我想問(wèn)一下能不能用Huvec_Co和Huvec_Expt計(jì)算出p.value。

這個(gè)是差異基因的截圖

這個(gè)是列名的截圖

這個(gè)問(wèn)題里面涉及到兩個(gè)問(wèn)題：

1、沒(méi)有生物學(xué)重復(fù)的時(shí)候 可以使用 FC 值 即倍數(shù)變化 篩選差異基因嗎？

2、沒(méi)有生物學(xué)重復(fù)的時(shí)候 還有算法可以做差異分析嗎？進(jìn)而得到一個(gè)統(tǒng)計(jì)學(xué)顯著性Pvalue值。

先看第一個(gè)：

毫無(wú)疑問(wèn)，F(xiàn)C值 是基因在兩組樣本或者這里的一對(duì)一樣本中的倍數(shù)變化值，在早期生物信息分析里面篩選差異基因的時(shí)候，常用的指標(biāo)就是這個(gè)FC值，是可以用來(lái)篩選差異基因的，如使用閾值：FC > 2 或者FC < 1/2。

但是FC值有一個(gè)比較大的缺點(diǎn)，就是容易受到較小數(shù)值的影響（部分基因）：

如：

genei 在 A 組表達(dá)均值為 0.1，在 B 組中表達(dá)均值為 0.5，他們的差值只有 0.4，是非常微小的，但是 FC 達(dá)到了 5 倍

又：

genei 在 A 組表達(dá)均值為 2，在 B 組中表達(dá)均值為 5，他們的差值達(dá)到了 3，但是 FC 只有 2.5 倍

繪圖看一下基因表達(dá)均值與FC值的散點(diǎn)圖：

all <- read.table("GSE144169_Huvec_Co_Vs_Expt_Gene_Diff.txt.gz",header=T,sep="\t")
head(all)
all$mean <- rowMeans(all[,c(3,4)])
qplot(log10(all$mean), all$log2.fold_change., ylab="Log2FC", xlab="Log10(Exp Mean)", size=I(0.7))

且沒(méi)有一個(gè)統(tǒng)計(jì)學(xué)指標(biāo)來(lái)判斷這個(gè)基因的倍數(shù)變化是不是具有統(tǒng)計(jì)學(xué)意義。

上面這個(gè)GEO數(shù)據(jù)：https://www.ncbi.nlm./geo/query/acc.cgi?acc=GSE144169，提供的結(jié)果 是可以使用FC來(lái)進(jìn)行差異基因篩選的。

對(duì)于第二個(gè)問(wèn)題

當(dāng)然，還是有算法可以對(duì) 這種只有一個(gè) 樣本分組的差異進(jìn)行分析，如轉(zhuǎn)錄組差異分析三大算法之一的edgeR，分析代碼如下：

首先下載count矩陣：

# 加載R包
library(data.table)
library(tinyarray)

# load counts table from GEO
#urld <- "https://www.ncbi.nlm./geo/download/?type=rnaseq_counts&acc=GSE224401&format=file&file=GSE224401_raw_counts_GRCh38.p13_NCBI.tsv.gz"
acc <- "GSE144169"
urld <- "https://www.ncbi.nlm./geo/download/?type=rnaseq_counts&"
path <- paste0(urld, "acc=", acc, "&format=file&file=",acc,"_raw_counts_GRCh38.p13_NCBI.tsv.gz");
file <- paste0(acc,"_raw_counts_GRCh38.p13_NCBI.tsv.gz")
path
file
download.file(url = path, destfile = file)

tbl <- as.matrix(data.table::fread(file, header=T, colClasses="integer"), rownames=1)
dim(tbl)
ensembl_matrix=as.data.frame(tbl) 

library(AnnoProbe)
library(org.Hs.eg.db)
keytypes(org.Hs.eg.db)
e2s<-AnnotationDbi::select(org.Hs.eg.db,
 keys= rownames(ensembl_matrix),
 columns="SYMBOL",
 keytype = "ENTREZID")
head(e2s)
ids = na.omit(e2s)
ids=ids[!duplicated(ids$SYMBOL),]
ids=ids[!duplicated(ids$ENTREZID),]
head(ids)
symbol_matrix= ensembl_matrix[match(ids$ENTREZID,rownames(ensembl_matrix)),]
rownames(symbol_matrix) = ids$SYMBOL
symbol_matrix[1:4,]

接著拿到樣本表型信息：

library(AnnoProbe)
library(GEOquery)
acc
gset = getGEO(acc, destdir = '.', getGPL = F) 
pd = pData(gset[[1]]) 
pd[,1:5]
colnames(pd)
rownames(pd)
colnames(symbol_matrix)
com <- intersect(rownames(pd), colnames(symbol_matrix))
symbol_matrix <- symbol_matrix[, com]

pd=pd[com,]
pd$title 
group_list=ifelse(grepl('Control',pd$title ), 'control','case' )

# 保存數(shù)據(jù)
group_list
table(group_list)
group_list = factor(group_list,levels = c('control','case' )) 
dat=log10(edgeR::cpm(symbol_matrix)+1)
save(symbol_matrix,dat,group_list,file = 'step1-output.Rdata')

使用edgeR做一對(duì)一的差異分析：

rm(list=ls())
library(edgeR)

load("step1-output.Rdata")
ls()

# filter low expression gene
count_data <- symbol_matrix
keep <- rowSums(count_data)>0
filter_count <- count_data[keep,]

############## begin deg analysis
y <- DGEList(counts=filter_count, group=group_list)
y <- calcNormFactors(y)
dis <- 0.16 # for human
et <- exactTest(y, dispersion=dis)
res <- as.data.frame(topTags(et,dim(et)[1]))

head(res)

結(jié)果如下：

可以比較一下edgeR算法的logFC與這個(gè)數(shù)據(jù)集給出的FC值

## compare
head(res)
head(all)
com_gene <- intersect(rownames(res), all$gene)
dat_plot <- data.frame(com_gene,as.numeric(res[com_gene, "logFC"]), as.numeric(all[match(com_gene,all$gene), 5]))
head(dat_plot)
colnames(dat_plot) <- c("gene", "logFC", "edgeR_logFC")
head(dat_plot)

library(ggpubr)
p=ggscatter(dat_plot, x = "logFC", y = "edgeR_logFC",
   color = "black", shape = 21, size = 2, # Points color, shape and size
   add = "reg.line",  # Add regressin line
   add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
   conf.int = TRUE, # Add confidence interval
   cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
   cor.coeff.args = list(method = "pearson",  label.sep = "\n")
) 
p

結(jié)果如下：

edgeR 算法給出的建議：What to do if you have no replicates

他們公出了四點(diǎn)建議：但是任何一點(diǎn)都不是可以替代 有生物學(xué)重復(fù)的好方案 （千萬(wàn)要有組內(nèi)重復(fù)樣品設(shè)計(jì)）

第一條也是最好的一條，直接使用FC值來(lái)篩選基因進(jìn)行后續(xù)的研究

不要試圖去找一個(gè)統(tǒng)計(jì)學(xué)顯著性。

1. Be satisfied with a descriptive analysis, that might include an MDS plot and an analysis of fold changes. Do not attempt a significance analysis. This may be the best advice

第二條是我上面的代碼，給了一個(gè) dispersion值

2. Simply pick a reasonable dispersion value, based on your experience with similar data, and use that for exactTest or glmFit. Typical values for the common BCV (square-root dispersion) for datasets arising from well-controlled experiments are 0.4 for human data,0.1 for data on genetically identical model organisms or 0.01 for technical replicates

第三條：

第四條：

如果你也有類似的生物信息學(xué)細(xì)節(jié)問(wèn)題

不妨考慮我們的生物信息學(xué)馬拉松授課（買一得五） ，你的生物信息學(xué)入門課。，貼心答疑解惑你值得擁有！接下來(lái)的直播互動(dòng)授課開始時(shí)間是2024年12月30日， 晚上跟直播白天可以很方便的實(shí)踐操作！

如果你已經(jīng)熟悉了我們的課程，就聯(lián)系我們報(bào)名吧~
（添加好友務(wù)必備注 高?；蛘吖ぷ鲉挝?姓名+馬拉松，方便后續(xù)認(rèn)識(shí)）