iSeq/edgeR.R at master · ChengLiLab/iSeq · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93

### NOT READY !!! ####
# call differential Expression genes using "DESeq" r package

call_DEG = function(condition1, condition2){
      library( "edgeR" )
      condition.vector = as.vector(as.matrix(condition)[1,colnames(Expression)])
      Expression.int = sapply(Expression, as.integer)
      rownames(Expression.int) = rownames(Expression)

      dgel = DGEList(counts=Expression.int, group=factor(condition.vector))
      #dgel = calcNormFactors(dgel)
      dgel = estimateCommonDisp(dgel)
      dgel = estimateTagwiseDisp(dgel)


      # No replications
      if(length(colnames(condition[1,condition[1,] == input$Condition1])) == 1 & length(colnames(condition[1,condition[1,] == input$Condition2])) == 1){
        cds = estimateDispersions( cds, method = "blind" )
      } else {
        cds = estimateDispersions( cds, method = "pooled" )
      }

      res = nbinomTest( cds, condition1, condition2 )
      #DESeq::plotMA(res)
      #resSig = res[ res$padj < 0.01 & (res$log2FoldChange > 1 | res$log2FoldChange < -1) & (res$baseMeanA > 1 | res$baseMeanB > 1), ]
      resSig = na.omit(res)
      resSig <- resSig[ order(resSig$pval), ]
      return(resSig)
    }

    output$conditionSelecter1 <- renderUI({
      selectInput("Condition1", "Choose condition1:", unique(t(condition)[,1]))
    })
    output$conditionSelecter2 <- renderUI({
      selectInput("Condition2", "Choose condition2:", unique(t(condition)[,1]), selected = unique(t(condition)[,1])[2])
    })

    #samplesInC1 = reactive({colnames(condition[1,condition[1,] == input$Condition1])})
    #samplesInC2 = reactive({colnames(condition[1,condition[1,] == input$Condition2])})

    DEseqButtonAct <- eventReactive(input$DEseqButton, {
      DESeqOut <- call_DEG(input$Condition1, input$Condition2)
      padj = filterDEseqP()
      foldChange = filterDEseqF()
      baseMean = filterDEseqB()
      DESeqOut = DESeqOut[DESeqOut$baseMean >= baseMean & (DESeqOut$foldChange >= foldChange | DESeqOut$foldChange <= 1/foldChange) & DESeqOut$padj <= padj,]
      #assign('DESeqOut',DESeqOut,envir=.GlobalEnv)
      DT::datatable(DESeqOut)
    })


    output$DESeqOutput <- DT::renderDataTable({
          DEseqButtonAct()
    })

    ## filter
    filterDEseqP = reactive({
        padj = switch(input$DESeqPadjCutoff,
          "0.01" = 0.01,
          "0.05" = 0.05,
          "0.1" = 0.1,
          "none" = 1
          )
        })
    filterDEseqF = reactive({
        foldChange = switch(input$DESeqfoldChangeCutoff,
          "less than 0.5 or greater than 2" = 2,     # 0.5;2 FC
          "less than 1/3 or greater than 3" = 3,     # 1/3;3
          "less than 0.25 or greater than 4" = 4,    # 0.25;4
          "none" = 1                                 # all
          )
        })
    filterDEseqB = reactive({
        baseMean = switch(input$DESeqbaseMeanCutoff,
          "1" = 1,
          "10" = 10,
          "20" = 20,
          "none" = 0
          )
        })

    # download data
    output$downloadDESeqOutput <- downloadHandler(
    filename = function() {
      paste(input$file1, 'DESeq.iSeq.csv', sep='')
    },
    content = function(file) {
      write.csv(DESeqOut, file, quote = F)
    }
    )