2019-08-14 GEO Meta

rm(list = ls())

#install.packages("ggplot2", repos = 'https://mirror.lzu.edu.cn/CRAN/')

library(limma)

library(ggplot2)

projectPath = "E:\\201906GBM"

data4analysisPath = paste(projectPath, "DataForAnalysis", sep = "/")

deaPath = paste(projectPath, "DEA", sep = "/")

if(!dir.exists(deaPath)) dir.create(deaPath)  # 新建一个DEA文件夹，存放差异表达分析结果

datasets = c("GSE18842", "GSE19804", "GSE43458", "GSE62113")

for (dataset in datasets){

    setwd(data4analysisPath)

    dataset = "GSE15824"

    #GenesExp = read.table(paste(dataset, "_GenesExpData.txt", sep = ""), header = T, row.names = 1, sep = "\t", quote = "")




    # 对每组数据，提取样本分组信息

    if(dataset == "GSE108474"){  # 28个Normal, 221个Tumor

      SampleInfo = read.table(paste(dataset, "_SampleInfo.txt", sep = ""), header = T, row.names = 1, sep = "\t", quote = "")

      SampleInfo = SampleInfo[-c(543:550),]

      group_list = as.character( SampleInfo[, 11] )

      load("GSE108474_exprSet.Rdata")

      exprSet <- GSE108474_exprSet[,-c(543:550)]

      n_expr = exprSet[ , grep( "normal",        group_list )]

      g_expr = exprSet[ , grep( "glioblastoma",      group_list )]

      exprSet = cbind( n_expr, g_expr )

      group_list = c(rep( 'Normal', ncol( n_expr ) ),

                    rep( 'Tumor',    ncol( g_expr ) ) )

      table(group_list)

    }

    if(dataset == "GSE4290"){  # 23个Normal, 77个Tumor

      SampleInfo = read.table(paste(dataset, "_SampleInfo.txt", sep = ""), header = T, row.names = 1, sep = "\t", quote = "")

      group_list = as.character( SampleInfo[, 35] )

      load("exprSet_GSE4290.Rdata")

      exprSet <- exprSet_GSE4290

      n_expr = exprSet[ , grep( "non-tumor",        group_list )]

      g_expr = exprSet[ , grep( "glioblastoma",      group_list )]

      exprSet = cbind( n_expr, g_expr )

      group_list = c(rep( 'Normal', ncol( n_expr ) ),

                    rep( 'Tumor',    ncol( g_expr ) ) )

      table(group_list)

    }

    if(dataset == "GSE15824"){  # 5个Normal, 12个Tumor

      SampleInfo = read.table(paste(dataset, "_SampleInfo.txt", sep = ""), header = T, row.names = 1, sep = "\t", quote = "")

      group_list = as.character( SampleInfo[, 1] )

      load("exprSet_GSE15824.Rdata")

      exprSet <- exprSet_GSE15824

      ex <- exprSet

      qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))

      LogC <- (qx[5] > 100) ||

        (qx[6]-qx[1] > 50 && qx[2] > 0) ||

        (qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2)


      if (LogC) { ex[which(ex <= 0)] <- NaN

      exprSet <- log2(ex)

      print("log2 transform finished")}else{print("log2 transform not needed")}



      n_expr = exprSet[ , grep( "normal",        group_list )]

      g_expr = exprSet[ , grep( "Glioblastoma",      group_list )]

      exprSet = cbind( n_expr, g_expr )

      group_list = c(rep( 'Normal', ncol( n_expr ) ),

                    rep( 'Tumor',    ncol( g_expr ) ) )

      table(group_list)

    }

    if(dataset == "GSE50161"){  # 13个Normal, 34个Tumor

      SampleInfo = read.csv("GSE50161_SampleInfo.csv")


      group_list = as.character( SampleInfo[, 9] )

      load("exprSet_GSE50161.Rdata")

      exprSet <- exprSet_GSE50161

      n_expr = exprSet[ , grep( "normal",        group_list )]

      g_expr = exprSet[ , grep( "glioblastoma",      group_list )]

      exprSet = cbind( n_expr, g_expr )

      group_list = c(rep( 'Normal', ncol( n_expr ) ),

                    rep( 'Tumor',    ncol( g_expr ) ) )

      table(group_list)

    }

    if(dataset == "GSE116520"){  # 8个Normal, 17个Tumor

      SampleInfo = read.table(paste(dataset, "_SampleInfo.txt", sep = ""), header = T, row.names = 1, sep = "\t", quote = "")

      group_list = as.character( SampleInfo[, 1] )

      load("GSE116520_exprSet.Rdata")

      exprSet <- GSE116520_exprSet

      n_expr = exprSet[ , grep( "Control",        group_list )]

      g_expr = exprSet[ , grep( "Tumour",      group_list )]

      exprSet = cbind( n_expr, g_expr )

      group_list = c(rep( 'Normal', ncol( n_expr ) ),

                    rep( 'Tumor',    ncol( g_expr ) ) )

      table(group_list)

    }

    if(dataset == "GSE90598"){  # 7个Normal, 16个Tumor

      SampleInfo = read.table(paste(dataset, "_SampleInfo.txt", sep = ""), header = T, row.names = 1, sep = "\t", quote = "")

      group_list = as.character( SampleInfo[, 1] )

      load("GSE90598_exprSet.Rdata")

      exprSet <- GSE90598_exprSet

      n_expr = exprSet[ , grep( "Healthy",        group_list )]

      g_expr = exprSet[ , grep( "Glioblastoma",      group_list )]

      exprSet = cbind( n_expr, g_expr )

      group_list = c(rep( 'Normal', ncol( n_expr ) ),

                    rep( 'Tumor',    ncol( g_expr ) ) )

      table(group_list)

    }


    DataForDEA = exprSet

    SG = factor(group_list)

    # limma 差异表达分析

    if(dataset == "GSE19804"){      # 配对的

        SI = factor(c(1:9,1:9))  # SampleInfo$title

        design = model.matrix(~SI+SG)

        fit = lmFit(DataForDEA, design)

        fit2 = eBayes(fit)

        result = topTable(fit2, coef = "SGTumor", number = nrow(DataForDEA))    # 默认number只输出10个

    }else{  # 非配对的

        design = model.matrix(~0+SG)

        colnames(design) = levels(SG)

        contrast.matrix = makeContrasts(contrasts = "Tumor-Normal", levels = design)

        fit = lmFit(DataForDEA, design)

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

        fit2 = eBayes(fit2)

        result = topTable(fit2, coef = 1, number = nrow(DataForDEA))    # 默认number只输出10个

    }

    sig_result = subset(result, abs(logFC) >= 1 & adj.P.Val < 0.05)

    # 输出结果

    setwd(deaPath)

    write.csv(result[,c(1,2,4,5)], paste(dataset, "_all_result.csv", sep = ""))

    write.csv(sig_result[,c(1,2,4,5)], paste(dataset, "_sig_result.csv", sep = ""))


    #### vacano plot

    va_info = result[, c("logFC", "adj.P.Val")]

    va_info$Significant = 0

    for(i in 1:nrow(va_info)){

        va_info[is.na(va_info)] = 0

        if(va_info[i, 1] >= 1 & va_info[i, 2] < 0.05){

            va_info[i, 3] = "Up"

        }else if(va_info[i, 1] <= (-1) & va_info[i, 2] < 0.05){

            va_info[i, 3] = "Down"

        }else{

            va_info[i, 3] = "Normal"

        }

    }

    max_y = max(abs(c(max(va_info$logFC), min(va_info$logFC))))+1

    va_info$Significant = factor(va_info$Significant, levels = c("Up", "Normal", "Down"))

    # set color  http://www.cookbook-r.com/Graphs/Colors_(ggplot2)/

    cols = c("Up" = "#ff0000", "Normal" = "gray30", "Down" = "#32cd32")

    pdf(paste(dataset, "_vacanoplot.pdf", sep = ""), width = 6, height = 5)

    p = ggplot(data = va_info, aes(x = -log10(adj.P.Val), y = logFC, colour = Significant)) + ylim(-max_y, max_y) + scale_colour_manual(values = cols) + geom_point() + geom_hline(yintercept = 0, linetype = "dashed", color = "gray20", size = 1) + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + theme(legend.box.margin = margin(0,0,0,-10)) # + theme(legend.position="none")    # 删除所有legend

    print(p)

    dev.off()

}

最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 216,496评论 6赞 501
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 92,407评论 3赞 392
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 162,632评论 0赞 353
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 58,180评论 1赞 292
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 67,198评论 6赞 388
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 51,165评论 1赞 299
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 40,052评论 3赞 418
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 38,910评论 0赞 274
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 45,324评论 1赞 310
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,542评论 2赞 332
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,711评论 1赞 348
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 35,424评论 5赞 343
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 41,017评论 3赞 326
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,668评论 0赞 22
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,823评论 1赞 269
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 47,722评论 2赞 368
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 44,611评论 2赞 353