setwd("/Users/Matsui/hgc2017/day4") #発現差異解析 library(data.table) input_file <- "data/TCGA_BRCA_mRNAseq.tsv" expr <- as.matrix(fread(input_file, header = TRUE, drop = 1,sep="\t")) rownames(expr) <- as.vector(as.matrix(fread(input_file,header=TRUE,select=1,sep="\t"))) #?という遺伝子シンボルの行はデータから除く select_idx <- which("?" != rownames(expr)) expr <- expr[select_idx,] #tumor vs normalのサンプルのみに絞る barcode <- colnames(expr) samplecode <- as.numeric(substr(barcode,14,15)) tissuecode <- substr(barcode,9,12) cidx_tumor <- which(samplecode == 1) cidx_normal <- which(samplecode == 11) tumor <- tissuecode[cidx_tumor] normal <- tissuecode[cidx_normal] common_code <- intersect(tumor,normal) tumor_idx <- match(common_code,tumor) normal_idx <- match(common_code,normal) cidx_tumor <- cidx_tumor[tumor_idx] cidx_normal <- cidx_normal[normal_idx] #ケース(腫瘍組織)とコントロール(対応付きの正常組織)をexprから取り出して #caseとctrに代入しよう case <- expr[,cidx_tumor] ctr <- expr[,cidx_normal] #caseとcontrol群でゼロカウントの発現が極端に多いサンプルは除いておく #caseとcontrolそれぞれで、遺伝子ごとにゼロカウントの発現データが何サンプルあるかを数えよう nzero_case <- apply(case,1,function(x)sum(x == 0)) nzero_ctr <- apply(ctr,1,function(x)sum(x == 0)) #ゼロカウントが含まれていない行の添え字を計算しましょう nonzero_case <- which(nzero_case == 0) nonzero_ctr <- which(nzero_ctr == 0) #ケースとコントロールで共通した行の添え字のみを取り出す common_ridx <- intersect(nonzero_case,nonzero_ctr) #common_idxを用いて、caseとctrをフィルタリング case <- case[common_ridx,] ctr <- ctr[common_ridx,] #一つのデータに列方向につなげる expr <- cbind(case,ctr) #重複遺伝子の発現プロファイルを一つにまとめる #分散最大の遺伝子プロファイルを選択する unique_gene <- unique(rownames(expr)) unique_expr <- matrix(NA,nrow = length(unique_gene), ncol = ncol(expr)) rownames(unique_expr) <- unique_gene colnames(unique_expr) <- colnames(expr) for(i in seq_along(unique_gene)){ idx <- which(unique_gene[i] == rownames(expr)) score <- apply(expr[idx,,drop = F], 1, var) select_idx <- idx[which.max(score)] unique_expr[i,] <- expr[select_idx,] cat(i,"\n") } expr <- unique_expr case <- expr[,match(colnames(case),colnames(expr))] ctr <- expr[,match(colnames(ctr),colnames(expr))] #log fold change case_mean <- apply(case,1,mean) ctr_mean <- apply(ctr,1,mean) log2FC <- log2((case_mean + 1) / (ctr_mean + 1)) #logFCでがん特的な発現上昇・減少である上位100遺伝子 ntop <- 100 toplogFC_up <- log2FC[order(log2FC,decreasing = TRUE)][1:ntop] names(toplogFC_up) toplogFC_down <- log2FC[order(log2FC)][1:ntop] names(toplogFC_down) #t検定による発現差異解析を実行してみよう #対数変換しておく logcase <- log2(case + 1) #logcase <- log2(case + 1) logctr <- log2(ctr + 1) #logctr <- log2(ctr + 1) #ケース群特異的に発現上昇あるいは減少している遺伝子群を同定してみよう。 pvals_up <- pvals_down <- rep(NA,length(nrow(expr))) #両側検定の結果も用意しておく pvals_twoside <- rep(NA,length(nrow(expr))) tstat <- rep(NA, length(nrow(expr))) #t統計量用の変数 n <- nrow(expr) for(i in 1:n){ t_test <- t.test(logcase[i,],logctr[i,],paired = TRUE, alternative = "greater") pvals_up[i] <- t_test$p.value t_test <- t.test(logcase[i,],logctr[i,],paired = TRUE, alternative = "less") pvals_down[i] <- t_test$p.value t_test <- t.test(logcase[i,],logctr[i,],paired = TRUE, alternative = "two.sided") pvals_twoside[i] <- t_test$p.value tstat[i] <- t_test$statistic cat(i , "\n") } #多重検定補正 qvals_up <- p.adjust(pvals_up,method = "BH") qvals_down <- p.adjust(pvals_down,method = "BH") qvals_twoside <- p.adjust(pvals_twoside, method = "BH") #-log10(FDR)とlogFCのプロット pdf("data/ttest_result.pdf",10,10) plot(log2FC,-log10(qvals_twoside),pch = 20,cex=0.15,ylab="-log10(FDR)",xlab="log fold change") abline(v = -1.5, col = "red") abline(v = 1.5, col = "red") abline(h = 2, col = "blue") #FDR <= 0.01 grid() #|log2FC| >= 2 かつ FDR <= 0.01 の遺伝子群を表示 disp_gene <- which(abs(log2FC) >= 2 & qvals_twoside <= 0.01) text(log2FC[disp_gene],-log10(qvals_twoside)[disp_gene],labels = rownames(expr)[disp_gene],cex=0.4) dev.off() #有意な遺伝子群の名前を取り出す upgenes <- rownames(case)[qvals_up <= 0.05] downgenes <- rownames(ctr)[qvals_down <= 0.05] #limmaによる解析 library(limma) library(edgeR) dge <- DGEList(expr) logCPM <- cpm(dge,log=TRUE,prior.count = 1) group <- c(rep("T",ncol(case)),rep("N",ncol(ctr))) group <- factor(group,levels = c("N","T")) design <- model.matrix(~group) fit <- lmFit(logCPM,design) fit <- eBayes(fit,trend=TRUE) topTable(fit,coef = ncol(design)) #まとめて結果をファイルに出力 #出力用の結果(発現上昇・減少の両方の結果ファイルを用意する) output <- data.frame(gene=rownames(expr),adj.p.up = qvals_up, adj.p.down = qvals_down, log2FC = log2FC, tstat = tstat) #adj.p.upで昇順に並べ替えてエクセルで保存(tsvかcsvで拡張子だけxlsにする) output <- output[order(output$adj.p.up),] write.table(output,file="data/DGE_ttest.xls",row.names=FALSE,col.names=TRUE,quote=FALSE,sep="\t") #遺伝子セット解析 #遺伝子セットの読み込み gmt <- scan("data/c5.all.v6.0.symbols.gmt.txt",what = "character",sep="\n") gmt <- strsplit(gmt,"\t") setname <- sapply(gmt,function(x)x[1]) seturl <- sapply(gmt,function(x)x[2]) setlist <- sapply(gmt,function(x)x[-c(1,2)]) names(setlist) <- setname #fisherの正確度検定 #upgeneに対して実行 ngene <- nrow(expr) all_genes <- rownames(expr) no_upgenes <- setdiff(all_genes,upgenes) FETup_pval <- rep(NA, length(setlist)) for(i in seq_along(setlist)){ refset <- setlist[[i]] A <- length(intersect(refset,upgenes)) B <- length(setdiff(upgenes,refset)) C <- length(intersect(refset,no_upgenes)) D <- length(setdiff(no_upgenes,refset)) ctb <- matrix(c(A,B,C,D),nrow=2) FET <- fisher.test(ctb,alternative = "greater") FETup_pval[i] <- FET$p.value cat(i , "\n") } FETup_qval <- p.adjust(FETup_pval,method="BH") #downgeneに対して実行 no_downgenes <- setdiff(all_genes,downgenes) FETdown_pval <- rep(NA, length(setlist)) for(i in seq_along(setlist)){ refset <- setlist[[i]] A <- length(intersect(refset,downgenes)) B <- length(setdiff(downgenes,refset)) C <- length(intersect(refset,no_downgenes)) D <- length(setdiff(no_downgenes,refset)) ctb <- matrix(c(A,B,C,D),nrow=2) FET <- fisher.test(ctb,alternative = "greater") FETdown_pval[i] <- FET$p.value cat(i , "\n") } FETdown_qval <- p.adjust(FETdown_pval,method="BH") #結果ファイルをupgeneとdowngeneで用意 FETup_ouput <- data.frame(pathway=setname,url=seturl,adj.p = FETup_qval) FETdown_ouput <- data.frame(pathway=setname,url=seturl,adj.p = FETdown_qval) #FDRで昇順に並び替えて保存 FETup_ouput <- FETup_ouput[order(FETup_ouput$adj.p),] FETdown_ouput <- FETdown_ouput[order(FETdown_ouput$adj.p),] write.table(FETup_ouput,file="data/FETup.xls",row.names=FALSE,col.names=TRUE, sep="\t",quote=F) write.table(FETdown_ouput,file="data/FETdown.xls",row.names=FALSE,col.names=TRUE, sep="\t",quote=F) #GSEA #入力データは #1.tstat #2.遺伝子セット library(fgsea) library(ggplot2) names(tstat) <- rownames(expr) geneset <- gmtPathways("data/c5.all.v6.0.symbols.gmt.txt") fgsea_result <- fgsea(geneset,tstat,nperm = 10000) #FDRで昇順に並べ替える fgsea_result <- fgsea_result[order(fgsea_result$padj),] output <- as.matrix(fgsea_result) write.table(output,file="data/GSEA_result.xls",row.names=F,col.names=T,quote=F,sep="\t") signif_result <- fgsea_result[signif_idx,] p <- plotEnrichment(geneset[["GO_DNA_REPAIR"]],tstat) + labs(title="GO_DNA_REPAIR") ggsave(p,filename = "data/GSEA_ESplot.pdf") ntop <- 10 top_pathway_up <- fgsea_result[fgsea_result$ES > 0 & fgsea_result$padj <= 0.05,] top_pathway_up <- top_pathway_up[order(top_pathway_up$padj),] top_pathway_down <- fgsea_result[fgsea_result$ES < 0 & fgsea_result$padj <= 0.05,] top_pathway_down <- top_pathway_down[order(top_pathway_down$padj,decreasing=TRUE),] top_pathway <- rbind(top_pathway_up[1:ntop,], top_pathway_down[1:ntop,]) p <- plotGseaTable(geneset[top_pathway$pathway],tstat,fgsea_result) ggsave(p,filename="data/GSEA_table.pdf",width = 20,height = 20)