1. Load Count Data

data <- data.frame(read.table("counts.txt"))

or make example count data

GenerateData <- function(nRow) {
    expr_mean <- runif(1,10,100)
    expr_size <- runif(1,1,10)
    expr <- rnbinom(20,mu=expr_mean,size=expr_size)
    return(expr)
}

data <- data.frame(t(sapply(1:10000,GenerateData)))
colnames(data) <- c(paste0("A",1:10),paste0("B",1:10))
rownames(data) <- paste0("gene",1:10000)
head(data)

##        A1  A2  A3 A4  A5  A6  A7  A8  A9 A10  B1 B2  B3  B4 B5  B6 B7 B8  B9
## gene1  72 255  63 70  33 124 190 147  78  59 150 54  75 119 95  63 63 73 122
## gene2  80 103 173 79  48 140  95  97 114 147 110 84  89  22 29 109 88 62 103
## gene3  47  50  42 70  29  57  71  69  42  16  47 35  37  42 27  58 47 53  21
## gene4  43  23  80 67 112 100 207  78  41  61  40 34  54  27 93  72 25 82  25
## gene5 111  70 119 52 114 101  73  56  76  62 107 59 118  78 70 141 55 56  89
## gene6  66 146 144 87  29  83  93  89  67 154  70 59  68  90 38  65 93 85 121
##       B10
## gene1  96
## gene2 110
## gene3  49
## gene4 105
## gene5 105
## gene6  47

2. Designate Group Information and Make Design Matrix

Group <- c(rep("A",10),rep("B",10))
Group

##  [1] "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "B" "B" "B" "B" "B" "B" "B" "B" "B"
## [20] "B"

design <- model.matrix(~Group, data = data)
head(design)

##   (Intercept) GroupB
## 1           1      0
## 2           1      0
## 3           1      0
## 4           1      0
## 5           1      0
## 6           1      0

3. Normalize Count Data

dge <- DGEList(counts=data, group=Group)
dge <- calcNormFactors(dge)

## calcNormFactors has been renamed to normLibSizes

dge

## An object of class "DGEList"
## $counts
##        A1  A2  A3 A4  A5  A6  A7  A8  A9 A10  B1 B2  B3  B4 B5  B6 B7 B8  B9
## gene1  72 255  63 70  33 124 190 147  78  59 150 54  75 119 95  63 63 73 122
## gene2  80 103 173 79  48 140  95  97 114 147 110 84  89  22 29 109 88 62 103
## gene3  47  50  42 70  29  57  71  69  42  16  47 35  37  42 27  58 47 53  21
## gene4  43  23  80 67 112 100 207  78  41  61  40 34  54  27 93  72 25 82  25
## gene5 111  70 119 52 114 101  73  56  76  62 107 59 118  78 70 141 55 56  89
##       B10
## gene1  96
## gene2 110
## gene3  49
## gene4 105
## gene5 105
## 9995 more rows ...
## 
## $samples
##    group lib.size norm.factors
## A1     A   544973    0.9930964
## A2     A   544981    0.9984063
## A3     A   549961    1.0016685
## A4     A   548309    1.0043131
## A5     A   547009    1.0040734
## 15 more rows ...

4. Estimate Technical Variance

tV <- tecVarEstim(dge,design)
tV

## An object of class "TecVarList"
## $targets
##    group lib.size norm.factors
## A1     A 541210.7    0.9930964
## A2     A 544112.5    0.9984063
## A3     A 550878.6    1.0016685
## A4     A 550673.9    1.0043131
## A5     A 549237.2    1.0040734
## 15 more rows ...
## 
## $design
##   (Intercept) GroupB
## 1           1      0
## 2           1      0
## 3           1      0
## 4           1      0
## 5           1      0
## 15 more rows ...
## 
## $logcpm
##             A1       A2       A3       A4       A5       A6       A7       A8
## gene1 7.094603 8.883530 6.882674 7.030741 5.993754 7.826846 8.437269 8.081949
## gene2 7.242758 7.591996 8.311440 7.200655 6.508332 7.999286 7.452342 7.492147
## gene3 6.499564 6.577923 6.319797 7.030741 5.818665 6.732478 7.042335 7.012549
## gene4 6.376630 5.520767 7.217833 6.969344 7.697373 7.522047 8.559655 7.184725
## gene5 7.705528 7.047555 7.779113 6.615735 7.722464 7.536119 7.081321 6.720792
##             A9      A10       B1       B2       B3       B4       B5       B6
## gene1 7.173945 6.807819 8.123123 6.683449 7.127431 7.759248 7.468730 6.885365
## gene2 7.709901 8.096507 7.682498 7.302512 7.368425 5.426976 5.822806 7.657397
## gene3 6.311705 5.045966 6.489598 6.085322 6.146139 6.300522 5.726570 6.769891
## gene4 6.278551 6.854375 6.267130 6.045775 6.668035 5.699644 7.438671 7.072446
## gene5 7.137427 6.877102 7.643322 6.806867 7.767502 7.162519 7.038709 8.022848
##             B7       B8       B9      B10
## gene1 6.899075 7.078694 7.808891 7.491100
## gene2 7.368670 6.849943 7.568968 7.683776
## gene3 6.491288 6.631379 5.378942 6.548582
## gene4 5.631037 7.242151 5.610143 7.617878
## gene5 6.709541 6.707972 7.362546 7.617878
## 9995 more rows ...
## 
## $tecVar
##                A1          A2          A3          A4          A5          A6
## gene1 0.010114034 0.010063320 0.009947076 0.009950552 0.009975021 0.009876954
## gene2 0.009207474 0.009161213 0.009055154 0.009058326 0.009080654 0.008991168
## gene3 0.019638222 0.019536422 0.019303137 0.019310112 0.019359211 0.019162476
## gene4 0.013219630 0.013153352 0.013001406 0.013005951 0.013037939 0.012909958
## gene5 0.011516535 0.011458966 0.011326973 0.011330921 0.011358709 0.011247338
##                A7          A8          A9         A10         B1         B2
## gene1 0.009895332 0.009963804 0.009892254 0.010060067 0.01063488 0.01064313
## gene2 0.009007938 0.009070419 0.009005129 0.009158246 0.01266353 0.01267334
## gene3 0.019199336 0.019336703 0.019193163 0.019529892 0.02210769 0.02212548
## gene4 0.012933887 0.013023275 0.012929880 0.013149101 0.01832997 0.01834470
## gene5 0.011268210 0.011345971 0.011264714 0.011455273 0.01098265 0.01099117
##               B3         B4         B5         B6         B7         B8
## gene1 0.01053299 0.01039287 0.01059049 0.01054952 0.01064708 0.01045736
## gene2 0.01254253 0.01237634 0.01261075 0.01256213 0.01267803 0.01245283
## gene3 0.02188820 0.02158672 0.02201200 0.02192378 0.02213400 0.02172545
## gene4 0.01814802 0.01789798 0.01825068 0.01817753 0.01835175 0.01801304
## gene5 0.01087749 0.01073289 0.01093683 0.01089454 0.01099524 0.01079944
##               B9        B10
## gene1 0.01049185 0.01064389
## gene2 0.01249373 0.01267424
## gene3 0.02179966 0.02212713
## gene4 0.01807458 0.01834606
## gene5 0.01083503 0.01099195
## 9995 more rows ...

5. Fit BALLI and See Top Significant Genes

fit <- balli(tV,intV=2)
fit

## An object of class "Balli"
## $Result
##       log2FC_GroupB       lLLI     lBALLI      pLLI    pBALLI       BCF
## gene1   -0.08879523 0.08944797 0.07950932 0.7648801 0.7779635 0.1249999
## gene2   -0.48745204 2.68228922 2.38426500 0.1014698 0.1225625 0.1249963
## gene3   -0.18252535 0.62845747 0.55863197 0.4279217 0.4548120 0.1249937
## gene4   -0.48894684 1.91830326 1.70516493 0.1660438 0.1916139 0.1249957
## gene5    0.06165777 0.10305350 0.09160314 0.7481956 0.7621489 0.1249996
## 9995 more rows ...
## 
## $topGenes
##          log2FC_GroupB         pLLI       pBALLI   adjpLLI adjpBALLI
## gene5528     1.4443636 8.479278e-05 0.0002107813 0.4744766 0.8932819
## gene4740    -1.1265520 1.242018e-04 0.0002953372 0.4744766 0.8932819
## gene8600    -0.8356578 2.038886e-04 0.0004613161 0.4744766 0.8932819
## gene6619     0.8274648 3.278305e-04 0.0007070187 0.4744766 0.8932819
## gene201      1.1107640 3.430630e-04 0.0007355785 0.4744766 0.8932819
## 9995 more rows ...