setwd("/Users/Matsui/hgc2017/day4") #発現差異解析 #データの読み込み library(data.table) input_file <- "data/TCGA_BRCA_mRNAseq.tsv" #####exercise###### expr <- as.matrix(fread()) rownames(expr) <- as.vector(as.matrix(fread())) #####データを整形・フィルタリング##### # "?"という遺伝子シンボルの行はデータから除く #####exercise###### # exprの行名ベクトルに対して"?"に対応する添え字? select_idx <- which() ################## expr <- expr[select_idx,] #同一検体からの正常組織と腫瘍組織のペアのみの発現データをつくる barcode <- colnames(expr) #####exercise###### #barcodeの14文字目から15文字までがサンプルの種別(腫瘍組織、正常組織)を表している #barcodeの各文字列から14文字目から15文字目までを取り出してください #ヒント:substr関数 samplecode <- #9文字目から12文字目までが組織(患者)IDを表している #barcodeの各文字列から9文字目から12文字目までを取り出してください tissuecode <- ################## #samplecodeの中は、1が腫瘍組織, 11が正常組織を表している #データexprの列名ベクトルに対して腫瘍組織と正常組織を表す添え字を作成 #####exercise###### ##samplecodeから"1"に対応する添え字を取り出す cidx_tumor <- which() #samplecodeから"11"に対応する添え字を取り出す cidx_normal <- which() ################## #腫瘍組織・正常組織がそろっているサンプルIDを取り出す tumor <- tissuecode[cidx_tumor] normal <- tissuecode[cidx_normal] #同一検体からの正常組織と腫瘍組織のペアをつくる #まずtumorベクトルとnormalベクトルで共通している要素 common_code <- intersect(tumor,normal) #####exercise###### #共通したIDがtumorベクトルのどの要素に対応するかを計算する tumor_idx <- match(,) #共通したIDがnormalベクトルのどの要素に対応するかを計算する normal_idx <- match(,) ################## #腫瘍・正常が揃っている組織IDに対応した添え字のみを取り出す cidx_tumor <- cidx_tumor[tumor_idx] cidx_normal <- cidx_normal[normal_idx] #ケース(腫瘍組織)とコントロール(対応付きの正常組織)の変数をつくる case <- expr[,cidx_tumor] ctr <- expr[,cidx_normal] #caseとcontrol群でゼロカウントの発現が極端に多いサンプルは除いておく #caseとcontrolそれぞれで、遺伝子ごとにゼロカウントの発現データが何サンプルあるかを数えよう #####exercise###### #caseの各行ごとにゼロの数を数える nzero_case <- apply(case,,) #ctrの各行ごとにゼロの数を数える nzero_ctr <- apply(ctr,,) ################## #発現量ゼロが一つも含まれていない遺伝子のみを取り出す #####exercise###### #nzero_caseベクトルのうちゼロに等しいベクトルの添え字を取り出す nonzero_case <- which() #nzero_ctrベクトルのうちゼロに等しいベクトルの添え字を取り出す nonzero_ctr <- which() ################## #ケースとコントロールで共通した行の添え字のみを取り出す common_ridx <- intersect(nonzero_case,nonzero_ctr) #common_idxを用いて、発現量ゼロが含まれている行をフィルタリング case <- case[common_ridx,] ctr <- ctr[common_ridx,] #caseとctrを列方向につなげて一つのデータにする expr <- cbind(,) #重複遺伝子の発現プロファイルを一つにまとめる #重複している遺伝子名に対しては #同じ遺伝子名のプロファイルのうち分散が最大である遺伝子プロファイルを選択する #ユニークな遺伝子名を取り出す #####exercise###### unique_gene <- unique() ################## #重複なしの遺伝子発現データを入れる変数をつくる 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)){ # i番目の遺伝子名に対応するexprの行添え字を取り出す idx <- which(unique_gene[i] == rownames(expr)) #####exercise###### #(重複)遺伝子ごとに分散を計算 # 重複していない場合はそのものの分散 # expr[idx,,drop = F]のdrop = Fは重複していないときに # expr[idx, ]がベクトルにならないように指定している(applyは行列にしか適用できないため) score <- apply(expr[idx,,drop = F], , ) ################## # 分散が最大の遺伝子名に対応する添え字を取り出す select_idx <- idx[which.max(score)] unique_expr[i,] <- expr[select_idx,] cat(i,"\n") } #exprを上書き expr <- unique_expr #caseとctrに対応する列を取り出す case <- expr[,match(colnames(case),colnames(expr))] ctr <- expr[,match(colnames(ctr),colnames(expr))] ######発現差異解析##### #1. log fold change #2. t検定 #1.log fold change #遺伝子ごとのlog fold changeを計算 #####exercise###### #遺伝子ごとに平均を計算しよう #caseの各行に対して平均を計算 case_mean <- #ctrの各行に対して平均を計算 ctr_mean <- # log fold change (LFC) を計算 log2FC <- ################## #LFC上位100遺伝子は何かを調べてみよう #####exercise###### ################## #2. t検定 #対数変換 #####exercise###### logcase <- logctr <- ################## #ケース群特異的に発現上昇あるいは減少している遺伝子群を同定してみよう。 #####exercise###### for(i in 1:n){ #各遺伝子ごとにt検定を実行する cat(i , "\n") } ################## #多重検定補正 #####exercise###### #t検定で計算したp値に対して多重検定補正してみよう ################## # -log10(FDR)とlogFCのプロット(volcano plot) pdf("data/ttest_result.pdf",10,10) plot(log2FC,-log10(qvals_twoside),pch = 20,cex=0.15,ylab="-log10(FDR)",xlab="log fold change") grid() abline(v = -1.5, col = "red") abline(v = 1.5, col = "red") abline(h = 2, col = "blue") #FDR <= 0.01 #|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 <- downgenes <- #まとめて結果をファイルに出力 #出力用の結果(発現上昇・減少の両方の結果ファイルを用意する) #####exercise###### write.table() ################## #(参考)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)) #遺伝子セット解析 #遺伝子セットの読み込み 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]] #####exercise###### #二日目で学んだ2×2の分割表を作成してみよう A <- B <- C <- D <- ################## ctb <- matrix(c(A,B,C,D),nrow=2) #####exercise###### #ctbに対してフィッシャの正確度検定を実行しよう FET <- ################## 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]] #####exercise###### #二日目で学んだ2×2の分割表を作成してみよう A <- B <- C <- D <- ################## 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) #####exercise###### #t統計量と遺伝子セットファイルを用意して #fgsea関数を実行してください fgsea_result <- ################## #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_idx <- which(fsea_result$padj <= 0.05) signif_result <- fgsea_result[signif_idx,] #####exercise###### #GO_DNA_REPAIRについてplotEnrichmentを実行してみよう p <- plotEnrichment() ################## ggsave(p,filename = "data/GSEA_ESplot.pdf") #上位FDRが上位10個についてplotGweaTable() 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,]) #####exercise###### p <- plotGseaTable() ################## ggsave(p,filename="data/GSEA_table.pdf",width = 20,height = 20)