关于课程中TCGA差异基因表达代码-问题

按老师给的代码,提取差异基因,遇到问题(黑体字已标出)

fdr = 0.05 
fold_change = 2
# 读取基因表达的矩阵
rawdata<-read.table("Gene_HTSeq_Counts.txt", header=TRUE, stringsAsFactors=FALSE, row.names=1)
# 针对基因的表达量进行过滤,过滤标准设置为:至少有25% 的样本,基因表达量大于2
quant <- apply(rawdata,1,quantile,0.75)    
keep <- which((quant >= 2) == 1)
rawdata <- rawdata[keep,]
dim(rawdata)
#基于样本的编号,判定是癌症还是正常样本,第14位1为癌症,0为正常组织。
group <- factor(substr(colnames(rawdata),14,14))    
summary(group)
  0   1 
473  41 
# 获得基因名称列表
genes<-rownames(rawdata)
# 制作DEGlist对象
y <- DGEList(counts=rawdata, genes=genes, group=group)
nrow(y)
# TMM 标准化
y <- calcNormFactors(y)
# 分布估计
y <- estimateCommonDisp(y, verbose=TRUE)
y <- estimateTagwiseDisp(y)
# 差异分析检验
et <- exactTest(y)           
# 针对FDR, fold change 对基因的显著性表达进行“多重检验”
et$table$FDR <- p.adjust(et$table$PValue, method="BH")
summary(de <- decideTestsDGE(et, adjust.method="BH", p.value=fdr, lfc=log2(fold_change)))
        1-0
Down    2613
NotSig 12337
Up      2413     

######################################1为正常,0为癌症,此处看起来是正常/癌症组织,那么上调或者下调是不是与我们理解的刚好相反?,那么这里以及下面的代码怎么改?

#标识基因的上下调情况
et$table$regulate = as.factor(ifelse(et$table$FDR < fdr & abs(et$table$logFC) >=log2(fold_change), ifelse(et$table$logFC>log2(fold_change),'Up','Down'),'Normal'))


请先 登录 后评论

1 个回答

omicsgene - 生物信息
擅长:重测序,遗传进化,转录组,GWAS

上下调关系:可以查看最后的结果,也就是基因表达量与fold_change的表格,哪一组当中基因整体表达量是否高于另一组当中的表达量,自然就知道上下调关系了; 

这个代码不是很重要,自己check一下结果数据就知道了; 

请先 登录 后评论