Wednesday, February 1, 2017

Subsetting and ordering gene lists in R

In a differential expression (DE) experiment you obtained a gene list with several output columns:

> head(res)
          id baseMean baseMeanA baseMeanB foldChange log2FoldChange      pval
1 SEC24B-AS1 3.837647 3.0343409  4.640954   1.529477      0.6130381 0.4469177
2       A1BG 1.769987 0.7314147  2.808559   3.839900      1.9410687 0.3359994
3       A1CF 0.000000 0.0000000  0.000000        NaN            NaN        NA
4      GGACT 5.722972 3.5491983  7.896745   2.224938      1.1537650 0.1943902
5        A2M 0.000000 0.0000000  0.000000        NaN            NaN        NA
6      A2ML1 0.000000 0.0000000  0.000000        NaN            NaN        NA
  padj
1    1
2    1
3   NA
4    1
5   NA
6   NA
You have an ordered list of genes, for example, a list of genes in order of appearance on chromosome 7.

> cisgenes<-read.table("cis-genes")
> cisgenes
             V1
1        OSBPL3
2          CYCS
3       C7orf31
4          NPVF
5      RNU6-16P
6       MIR148A
7        NFE2L3
8     HNRNPA2B1
9          CBX3
10        SNX10
11    LOC441204
12     KIAA0087
13      C7orf71
14        SKAP2
15        HOXA1
16        HOXA2
17     HOTAIRM1
18        HOXA3
19        HOXA4
20        HOXA5
21        HOXA6
22     HOXA-AS3
23        HOXA7
24 HOXA10-HOXA9
25        HOXA9
26     HOXA-AS4
27      MIR196B
28       HOXA10
29       HOXA11
30    HOXA11-AS
31       HOXA13
32       HOTTIP
33         EVX1
34       HIBADH
35      TAX1BP1
36        JAZF1
You need to subset DE output with this gene list while preserving the order of genes.

First, using logical vector denoting all positions in the DE output that contain genes from the cisgenes list, subset the DE output:

> cisgenestable<-res[res$id %in% cisgenes$V1,]
> cisgenestable
             id     baseMean    baseMeanA   baseMeanB foldChange log2FoldChange
2320    C7orf31 2.192879e+01   18.8046996   25.052883  1.3322671     0.41388339
2343    C7orf71 0.000000e+00    0.0000000    0.000000        NaN            NaN
2655       CBX3 1.759203e+03 1825.7736950 1692.632838  0.9270770    -0.10923889
4116       CYCS 1.296481e+03 1029.4350096 1563.527666  1.5188212     0.60295200
5312       EVX1 0.000000e+00    0.0000000    0.000000        NaN            NaN
7270     HIBADH 1.627873e+03 1682.7986406 1572.946466  0.9347205    -0.09739299
7434  HNRNPA2B1 5.182713e+03 5138.0376662 5227.387459  1.0173899     0.02487263
7464      HOXA1 5.655727e+01   59.2621971   53.852338  0.9087132    -0.13810314
7465     HOXA10 5.852647e+01   53.5037521   63.549197  1.1877522     0.24823384
7466     HOXA11 3.429798e+01   31.0907859   37.505171  1.2063115     0.27060245
7467     HOXA13 0.000000e+00    0.0000000    0.000000        NaN            NaN
7468      HOXA2 3.410123e+00    2.6565888    4.163658  1.5672947     0.64827651
7469      HOXA3 3.221184e+01   29.5762769   34.847408  1.1782216     0.23661089
7470      HOXA4 3.436599e+01   32.5189782   36.213007  1.1135961     0.15522606
7471      HOXA5 4.144621e+02  382.1294264  446.794798  1.1692237     0.22555103
7472      HOXA6 4.799613e-01    0.9599225    0.000000  0.0000000           -Inf
7473      HOXA7 6.213861e+01   66.3285481   57.948674  0.8736611    -0.19485431
7474      HOXA9 6.145887e+01   59.5366534   63.381083  1.0645725     0.09027418
8140      JAZF1 4.351783e+02  414.6587972  455.697835  1.0989706     0.13615282
8350   KIAA0087 0.000000e+00    0.0000000    0.000000        NaN            NaN
11452    NFE2L3 1.338157e+01   15.3937337   11.369414  0.7385742    -0.43718525
11682      NPVF 0.000000e+00    0.0000000    0.000000        NaN            NaN
12362    OSBPL3 2.267887e+03 2432.8657172 2102.908879  0.8643752    -0.21027035
15469     SKAP2 5.550645e+02  600.1091497  510.019934  0.8498786    -0.23467129
16074     SNX10 0.000000e+00    0.0000000    0.000000        NaN            NaN
16809   TAX1BP1 1.335639e+02  132.2648988  134.862901  1.0196424     0.02806329
19566  HOTAIRM1 3.643398e+01   35.9232515   36.944708  1.0284344     0.04044978
19583  HOXA-AS3 1.585353e+00    2.0919464    1.078760  0.5156730    -0.95547163
19993    HOTTIP 9.142684e-02    0.1828537    0.000000  0.0000000           -Inf
20308 LOC441204 0.000000e+00    0.0000000    0.000000        NaN            NaN
              pval       padj
2320  0.1849568634 1.00000000
2343            NA         NA
2655  0.6771715212 1.00000000
4116  0.0003720528 0.04280161
5312            NA         NA
7270  0.7867184492 1.00000000
7434  0.8618323432 1.00000000
7464  0.6034579494 1.00000000
7465  0.5036170761 1.00000000
7466  0.5783529928 1.00000000
7467            NA         NA
7468  0.5237115868 1.00000000
7469  0.7841492733 1.00000000
7470  0.7264099930 1.00000000
7471  0.1572548872 1.00000000
7472  0.4054614692 1.00000000
7473  0.3100149829 1.00000000
7474  0.7576922017 1.00000000
8140  0.3905309791 1.00000000
8350            NA         NA
11452 0.3423512037 1.00000000
11682           NA         NA
12362 0.3125796571 1.00000000
15469 0.2314128577 1.00000000
16074           NA         NA
16809 0.9664970152 1.00000000
19566 0.9024484685 1.00000000
19583 0.6002077334 1.00000000
19993 0.9516100747 1.00000000
20308           NA         NA
Next, use match to create a position vector, for example, showing here that on third position of the cisgenes there is first entry in cisgenestable$id:

match(cisgenes$V1, cisgenestable$id)
 [1] 23  4  1 22 NA NA 21  7  3 25 30 20  2 24  8 12 27 13 14 15 16 28 17 NA 18
[26] NA NA  9 10 NA 11 29  5  6 26 19
Use this vector to re-order cisgenestable:

> cisgenestable2<-cisgenestable[match(cisgenes$V1, cisgenestable$id),]
> cisgenestable2
             id     baseMean    baseMeanA   baseMeanB foldChange log2FoldChange
12362    OSBPL3 2.267887e+03 2432.8657172 2102.908879  0.8643752    -0.21027035
4116       CYCS 1.296481e+03 1029.4350096 1563.527666  1.5188212     0.60295200
2320    C7orf31 2.192879e+01   18.8046996   25.052883  1.3322671     0.41388339
11682      NPVF 0.000000e+00    0.0000000    0.000000        NaN            NaN
NA         <NA>           NA           NA          NA         NA             NA
NA.1       <NA>           NA           NA          NA         NA             NA
11452    NFE2L3 1.338157e+01   15.3937337   11.369414  0.7385742    -0.43718525
7434  HNRNPA2B1 5.182713e+03 5138.0376662 5227.387459  1.0173899     0.02487263
2655       CBX3 1.759203e+03 1825.7736950 1692.632838  0.9270770    -0.10923889
16074     SNX10 0.000000e+00    0.0000000    0.000000        NaN            NaN
20308 LOC441204 0.000000e+00    0.0000000    0.000000        NaN            NaN
8350   KIAA0087 0.000000e+00    0.0000000    0.000000        NaN            NaN
2343    C7orf71 0.000000e+00    0.0000000    0.000000        NaN            NaN
15469     SKAP2 5.550645e+02  600.1091497  510.019934  0.8498786    -0.23467129
7464      HOXA1 5.655727e+01   59.2621971   53.852338  0.9087132    -0.13810314
7468      HOXA2 3.410123e+00    2.6565888    4.163658  1.5672947     0.64827651
19566  HOTAIRM1 3.643398e+01   35.9232515   36.944708  1.0284344     0.04044978
7469      HOXA3 3.221184e+01   29.5762769   34.847408  1.1782216     0.23661089
7470      HOXA4 3.436599e+01   32.5189782   36.213007  1.1135961     0.15522606
7471      HOXA5 4.144621e+02  382.1294264  446.794798  1.1692237     0.22555103
7472      HOXA6 4.799613e-01    0.9599225    0.000000  0.0000000           -Inf
19583  HOXA-AS3 1.585353e+00    2.0919464    1.078760  0.5156730    -0.95547163
7473      HOXA7 6.213861e+01   66.3285481   57.948674  0.8736611    -0.19485431
NA.2       <NA>           NA           NA          NA         NA             NA
7474      HOXA9 6.145887e+01   59.5366534   63.381083  1.0645725     0.09027418
NA.3       <NA>           NA           NA          NA         NA             NA
NA.4       <NA>           NA           NA          NA         NA             NA
7465     HOXA10 5.852647e+01   53.5037521   63.549197  1.1877522     0.24823384
7466     HOXA11 3.429798e+01   31.0907859   37.505171  1.2063115     0.27060245
NA.5       <NA>           NA           NA          NA         NA             NA
7467     HOXA13 0.000000e+00    0.0000000    0.000000        NaN            NaN
19993    HOTTIP 9.142684e-02    0.1828537    0.000000  0.0000000           -Inf
5312       EVX1 0.000000e+00    0.0000000    0.000000        NaN            NaN
7270     HIBADH 1.627873e+03 1682.7986406 1572.946466  0.9347205    -0.09739299
16809   TAX1BP1 1.335639e+02  132.2648988  134.862901  1.0196424     0.02806329
8140      JAZF1 4.351783e+02  414.6587972  455.697835  1.0989706     0.13615282
              pval       padj
12362 0.3125796571 1.00000000
4116  0.0003720528 0.04280161
2320  0.1849568634 1.00000000
11682           NA         NA
NA              NA         NA
NA.1            NA         NA
11452 0.3423512037 1.00000000
7434  0.8618323432 1.00000000
2655  0.6771715212 1.00000000
16074           NA         NA
20308           NA         NA
8350            NA         NA
2343            NA         NA
15469 0.2314128577 1.00000000
7464  0.6034579494 1.00000000
7468  0.5237115868 1.00000000
19566 0.9024484685 1.00000000
7469  0.7841492733 1.00000000
7470  0.7264099930 1.00000000
7471  0.1572548872 1.00000000
7472  0.4054614692 1.00000000
19583 0.6002077334 1.00000000
7473  0.3100149829 1.00000000
NA.2            NA         NA
7474  0.7576922017 1.00000000
NA.3            NA         NA
NA.4            NA         NA
7465  0.5036170761 1.00000000
7466  0.5783529928 1.00000000
NA.5            NA         NA
7467            NA         NA
19993 0.9516100747 1.00000000
5312            NA         NA
7270  0.7867184492 1.00000000
16809 0.9664970152 1.00000000
8140  0.3905309791 1.00000000
Convert Inf entries to NA and plot a barplot:

is.na(cisgenestable2)<-sapply(cisgenestable2, is.infinite)
pdf("Rplots.pdf")
barplot(cisgenestable2$log2FoldChange,names.arg=cisgenestable2$id, ylab="-logFC", border="red", las=2, cex.names=0.85)
dev.off()

No comments:

Post a Comment