search

ixxmu / mp_duty

抓取网络文章到github issues保存
https://archives.duty-machine.now.sh/
122 stars 30 forks source link

送你一篇TCGA数据挖掘文章 #27

Closed ixxmu closed 4 years ago

ixxmu commented 4 years ago

https://mp.weixin.qq.com/s/t48BBhj19LZtLNcuHM3QEQ

github-actions[bot] commented 4 years ago

送你一篇TCGA数据挖掘文章 by 生信技能树



下载样本信息

首先从UCSC Xena数据库下载样本信息

UCSC Xena网址：https://xena.ucsc.edu/public-hubs/

根据上图的链接可以下载到样本信息，文件为

“TCGA-BRCA.GDC_phenotype.tsv.gz”

通过上图得到的三个重要信息：

  1. 数据下载链接

  2. 样本数：1217

  3. count值计算方法：log2(count+1)

下载好样本信息及表达矩阵的数据之后，我们就可以开始处理数据了。


挑选出TNBC样本

首先根据样本信息找到三阴性乳腺癌的样本。

载入下载好的样本信息文件

phenotype_file <- read.table('TCGA-BRCA.GDC_phenotype.tsv.gz',header = T,sep = '\t',quote = '')

## 检查一下表头，其实Xena上有两个样本信息的文件，选择'TCGA-BRCA.GDC_phenotype.tsv.gz'的原因就在于另一个样本信息文件所包含的内容过少。大家可以下载下来看一下。

(phenotype_colnames <- asN.data.frame(colnames(phenotype_file)))

## 三阴性乳腺癌的患者不表达ER，PR，Her2，所以先检查一下样本信息中的这三列

table(phenotype_file$breast_carcinoma_estrogen_receptor_status) ## 雌激素受体(ER)

table(phenotype_file$breast_carcinoma_progesterone_receptor_status) ## 黄体酮受体(PR)

table(phenotype_file$lab_proc_her2_neu_immunohistochemistry_receptor_status) ## 人类表皮生长因子受体(HER2)

## 取出含有'receptor_status'信息的列

colnames_num <- grep('receptor_status',colnames(phenotype_file))

phenotype_colnames <- colnames(phenotype_file)[colnames_num]

eph <- phenotype_file[,colnames_num[1:3]]


## 之后的代码用到apply函数，会用跳过这部分就好

## apply函数需要三个参数，第一个参数是matrix

## 第二个参数如果是1，说明是按行取；第二个参数如果是2，说明是按列取

## 第三个参数是方法

## example

## > x <- cbind(x1 = 3, x2 = c(4:1))

## > dimnames(x)[[1]] <- letters[1:4]

## > x

## x1 x2

## a  3  4

## b  3  3

## c  3  2

## d  3  1

## > apply(x, 2, function(x) sum(x =='3'))

## x1 x2 

## 4  1 

## > apply(x, 2, function(x) sum(x =='3'))

## x1 x2 

## 4  1


## 根据ER,PR,HER2将样本分组

tnbc_rownum <- apply(eph, 1, function(x) sum(x =='Negative'))

tnbc_sample <- phenotype_file[tnbc_rownum == 3, 1]

save(tnbc_s,file = 'tnbc_sample.Rdata')


到这里，我们就从1217个样本中挑出了118个tnbc样本，接下来就可以用在表达矩阵中取出这些样本了

从Xena下载到的矩阵不是可以直接用的，我们要先把它处理一下

library(data.table)

a=fread('TCGA-BRCA.htseq_counts.tsv',sep = '\t',header = T)

a=as.data.frame(a)

a[1:4,1:4] 

rownames(a)=a[,1]

a=a[,-1]

genes=rownames(a)

a[1:4,1:4]

## 在数据的介绍页面上我们已经得知了数据的计算方法现在我们只要把它还原回去就可以了

a=2^a-1

a[1:4,1:4] 

## 接下来我们就要取出表达数据，但是我们只要这118个tnbc样本中成对的样本，## 即，同一个tnbc病人既有正常样本，又有肿瘤样本的表达信息

load('tnbc_s.Rdata')

tnbc_p=substring(tnbc_s,1,12)

all_p=substring(colnames(a),1,12)

paired_p=names(table(all_p)[table(all_p)==2])

## 这样挑选过后，符合要求的样本就只有9个了

need_p=intersect(tnbc_p,paired_p)

exprSet=a[,all_p %in% need_p]

tmp=apply(exprSet,1,function(x){

  sum(x==0) < 10

})

exprSet=exprSet[tmp,]

save(exprSet,file = 'tnbc_paired_exprSet.Rdata')


差异分析

## 接下来做差异分析

Rdata_dir='.'

load( file = 

        file.path(Rdata_dir,'tnbc_paired_exprSet.Rdata')

)

dim(exprSet)

## 根据TCGA样本的命名可以区分正常组织和肿瘤样本的测序结果

group_list=factor(ifelse(as.numeric(substr(colnames(exprSet),14,15)) < 10,'tumor','normal'))

table(group_list)

## 对表达差异结果数据取整

exprSet=ceiling(exprSet)

这个时候exprSet结果就可以进行差异分析了


## 方法一：DESeq2

if(T){

  library(DESeq2)

  

  (colData <- data.frame(row.names=colnames(exprSet), 

                         group_list=group_list) )

  dds <- DESeqDataSetFromMatrix(countData = exprSet,

                                colData = colData,

                                design = ~ group_list)

  tmp_f=file.path(Rdata_dir,'TCGA-KIRC-miRNA-DESeq2-dds.Rdata')

  if(!file.exists(tmp_f)){

    dds <- DESeq(dds)

    save(dds,file = tmp_f)

  }

  load(file = tmp_f)

  res <- results(dds, 

                 contrast=c("group_list","tumor","normal"))

  resOrdered <- res[order(res$padj),]

  head(resOrdered)

  DEG =as.data.frame(resOrdered)

  DESeq2_DEG = na.omit(DEG)

  

  nrDEG=DESeq2_DEG[,c(2,6)]

  colnames(nrDEG)=c('log2FoldChange','pvalue')  

}


## 方法二：edgeR

if(T){

  library(edgeR)

  d <- DGEList(counts=exprSet,group=factor(group_list))

  keep <- rowSums(cpm(d)>1) >= 2

  table(keep)

  d <- d[keep, , keep.lib.sizes=FALSE]

  d$samples$lib.size <- colSums(d$counts)

  d <- calcNormFactors(d)

  d$samples

  dge=d

  design <- model.matrix(~0+factor(group_list))

  rownames(design)<-colnames(dge)

  colnames(design)<-levels(factor(group_list))

  dge=d

  dge <- estimateGLMCommonDisp(dge,design)

  dge <- estimateGLMTrendedDisp(dge, design)

  dge <- estimateGLMTagwiseDisp(dge, design)

  

  fit <- glmFit(dge, design)

  # https://www.biostars.org/p/110861/

  lrt <- glmLRT(fit,  contrast=c(-1,1)) 

  nrDEG=topTags(lrt, n=nrow(dge))

  nrDEG=as.data.frame(nrDEG)

  head(nrDEG)

  edgeR_DEG =nrDEG 

  nrDEG=edgeR_DEG[,c(1,5)]

  colnames(nrDEG)=c('log2FoldChange','pvalue') 

}


## 方法三：limma

if(T){

  suppressMessages(library(limma))

  design <- model.matrix(~0+factor(group_list))

  colnames(design)=levels(factor(group_list))

  rownames(design)=colnames(exprSet)

  design

  

  dge <- DGEList(counts=exprSet)

  dge <- calcNormFactors(dge)

  logCPM <- cpm(dge, log=TRUE, prior.count=3)

  

  v <- voom(dge,design,plot=TRUE, normalize="quantile")

  fit <- lmFit(v, design)

  

  group_list

  cont.matrix=makeContrasts(contrasts=c('tumor-normal'),levels = design)

  fit2=contrasts.fit(fit,cont.matrix)

  fit2=eBayes(fit2)

  

  tempOutput = topTable(fit2, coef='tumor-normal', n=Inf)

  DEG_limma_voom = na.omit(tempOutput)

  head(DEG_limma_voom)

  nrDEG=DEG_limma_voom[,c(1,4)]

  colnames(nrDEG)=c('log2FoldChange','pvalue') 

}


## 比较一下这三个差异分析的结果

nrDEG1=DEG_limma_voom[,c(1,4)]

colnames(nrDEG1)=c('log2FoldChange','pvalue') 

nrDEG2=edgeR_DEG[,c(1,5)]

colnames(nrDEG2)=c('log2FoldChange','pvalue') 

nrDEG3=DESeq2_DEG[,c(2,6)]

colnames(nrDEG3)=c('log2FoldChange','pvalue')  

mi=unique(c(rownames(nrDEG1),rownames(nrDEG1),rownames(nrDEG1)))

lf=data.frame(limma=nrDEG1[mi,1],

              edgeR=nrDEG2[mi,1],

              DESeq2=nrDEG3[mi,1])

cor(na.omit(lf))

##            limma     edgeR    DESeq2

## limma  1.0000000 0.9134287 0.9402436

## edgeR  0.9134287 1.0000000 0.9568595

## DESeq2 0.9402436 0.9568595 1.0000000


注释

library( "clusterProfiler" )

library( "org.Hs.eg.db" )

nrDEG = DEG_limma_voom


colnames(DESeq2_DEG)[2] = 'logFC'

colnames(DESeq2_DEG)[5] = 'P.Value'

nrDEG = DESeq2_DEG


colnames(edgeR_DEG)[1] = 'logFC'

colnames(edgeR_DEG)[4] = 'P.Value'

nrDEG = edgeR_DEG


logFC_cutoff <- with( nrDEG, mean( abs( logFC ) ) + 2 * sd( abs( logFC ) ) )

logFC_cutoff

logFC_cutoff = 1.2

nrDEG$change = as.factor( ifelse( nrDEG$P.Value < 0.05 & abs(nrDEG$logFC) > logFC_cutoff,

                                  ifelse( nrDEG$logFC > logFC_cutoff , 'UP', 'DOWN' ), 'NOT' ) )

table(nrDEG$change)

keytypes(org.Hs.eg.db)

library("stringr")

rownames( nrDEG ) <- str_sub(rownames( nrDEG ), start = 1, end = 15)

nrDEG$ENSEMBL <- rownames( nrDEG )

df <- bitr( rownames( nrDEG ), fromType = "ENSEMBL", toType = c( "ENTREZID" ), OrgDb = org.Hs.eg.db )

head( df )

{

  nrDEG$SYMBOL = rownames( nrDEG )

  nrDEG = merge( nrDEG, df, by='ENSEMBL' )

}

head( nrDEG )

{

  gene_up = nrDEG[ nrDEG$change == 'UP', 'ENTREZID' ] 

  gene_down = nrDEG[ nrDEG$change == 'DOWN', 'ENTREZID' ]

  gene_diff = c( gene_up, gene_down )

  gene_all = as.character(nrDEG[ ,'ENTREZID'] )

}

{

  geneList = nrDEG$logFC

  names( geneList ) = nrDEG$ENTREZID

  geneList = sort( geneList, decreasing = T )

}

{

  ## KEGG pathway analysis

  kk.up <- enrichKEGG(   gene          =  gene_up    ,

                         organism      =  'hsa'      ,

                         universe      =  gene_all   ,

                         pvalueCutoff  =  0.8        ,

                         qvalueCutoff  =  0.8        )

  kk.down <- enrichKEGG( gene          =  gene_down  ,

                         organism      =  'hsa'      ,

                         universe      =  gene_all   ,

                         pvalueCutoff  =  0.05       )

}

head( kk.up )[ ,1:6 ]

head( kk.down )[ ,1:6 ]


library( "ggplot2" )

{

  kegg_down_dt <- as.data.frame( kk.down )

  kegg_up_dt <- as.data.frame( kk.up )

  down_kegg <- kegg_down_dt[ kegg_down_dt$pvalue < 0.05, ]

  down_kegg$group = -1

  up_kegg <- kegg_up_dt[ kegg_up_dt$pvalue < 0.05, ]

  up_kegg$group = 1

  

  dat = rbind( up_kegg, down_kegg )

  dat$pvalue = -log10( dat$pvalue )

  dat$pvalue = dat$pvalue * dat$group

  

  dat = dat[ order( dat$pvalue, decreasing = F ), ]

  

  g_kegg <- ggplot( dat, 

                    aes(x = reorder( Description, order( pvalue, decreasing=F ) ), y = pvalue, fill = group)) + 

    geom_bar( stat = "identity" ) + 

    scale_fill_gradient( low = "blue", high = "red", guide = FALSE ) + 

    scale_x_discrete( name = "Pathway names" ) +

    scale_y_continuous( name = "log10P-value" ) +

    coord_flip() + theme_bw() + theme( plot.title = element_text( hjust = 0.5 ) ) +

    ggtitle( "Pathway Enrichment" ) 

  print( g_kegg )

  ggsave( g_kegg, filename = 'kegg_up_down_edgeR.png' )

}

KEGG通路几乎没有差异，说明三个方法找到的差异基因很一致


# 挑选显著表达差异的基因热图可视化看看效果

load(file = 'tnbc_paired_exprSet.Rdata')

load('TCGA-BRCA-mRNA-DEG_results.Rdata')

exprSet[1:4,1:4]

exprSet=log2(edgeR::cpm(exprSet)+1)

pheatmap::pheatmap(exprSet[rownames(head(DESeq2_DEG)),])


## 挑选指定基因看它在同一个病人的配对表达情况

dim(exprSet)

group_list=factor(ifelse(as.numeric(substr(colnames(exprSet),14,15)) < 10,'tumor','normal'))

table(group_list)

library(ggpubr)

dat=data.frame(g=group_list,

               p=colnames(exprSet),

               v=as.numeric(exprSet[rownames(DESeq2_DEG[1,]),])) 

ggpaired(dat, x = "g", y = "v",

         color = "g", line.color = "gray", line.size = 0.4,

         palette = "npg")+stat_compare_means()

第一个基因相当有差异性


■   ■   ■

生信基础知识大全系列：生信基础知识100讲   

史上最强的生信自学环境准备课来啦！！ 7次改版，11节课程，14K的讲稿，30个夜晚打磨，100页PPT的课程。   

如果需要组装自己的服务器；代办生物信息学服务器

如果需要帮忙下载海外数据(GEO/TCGA/GTEx等等)，点我？

如果需要线下辅导及培训，看招学徒 

如果需要个人电脑：个人计算机推荐

如果需要置办生物信息学书籍，看：生信人必备书单

如果需要实习岗位：实习职位发布

如果需要售后：点我

如果需要入门资料大全：点我


点击下面的阅读原文直达我的github