Hands-On Tidyverse - with count data

Setup

[1]:

# Load required packages
library(tidyverse)
library(foreach)
library(stringr)
library(haven)

library(DESeq2)
library(tools)
library(limma)
library(qvalue)

library(ggplot2)
library(RColorBrewer)
library(gridExtra)
library(dendextend)

library(plotly)

Registered S3 methods overwritten by 'ggplot2':
  method         from
  [.quosures     rlang
  c.quosures     rlang
  print.quosures rlang
── Attaching packages ─────────────────────────────────────── tidyverse 1.2.1 ──
✔ ggplot2 3.1.1     ✔ purrr   0.3.2
✔ tibble  2.1.2     ✔ dplyr   0.8.1
✔ tidyr   0.8.3     ✔ stringr 1.4.0
✔ readr   1.3.1     ✔ forcats 0.4.0
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

Attaching package: ‘foreach’

The following objects are masked from ‘package:purrr’:

    accumulate, when

Loading required package: S4Vectors
Loading required package: stats4
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: ‘BiocGenerics’

The following objects are masked from ‘package:parallel’:

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from ‘package:dplyr’:

    combine, intersect, setdiff, union

The following objects are masked from ‘package:stats’:

    IQR, mad, sd, var, xtabs

The following objects are masked from ‘package:base’:

    anyDuplicated, append, as.data.frame, basename, cbind, colnames,
    dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
    grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
    order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
    rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
    union, unique, unsplit, which, which.max, which.min


Attaching package: ‘S4Vectors’

The following objects are masked from ‘package:dplyr’:

    first, rename

The following object is masked from ‘package:tidyr’:

    expand

The following object is masked from ‘package:base’:

    expand.grid

Loading required package: IRanges

Attaching package: ‘IRanges’

The following objects are masked from ‘package:dplyr’:

    collapse, desc, slice

The following object is masked from ‘package:purrr’:

    reduce

Loading required package: GenomicRanges
Loading required package: GenomeInfoDb
Loading required package: SummarizedExperiment
Loading required package: Biobase
Welcome to Bioconductor

    Vignettes contain introductory material; view with
    'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.

Loading required package: DelayedArray
Loading required package: matrixStats

Attaching package: ‘matrixStats’

The following objects are masked from ‘package:Biobase’:

    anyMissing, rowMedians

The following object is masked from ‘package:dplyr’:

    count

Loading required package: BiocParallel

Attaching package: ‘DelayedArray’

The following objects are masked from ‘package:matrixStats’:

    colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges

The following object is masked from ‘package:purrr’:

    simplify

The following objects are masked from ‘package:base’:

    aperm, apply, rowsum


Attaching package: ‘limma’

The following object is masked from ‘package:DESeq2’:

    plotMA

The following object is masked from ‘package:BiocGenerics’:

    plotMA


Attaching package: ‘gridExtra’

The following object is masked from ‘package:Biobase’:

    combine

The following object is masked from ‘package:BiocGenerics’:

    combine

The following object is masked from ‘package:dplyr’:

    combine


---------------------
Welcome to dendextend version 1.12.0
Type citation('dendextend') for how to cite the package.

Type browseVignettes(package = 'dendextend') for the package vignette.
The github page is: https://github.com/talgalili/dendextend/

Suggestions and bug-reports can be submitted at: https://github.com/talgalili/dendextend/issues
Or contact: <tal.galili@gmail.com>

        To suppress this message use:  suppressPackageStartupMessages(library(dendextend))
---------------------


Attaching package: ‘dendextend’

The following object is masked from ‘package:stats’:

    cutree


Attaching package: ‘plotly’

The following object is masked from ‘package:IRanges’:

    slice

The following object is masked from ‘package:S4Vectors’:

    rename

The following object is masked from ‘package:ggplot2’:

    last_plot

The following object is masked from ‘package:stats’:

    filter

The following object is masked from ‘package:graphics’:

    layout


[2]:

# set directories
DATDIR <- "/data/hts2018_pilot/star_counts"
CURDIR <- "/home/jovyan/work/scratch/analysis_output"
OUTDIR <- file.path(CURDIR, "out")
IMGDIR <- file.path(CURDIR, "img")

# Metadata (metadtfile)
METADTFILE <- "/home/jovyan/work/HTS2018-notebooks/josh/info/2018_pilot_metadata_anon.tsv"

Reading in count data

The gene counts from the pilot data from the 2018 course are in the directory /data/hts2018_pilot/start_counts.

Let’s have a look at them:

[3]:

list.files("/data/hts2018_pilot/star_counts/")
  1. '1_MA_J_S18_L001_ReadsPerGene.out.tab'
  2. '1_MA_J_S18_L002_ReadsPerGene.out.tab'
  3. '1_MA_J_S18_L003_ReadsPerGene.out.tab'
  4. '1_MA_J_S18_L004_ReadsPerGene.out.tab'
  5. '1_RZ_J_S26_L001_ReadsPerGene.out.tab'
  6. '1_RZ_J_S26_L002_ReadsPerGene.out.tab'
  7. '1_RZ_J_S26_L003_ReadsPerGene.out.tab'
  8. '1_RZ_J_S26_L004_ReadsPerGene.out.tab'
  9. '10_MA_C_S3_L001_ReadsPerGene.out.tab'
  10. '10_MA_C_S3_L002_ReadsPerGene.out.tab'
  11. '10_MA_C_S3_L003_ReadsPerGene.out.tab'
  12. '10_MA_C_S3_L004_ReadsPerGene.out.tab'
  13. '10_RZ_C_S11_L001_ReadsPerGene.out.tab'
  14. '10_RZ_C_S11_L002_ReadsPerGene.out.tab'
  15. '10_RZ_C_S11_L003_ReadsPerGene.out.tab'
  16. '10_RZ_C_S11_L004_ReadsPerGene.out.tab'
  17. '11_MA_J_S20_L001_ReadsPerGene.out.tab'
  18. '11_MA_J_S20_L002_ReadsPerGene.out.tab'
  19. '11_MA_J_S20_L003_ReadsPerGene.out.tab'
  20. '11_MA_J_S20_L004_ReadsPerGene.out.tab'
  21. '11_RZ_J_S28_L001_ReadsPerGene.out.tab'
  22. '11_RZ_J_S28_L002_ReadsPerGene.out.tab'
  23. '11_RZ_J_S28_L003_ReadsPerGene.out.tab'
  24. '11_RZ_J_S28_L004_ReadsPerGene.out.tab'
  25. '12_MA_P_S36_L001_ReadsPerGene.out.tab'
  26. '12_MA_P_S36_L002_ReadsPerGene.out.tab'
  27. '12_MA_P_S36_L003_ReadsPerGene.out.tab'
  28. '12_MA_P_S36_L004_ReadsPerGene.out.tab'
  29. '12_RZ_P_S44_L001_ReadsPerGene.out.tab'
  30. '12_RZ_P_S44_L002_ReadsPerGene.out.tab'
  31. '12_RZ_P_S44_L003_ReadsPerGene.out.tab'
  32. '12_RZ_P_S44_L004_ReadsPerGene.out.tab'
  33. '13_MA_J_S21_L001_ReadsPerGene.out.tab'
  34. '13_MA_J_S21_L002_ReadsPerGene.out.tab'
  35. '13_MA_J_S21_L003_ReadsPerGene.out.tab'
  36. '13_MA_J_S21_L004_ReadsPerGene.out.tab'
  37. '13_RZ_J_S29_L001_ReadsPerGene.out.tab'
  38. '13_RZ_J_S29_L002_ReadsPerGene.out.tab'
  39. '13_RZ_J_S29_L003_ReadsPerGene.out.tab'
  40. '13_RZ_J_S29_L004_ReadsPerGene.out.tab'
  41. '14_MA_C_S4_L001_ReadsPerGene.out.tab'
  42. '14_MA_C_S4_L002_ReadsPerGene.out.tab'
  43. '14_MA_C_S4_L003_ReadsPerGene.out.tab'
  44. '14_MA_C_S4_L004_ReadsPerGene.out.tab'
  45. '14_RZ_C_S12_L001_ReadsPerGene.out.tab'
  46. '14_RZ_C_S12_L002_ReadsPerGene.out.tab'
  47. '14_RZ_C_S12_L003_ReadsPerGene.out.tab'
  48. '14_RZ_C_S12_L004_ReadsPerGene.out.tab'
  49. '15_MA_C_S5_L001_ReadsPerGene.out.tab'
  50. '15_MA_C_S5_L002_ReadsPerGene.out.tab'
  51. '15_MA_C_S5_L003_ReadsPerGene.out.tab'
  52. '15_MA_C_S5_L004_ReadsPerGene.out.tab'
  53. '15_RZ_C_S13_L001_ReadsPerGene.out.tab'
  54. '15_RZ_C_S13_L002_ReadsPerGene.out.tab'
  55. '15_RZ_C_S13_L003_ReadsPerGene.out.tab'
  56. '15_RZ_C_S13_L004_ReadsPerGene.out.tab'
  57. '16_MA_P_S37_L001_ReadsPerGene.out.tab'
  58. '16_MA_P_S37_L002_ReadsPerGene.out.tab'
  59. '16_MA_P_S37_L003_ReadsPerGene.out.tab'
  60. '16_MA_P_S37_L004_ReadsPerGene.out.tab'
  61. '16_RZ_P_S45_L001_ReadsPerGene.out.tab'
  62. '16_RZ_P_S45_L002_ReadsPerGene.out.tab'
  63. '16_RZ_P_S45_L003_ReadsPerGene.out.tab'
  64. '16_RZ_P_S45_L004_ReadsPerGene.out.tab'
  65. '2_MA_C_S1_L001_ReadsPerGene.out.tab'
  66. '2_MA_C_S1_L002_ReadsPerGene.out.tab'
  67. '2_MA_C_S1_L003_ReadsPerGene.out.tab'
  68. '2_MA_C_S1_L004_ReadsPerGene.out.tab'
  69. '2_RZ_C_S9_L001_ReadsPerGene.out.tab'
  70. '2_RZ_C_S9_L002_ReadsPerGene.out.tab'
  71. '2_RZ_C_S9_L003_ReadsPerGene.out.tab'
  72. '2_RZ_C_S9_L004_ReadsPerGene.out.tab'
  73. '2_TOT_C_S17_L001_ReadsPerGene.out.tab'
  74. '2_TOT_C_S17_L002_ReadsPerGene.out.tab'
  75. '2_TOT_C_S17_L003_ReadsPerGene.out.tab'
  76. '2_TOT_C_S17_L004_ReadsPerGene.out.tab'
  77. '21_MA_C_S6_L001_ReadsPerGene.out.tab'
  78. '21_MA_C_S6_L002_ReadsPerGene.out.tab'
  79. '21_MA_C_S6_L003_ReadsPerGene.out.tab'
  80. '21_MA_C_S6_L004_ReadsPerGene.out.tab'
  81. '21_RZ_C_S14_L001_ReadsPerGene.out.tab'
  82. '21_RZ_C_S14_L002_ReadsPerGene.out.tab'
  83. '21_RZ_C_S14_L003_ReadsPerGene.out.tab'
  84. '21_RZ_C_S14_L004_ReadsPerGene.out.tab'
  85. '22_MA_C_S7_L001_ReadsPerGene.out.tab'
  86. '22_MA_C_S7_L002_ReadsPerGene.out.tab'
  87. '22_MA_C_S7_L003_ReadsPerGene.out.tab'
  88. '22_MA_C_S7_L004_ReadsPerGene.out.tab'
  89. '22_RZ_C_S15_L001_ReadsPerGene.out.tab'
  90. '22_RZ_C_S15_L002_ReadsPerGene.out.tab'
  91. '22_RZ_C_S15_L003_ReadsPerGene.out.tab'
  92. '22_RZ_C_S15_L004_ReadsPerGene.out.tab'
  93. '23_MA_J_S22_L001_ReadsPerGene.out.tab'
  94. '23_MA_J_S22_L002_ReadsPerGene.out.tab'
  95. '23_MA_J_S22_L003_ReadsPerGene.out.tab'
  96. '23_MA_J_S22_L004_ReadsPerGene.out.tab'
  97. '23_RZ_J_S30_L001_ReadsPerGene.out.tab'
  98. '23_RZ_J_S30_L002_ReadsPerGene.out.tab'
  99. '23_RZ_J_S30_L003_ReadsPerGene.out.tab'
  100. '23_RZ_J_S30_L004_ReadsPerGene.out.tab'
  101. '24_MA_J_S23_L001_ReadsPerGene.out.tab'
  102. '24_MA_J_S23_L002_ReadsPerGene.out.tab'
  103. '24_MA_J_S23_L003_ReadsPerGene.out.tab'
  104. '24_MA_J_S23_L004_ReadsPerGene.out.tab'
  105. '24_RZ_J_S31_L001_ReadsPerGene.out.tab'
  106. '24_RZ_J_S31_L002_ReadsPerGene.out.tab'
  107. '24_RZ_J_S31_L003_ReadsPerGene.out.tab'
  108. '24_RZ_J_S31_L004_ReadsPerGene.out.tab'
  109. '26_MA_C_S8_L001_ReadsPerGene.out.tab'
  110. '26_MA_C_S8_L002_ReadsPerGene.out.tab'
  111. '26_MA_C_S8_L003_ReadsPerGene.out.tab'
  112. '26_MA_C_S8_L004_ReadsPerGene.out.tab'
  113. '26_RZ_C_S16_L001_ReadsPerGene.out.tab'
  114. '26_RZ_C_S16_L002_ReadsPerGene.out.tab'
  115. '26_RZ_C_S16_L003_ReadsPerGene.out.tab'
  116. '26_RZ_C_S16_L004_ReadsPerGene.out.tab'
  117. '27_MA_P_S38_L001_ReadsPerGene.out.tab'
  118. '27_MA_P_S38_L002_ReadsPerGene.out.tab'
  119. '27_MA_P_S38_L003_ReadsPerGene.out.tab'
  120. '27_MA_P_S38_L004_ReadsPerGene.out.tab'
  121. '27_RZ_P_S46_L001_ReadsPerGene.out.tab'
  122. '27_RZ_P_S46_L002_ReadsPerGene.out.tab'
  123. '27_RZ_P_S46_L003_ReadsPerGene.out.tab'
  124. '27_RZ_P_S46_L004_ReadsPerGene.out.tab'
  125. '3_MA_J_S19_L001_ReadsPerGene.out.tab'
  126. '3_MA_J_S19_L002_ReadsPerGene.out.tab'
  127. '3_MA_J_S19_L003_ReadsPerGene.out.tab'
  128. '3_MA_J_S19_L004_ReadsPerGene.out.tab'
  129. '3_RZ_J_S27_L001_ReadsPerGene.out.tab'
  130. '3_RZ_J_S27_L002_ReadsPerGene.out.tab'
  131. '3_RZ_J_S27_L003_ReadsPerGene.out.tab'
  132. '3_RZ_J_S27_L004_ReadsPerGene.out.tab'
  133. '3_TOT_J_S34_L001_ReadsPerGene.out.tab'
  134. '3_TOT_J_S34_L002_ReadsPerGene.out.tab'
  135. '3_TOT_J_S34_L003_ReadsPerGene.out.tab'
  136. '3_TOT_J_S34_L004_ReadsPerGene.out.tab'
  137. '35_MA_P_S39_L001_ReadsPerGene.out.tab'
  138. '35_MA_P_S39_L002_ReadsPerGene.out.tab'
  139. '35_MA_P_S39_L003_ReadsPerGene.out.tab'
  140. '35_MA_P_S39_L004_ReadsPerGene.out.tab'
  141. '35_RZ_P_S47_L001_ReadsPerGene.out.tab'
  142. '35_RZ_P_S47_L002_ReadsPerGene.out.tab'
  143. '35_RZ_P_S47_L003_ReadsPerGene.out.tab'
  144. '35_RZ_P_S47_L004_ReadsPerGene.out.tab'
  145. '36_MA_J_S24_L001_ReadsPerGene.out.tab'
  146. '36_MA_J_S24_L002_ReadsPerGene.out.tab'
  147. '36_MA_J_S24_L003_ReadsPerGene.out.tab'
  148. '36_MA_J_S24_L004_ReadsPerGene.out.tab'
  149. '36_RZ_J_S32_L001_ReadsPerGene.out.tab'
  150. '36_RZ_J_S32_L002_ReadsPerGene.out.tab'
  151. '36_RZ_J_S32_L003_ReadsPerGene.out.tab'
  152. '36_RZ_J_S32_L004_ReadsPerGene.out.tab'
  153. '38_MA_P_S40_L001_ReadsPerGene.out.tab'
  154. '38_MA_P_S40_L002_ReadsPerGene.out.tab'
  155. '38_MA_P_S40_L003_ReadsPerGene.out.tab'
  156. '38_MA_P_S40_L004_ReadsPerGene.out.tab'
  157. '38_RZ_P_S48_L001_ReadsPerGene.out.tab'
  158. '38_RZ_P_S48_L002_ReadsPerGene.out.tab'
  159. '38_RZ_P_S48_L003_ReadsPerGene.out.tab'
  160. '38_RZ_P_S48_L004_ReadsPerGene.out.tab'
  161. '4_MA_P_S35_L001_ReadsPerGene.out.tab'
  162. '4_MA_P_S35_L002_ReadsPerGene.out.tab'
  163. '4_MA_P_S35_L003_ReadsPerGene.out.tab'
  164. '4_MA_P_S35_L004_ReadsPerGene.out.tab'
  165. '4_RZ_P_S43_L001_ReadsPerGene.out.tab'
  166. '4_RZ_P_S43_L002_ReadsPerGene.out.tab'
  167. '4_RZ_P_S43_L003_ReadsPerGene.out.tab'
  168. '4_RZ_P_S43_L004_ReadsPerGene.out.tab'
  169. '4_TOT_P_S51_L001_ReadsPerGene.out.tab'
  170. '4_TOT_P_S51_L002_ReadsPerGene.out.tab'
  171. '4_TOT_P_S51_L003_ReadsPerGene.out.tab'
  172. '4_TOT_P_S51_L004_ReadsPerGene.out.tab'
  173. '40_MA_J_S25_L001_ReadsPerGene.out.tab'
  174. '40_MA_J_S25_L002_ReadsPerGene.out.tab'
  175. '40_MA_J_S25_L003_ReadsPerGene.out.tab'
  176. '40_MA_J_S25_L004_ReadsPerGene.out.tab'
  177. '40_RZ_J_S33_L001_ReadsPerGene.out.tab'
  178. '40_RZ_J_S33_L002_ReadsPerGene.out.tab'
  179. '40_RZ_J_S33_L003_ReadsPerGene.out.tab'
  180. '40_RZ_J_S33_L004_ReadsPerGene.out.tab'
  181. '45_MA_P_S41_L001_ReadsPerGene.out.tab'
  182. '45_MA_P_S41_L002_ReadsPerGene.out.tab'
  183. '45_MA_P_S41_L003_ReadsPerGene.out.tab'
  184. '45_MA_P_S41_L004_ReadsPerGene.out.tab'
  185. '45_RZ_P_S49_L001_ReadsPerGene.out.tab'
  186. '45_RZ_P_S49_L002_ReadsPerGene.out.tab'
  187. '45_RZ_P_S49_L003_ReadsPerGene.out.tab'
  188. '45_RZ_P_S49_L004_ReadsPerGene.out.tab'
  189. '47_MA_P_S42_L001_ReadsPerGene.out.tab'
  190. '47_MA_P_S42_L002_ReadsPerGene.out.tab'
  191. '47_MA_P_S42_L003_ReadsPerGene.out.tab'
  192. '47_MA_P_S42_L004_ReadsPerGene.out.tab'
  193. '47_RZ_P_S50_L001_ReadsPerGene.out.tab'
  194. '47_RZ_P_S50_L002_ReadsPerGene.out.tab'
  195. '47_RZ_P_S50_L003_ReadsPerGene.out.tab'
  196. '47_RZ_P_S50_L004_ReadsPerGene.out.tab'
  197. '9_MA_C_S2_L001_ReadsPerGene.out.tab'
  198. '9_MA_C_S2_L002_ReadsPerGene.out.tab'
  199. '9_MA_C_S2_L003_ReadsPerGene.out.tab'
  200. '9_MA_C_S2_L004_ReadsPerGene.out.tab'
  201. '9_RZ_C_S10_L001_ReadsPerGene.out.tab'
  202. '9_RZ_C_S10_L002_ReadsPerGene.out.tab'
  203. '9_RZ_C_S10_L003_ReadsPerGene.out.tab'
  204. '9_RZ_C_S10_L004_ReadsPerGene.out.tab'

We can also see these in the terminal window (open a terminal window and use the bash command ‘ls’). Let’s quickly go to the terminal and do this. Also, we can look at the content of the files.

[ ]:



[4]:

# Save the names in variable
stardirs <- list.files(DATDIR)

# Look at the beginning of the first file from R
cmdstr <- paste("head", file.path(DATDIR, stardirs[1]))
cmdout <- system(cmdstr, intern = TRUE)
str_split(cmdout, pattern = "\t")
    1. 'N_unmapped'
    2. '2690'
    3. '2690'
    4. '2690'
    1. 'N_multimapping'
    2. '66100'
    3. '66100'
    4. '66100'
    1. 'N_noFeature'
    2. '10626'
    3. '2238382'
    4. '20347'
    1. 'N_ambiguous'
    2. '173170'
    3. '1622'
    4. '647'
    1. 'CNAG_04548'
    2. '0'
    3. '0'
    4. '0'
    1. 'CNAG_07303'
    2. '0'
    3. '0'
    4. '0'
    1. 'CNAG_07304'
    2. '8'
    3. '0'
    4. '8'
    1. 'CNAG_00001'
    2. '0'
    3. '0'
    4. '0'
    1. 'CNAG_07305'
    2. '0'
    3. '0'
    4. '0'
    1. 'CNAG_00002'
    2. '66'
    3. '0'
    4. '66'

There are several things to note: * There are four columns. We only want the first (gene name) and the fourth (count). * There are no headers. * This is a tab-delimited file (we can’t see this, but what we can see is that it is not a csv)

Exercise: 1. How many files are in the directory? 2. Print the first 10 filenames 3. Use the command read_tsv to read in the second file and save it in a tibble called “sample_2”. Use the note above to pass the correct options to read_tsv.

[5]:

sample_file <- paste0(DATDIR, "/", stardirs[2])
sample_2 <- read_tsv(sample_file, col_names = FALSE)

Parsed with column specification:
cols(
  X1 = col_character(),
  X2 = col_double(),
  X3 = col_double(),
  X4 = col_double()
)

[6]:

sample_2 %>% head
A tibble: 6 × 4
X1X2X3X4
<chr><dbl><dbl><dbl>
N_unmapped 2684 2684 2684
N_multimapping 65234 6523465234
N_noFeature 10340220418720004
N_ambiguous 169504 1523 652
CNAG_04548 0 0 0
CNAG_07303 0 0 0

Our goal is to read in all of these files and collect the first and fourth columns into one large tibble. Let’s first do this for two files.

Exercise:

1. Read in the first two files, one into a tibble called df1, the other into a  tibble called df2.
2. Remove the middle two columns using dplyr, and rename the remaining two columns 'gene' and the name of the count file.
3. Join the two tibbles together using 'gene' as the key.

[7]:

# Fancy way to read in columns 1 and 4 only

coltypes<- "?--i"

sample_file <- paste0(DATDIR, "/", stardirs[1])
df1 <- readr::read_tsv(sample_file, col_types = coltypes, col_names = c("gene", stardirs[1]))

sample_file <- paste0(DATDIR, "/", stardirs[2])
df2 <- readr::read_tsv(sample_file, col_types = coltypes, col_names = c("gene", stardirs[2]))

[8]:

df1 %>% head
df2 %>% head
A tibble: 6 × 2
gene1_MA_J_S18_L001_ReadsPerGene.out.tab
<chr><int>
N_unmapped 2690
N_multimapping66100
N_noFeature 20347
N_ambiguous 647
CNAG_04548 0
CNAG_07303 0
A tibble: 6 × 2
gene1_MA_J_S18_L002_ReadsPerGene.out.tab
<chr><int>
N_unmapped 2684
N_multimapping65234
N_noFeature 20004
N_ambiguous 652
CNAG_04548 0
CNAG_07303 0
 
[9]:

full_join(df1, df2, by = "gene") %>% head
A tibble: 6 × 3
gene1_MA_J_S18_L001_ReadsPerGene.out.tab1_MA_J_S18_L002_ReadsPerGene.out.tab
<chr><int><int>
N_unmapped 2690 2684
N_multimapping6610065234
N_noFeature 2034720004
N_ambiguous 647 652
CNAG_04548 0 0
CNAG_07303 0 0

Of course, we don’t want to do this manually for every file. We’ll use the foreach package in R to iterate over the files. This will require defining some of the steps above as functions, so first let’s review what a function is.

Functions

Functions are simply objects that do something. In the functional programming paradigm, functions should be self-contained, in that they receive as inputs all the variables they need and do not modify anything else. They ‘return’ an output.

Example

Good

[10]:

myfunction_add <- function(a,b){
    a + b   # In R, the last statement is what is returned
}

[11]:

myfunction_add(1,2)
3

Bad

[12]:

a <- 1
b <- 2

myfunction_add <- function(){
    a + b # We are using values from the 'global environment' instead of passing them in

}

[13]:

myfunction_add()
3

Exercise:

1. Write a function to multiply two numbers and return the result.
2. Write a function to join two dataframes

[14]:

mycombine <- function(df1, df2) {
    # Combine two data frames by gene names
    #
    # Args:
    #   df1 (Dataframe): the first count data
    #   df2 (Dataframe): the second count data
    #
    # Returns:
    #   (Dataframe) The combined data frame of df1 and df2
    full_join(df1, df2, by = "gene")
}

myfile <- function(filedir, filename) {
    # Get the absolute paths of a file
    #
    # Args:
    #   filedir  (Character): the directory of the folder
    #   filename (Character): the filename
    #
    # Returns:
    #   (Character) the directory of the input file
    file.path(filedir, filename)
}



[15]:

coltypes<- "?--i"

out <- foreach(stardir = stardirs, .combine = mycombine) %do% {

    # get a directory of each count file
    cntfile <- myfile(DATDIR, stardir)

    # read in the count file
    readr::read_tsv(cntfile, col_names = FALSE, col_types = coltypes) %>%
           dplyr::rename_(.dots=setNames(names(.), c("gene",stardir)))
           #dplyr::rename("gene" = "X1", `stardir` = "X4")
}

Warning message:
“rename_() is deprecated.
Please use rename() instead

The 'programming' vignette or the tidyeval book can help you
to program with rename() : https://tidyeval.tidyverse.org
This warning is displayed once per session.”

[16]:

out %>% head
A tibble: 6 × 205
gene1_MA_J_S18_L001_ReadsPerGene.out.tab1_MA_J_S18_L002_ReadsPerGene.out.tab1_MA_J_S18_L003_ReadsPerGene.out.tab1_MA_J_S18_L004_ReadsPerGene.out.tab1_RZ_J_S26_L001_ReadsPerGene.out.tab1_RZ_J_S26_L002_ReadsPerGene.out.tab1_RZ_J_S26_L003_ReadsPerGene.out.tab1_RZ_J_S26_L004_ReadsPerGene.out.tab10_MA_C_S3_L001_ReadsPerGene.out.tab⋯47_RZ_P_S50_L003_ReadsPerGene.out.tab47_RZ_P_S50_L004_ReadsPerGene.out.tab9_MA_C_S2_L001_ReadsPerGene.out.tab9_MA_C_S2_L002_ReadsPerGene.out.tab9_MA_C_S2_L003_ReadsPerGene.out.tab9_MA_C_S2_L004_ReadsPerGene.out.tab9_RZ_C_S10_L001_ReadsPerGene.out.tab9_RZ_C_S10_L002_ReadsPerGene.out.tab9_RZ_C_S10_L003_ReadsPerGene.out.tab9_RZ_C_S10_L004_ReadsPerGene.out.tab
<chr><int><int><int><int><int><int><int><int><int>⋯<int><int><int><int><int><int><int><int><int><int>
N_unmapped 2690 2684 2672 2585 7218 7022 7355 707638278⋯ 10036 9962 2245 2291 2276 2105 3386 3612 3283 4853
N_multimapping6610065234665386506639584838807940133839549064124⋯ 536339 52936376258744907637075176149388149618150874156664
N_noFeature 2034720004205492050576814675565477774977371228540⋯1055322104717127956276382837228459503625504801510186525524
N_ambiguous 647 652 697 616 1431 1337 1425 1322 147⋯ 1260 1236 903 848 943 838 1354 1333 1359 1357
CNAG_04548 0 0 0 1 0 0 0 1 0⋯ 0 0 0 0 0 0 1 0 0 0
CNAG_07303 0 0 0 0 0 0 0 0 0⋯ 0 0 0 0 0 0 0 0 0 0

Exercise: Create the ‘out’ tibble using a for loop instead of foreach.

Gather and spread

Now, we have a few other things to fix. To begin with, the first four rows are really summaries and not genes.

[17]:

### Gather and spread the first four rows
out %>%
    dplyr::slice(1:4) %>%
    gather(expid, value, -gene) %>%
    spread(gene, value) %>%
    rename_(.dots = setNames(names(.), c("expid", "namb", "nmulti", "nnofeat","nunmap"))) ->
    nmisc

[18]:

nmisc %>% head
A tibble: 6 × 5
expidnambnmultinnofeatnunmap
<chr><int><int><int><int>
1_MA_J_S18_L001_ReadsPerGene.out.tab 647 66100 203472690
1_MA_J_S18_L002_ReadsPerGene.out.tab 652 65234 200042684
1_MA_J_S18_L003_ReadsPerGene.out.tab 697 66538 205492672
1_MA_J_S18_L004_ReadsPerGene.out.tab 616 65066 205052585
1_RZ_J_S26_L001_ReadsPerGene.out.tab14313958487681467218
1_RZ_J_S26_L002_ReadsPerGene.out.tab13373880797556547022

Let’s break this down and see what each step does.

[19]:

out %>%
    dplyr::slice(1:4) -> temp1

temp1
A spec_tbl_df: 4 × 205
gene1_MA_J_S18_L001_ReadsPerGene.out.tab1_MA_J_S18_L002_ReadsPerGene.out.tab1_MA_J_S18_L003_ReadsPerGene.out.tab1_MA_J_S18_L004_ReadsPerGene.out.tab1_RZ_J_S26_L001_ReadsPerGene.out.tab1_RZ_J_S26_L002_ReadsPerGene.out.tab1_RZ_J_S26_L003_ReadsPerGene.out.tab1_RZ_J_S26_L004_ReadsPerGene.out.tab10_MA_C_S3_L001_ReadsPerGene.out.tab⋯47_RZ_P_S50_L003_ReadsPerGene.out.tab47_RZ_P_S50_L004_ReadsPerGene.out.tab9_MA_C_S2_L001_ReadsPerGene.out.tab9_MA_C_S2_L002_ReadsPerGene.out.tab9_MA_C_S2_L003_ReadsPerGene.out.tab9_MA_C_S2_L004_ReadsPerGene.out.tab9_RZ_C_S10_L001_ReadsPerGene.out.tab9_RZ_C_S10_L002_ReadsPerGene.out.tab9_RZ_C_S10_L003_ReadsPerGene.out.tab9_RZ_C_S10_L004_ReadsPerGene.out.tab
<chr><int><int><int><int><int><int><int><int><int>⋯<int><int><int><int><int><int><int><int><int><int>
N_unmapped 2690 2684 2672 2585 7218 7022 7355 707638278⋯ 10036 9962 2245 2291 2276 2105 3386 3612 3283 4853
N_multimapping6610065234665386506639584838807940133839549064124⋯ 536339 52936376258744907637075176149388149618150874156664
N_noFeature 2034720004205492050576814675565477774977371228540⋯1055322104717127956276382837228459503625504801510186525524
N_ambiguous 647 652 697 616 1431 1337 1425 1322 147⋯ 1260 1236 903 848 943 838 1354 1333 1359 1357
 
[20]:

temp1 %>% gather(expid, value, -gene) -> temp2

head(temp2)
A tibble: 6 × 3
geneexpidvalue
<chr><chr><int>
N_unmapped 1_MA_J_S18_L001_ReadsPerGene.out.tab 2690
N_multimapping1_MA_J_S18_L001_ReadsPerGene.out.tab66100
N_noFeature 1_MA_J_S18_L001_ReadsPerGene.out.tab20347
N_ambiguous 1_MA_J_S18_L001_ReadsPerGene.out.tab 647
N_unmapped 1_MA_J_S18_L002_ReadsPerGene.out.tab 2684
N_multimapping1_MA_J_S18_L002_ReadsPerGene.out.tab65234
 
[21]:

temp2 %>%  spread(gene, value) -> temp3

head(temp3)
A tibble: 6 × 5
expidN_ambiguousN_multimappingN_noFeatureN_unmapped
<chr><int><int><int><int>
1_MA_J_S18_L001_ReadsPerGene.out.tab 647 66100 203472690
1_MA_J_S18_L002_ReadsPerGene.out.tab 652 65234 200042684
1_MA_J_S18_L003_ReadsPerGene.out.tab 697 66538 205492672
1_MA_J_S18_L004_ReadsPerGene.out.tab 616 65066 205052585
1_RZ_J_S26_L001_ReadsPerGene.out.tab14313958487681467218
1_RZ_J_S26_L002_ReadsPerGene.out.tab13373880797556547022
 
[22]:

temp3 %>% rename_(.dots = setNames(names(.), c("expid", "namb", "nmulti", "nnofeat","nunmap"))) %>% head
A tibble: 6 × 5
expidnambnmultinnofeatnunmap
<chr><int><int><int><int>
1_MA_J_S18_L001_ReadsPerGene.out.tab 647 66100 203472690
1_MA_J_S18_L002_ReadsPerGene.out.tab 652 65234 200042684
1_MA_J_S18_L003_ReadsPerGene.out.tab 697 66538 205492672
1_MA_J_S18_L004_ReadsPerGene.out.tab 616 65066 205052585
1_RZ_J_S26_L001_ReadsPerGene.out.tab14313958487681467218
1_RZ_J_S26_L002_ReadsPerGene.out.tab13373880797556547022
 
[23]:

### Gather and spread the genes to get a count matrix
out %>%
    dplyr::slice(-(1:4)) %>%
    gather(expid, value, -gene) %>%
    spread(gene, value) -> genecounts

[24]:

genecounts[1:5,1:5]
A tibble: 5 × 5
expidCNAG_00001CNAG_00002CNAG_00003CNAG_00004
<chr><int><int><int><int>
1_MA_J_S18_L001_ReadsPerGene.out.tab0663874
1_MA_J_S18_L002_ReadsPerGene.out.tab0592579
1_MA_J_S18_L003_ReadsPerGene.out.tab0742779
1_MA_J_S18_L004_ReadsPerGene.out.tab0662269
1_RZ_J_S26_L001_ReadsPerGene.out.tab0501651
 
[25]:

out %>%
    dplyr::slice(-(1:4)) %>% t() -> check

[26]:

check[1:5,1:5]
A matrix: 5 × 5 of type chr
geneCNAG_04548CNAG_07303CNAG_07304CNAG_00001CNAG_07305
1_MA_J_S18_L001_ReadsPerGene.out.tab 0 0 8 0 0
1_MA_J_S18_L002_ReadsPerGene.out.tab 0 0 7 0 1
1_MA_J_S18_L003_ReadsPerGene.out.tab 0 0 10 0 0
1_MA_J_S18_L004_ReadsPerGene.out.tab 1 0 9 0 0 
[ ]: