DESeq2: Create count matrix and metadata data frame

Objective: create a count matrix and metadata frame for the 2019 pilot data

Load packages

[1]:

library(tidyverse)

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()

Import Counts

[2]:

### Due to time constraints, we skip the import of the count files (see HTS tidyverse notebook)

cntfile <- '/home/jovyan/work/scratch/analysis_output/out/hts-pilot-2019.RData'
attach(cntfile)
tools::md5sum(cntfile)
ls(2)
/home/jovyan/work/scratch/analysis_output/out/hts-pilot-2019.RData: 'f617195665ef950ffb36a89f2d18a9ca'
  1. 'genecounts'
  2. 'mapresults'
[3]:

### Look at the first five genes across the libraries

genecounts[1:8, 1:5]
dim(genecounts)
A tibble: 8 × 5
expidCNAG_00001CNAG_00002CNAG_00003CNAG_00004
<chr><int><int><int><int>
1_2019_P_M1_S1_L001_ReadsPerGene.out.tab 03548223
1_2019_P_M1_S1_L002_ReadsPerGene.out.tab 04346227
1_2019_P_M1_S1_L003_ReadsPerGene.out.tab 04649232
1_2019_P_M1_S1_L004_ReadsPerGene.out.tab 03458222
10_2019_P_M1_S10_L001_ReadsPerGene.out.tab03036130
10_2019_P_M1_S10_L002_ReadsPerGene.out.tab03737117
10_2019_P_M1_S10_L003_ReadsPerGene.out.tab02931135
10_2019_P_M1_S10_L004_ReadsPerGene.out.tab02327131
  1. 144
  2. 8500

Filter out 2018 libraries from count file

[4]:

genecounts %>%
    filter(str_detect(expid, "2019")) ->
        genecounts2019

dim(genecounts2019)
  1. 96
  2. 8500

Import metadata file

[5]:

metadtfile <- '/data/hts_2019_data/hts2019_pilot_rawdata/2019_pilot_metadata.tsv'
tools::md5sum(metadtfile)

mtdf<-readr::read_tsv(metadtfile)

mtdf
/data/hts_2019_data/hts2019_pilot_rawdata/2019_pilot_metadata.tsv: '757e63c2306c3c68a70b81dfcdbf7f22'
Parsed with column specification:
cols(
  .default = col_character(),
  sample_year = col_double(),
  enrich_rep = col_double(),
  RNA_sample_num = col_double(),
  library_num = col_double(),
  bio_replicate = col_double(),
  Nanodrop_260_280 = col_double(),
  Nanodrop_260_230 = col_double(),
  Nanodrop_concentration_ng_ul = col_double(),
  Bioanalyzer_concentration_ng_ul = col_double(),
  RIN_lowered_threshold = col_double()
)
See spec(...) for full column specifications.
A spec_tbl_df: 36 × 22
Labelsample_yeargroupenrich_repRNA_sample_numgenotypeconditionlibprep_personenrichment_methodenrichment_short⋯i5_primeri7_primerlibrary_numbio_replicateNanodrop_260_280Nanodrop_260_230Nanodrop_concentration_ng_ulBioanalyzer_concentration_ng_ulRIN_normal_thresholdRIN_lowered_threshold
<chr><dbl><chr><dbl><dbl><chr><chr><chr><chr><chr>⋯<chr><chr><dbl><dbl><dbl><dbl><dbl><dbl><chr><dbl>
1_2019_P_M1 2019P1 1WT pH4CmRNA M⋯i501i701 112.141.52293 197N/A 9.8
2_2019_P_M1 2019P1 2WT pH4CmRNA M⋯i502i701 222.121.79290 225N/A 9.9
3_2019_P_M1 2019P1 3WT pH4CmRNA M⋯i503i701 332.112.49302 241N/A 9.9
4_2019_P_M1 2019P1 4WT pH4PmRNA M⋯i504i701 442.131.15296 189N/A 9.7
5_2019_P_M1 2019P1 5WT pH4PmRNA M⋯i505i701 552.092.42337 26810 10.0
6_2019_P_M1 2019P1 6WT pH4PmRNA M⋯i506i701 662.082.40319 27610 10.0
7_2019_P_M1 2019P1 7sre1dpH4CmRNA M⋯i507i701 712.132.23232 127N/A 9.9
8_2019_P_M1 2019P1 8sre1dpH4CmRNA M⋯i508i701 822.082.24320 311N/A10.0
9_2019_P_M1 2019P1 9sre1dpH4CmRNA M⋯i501i702 932.091.46342 3269.6 9.6
10_2019_P_M12019P110sre1dpH4PmRNA M⋯i502i7021042.161.25262 16810 10.0
11_2019_P_M12019P111sre1dpH4PmRNA M⋯i503i7021152.122.28212 23610 10.0
12_2019_P_M12019P112sre1dpH4PmRNA M⋯i504i7021262.122.57246 31210 10.0
13_2019_P_M12019P113WT pH8CmRNA M⋯i505i7021311.932.28465 436N/A10.0
14_2019_P_M12019P114WT pH8CmRNA M⋯i506i7021422.131.88170 99N/A 9.9
15_2019_P_M12019P115WT pH8CmRNA M⋯i507i7021532.080.97300 190N/A10.0
16_2019_P_M12019P116WT pH8PmRNA M⋯i508i7021642.152.13236 129N/A10.0
17_2019_P_M12019P117WT pH8PmRNA M⋯i501i7031752.130.89233 111N/A10.0
18_2019_P_M12019P118WT pH8PmRNA M⋯i502i7031862.142.01213 134N/A10.0
19_2019_P_M12019P119sre1dpH8CmRNA M⋯i503i7031912.142.19236 80N/A10.0
20_2019_P_M12019P120sre1dpH8CmRNA M⋯i504i7032022.142.39172 131N/A10.0
21_2019_P_M12019P121sre1dpH8CmRNA M⋯i505i7032132.012.22397 295N/A10.0
22_2019_P_M12019P122sre1dpH8PmRNA M⋯i506i7032242.222.33567 49910 10.0
23_2019_P_M12019P123sre1dpH8PmRNA M⋯i507i7032352.132.56254 206N/A10.0
24_2019_P_M12019P124sre1dpH8PmRNA M⋯i508i7032462.112.16215 33610 10.0
2_2018_P_H1 2018P1 2WT YPDJRNaseHH⋯i501i704252 NA NA NA212010 NA
3_2018_P_H1 2018P1 3WT YPDJRNaseHH⋯i502i704263 NA NA NA150010 NA
4_2018_P_H1 2018P1 4WT YPDJRNaseHH⋯i503i704274 NA NA NA149010 NA
2_2018_P_H2 2018P2 2WT YPDJRNaseHH⋯i504i704282 NA NA NA212010 NA
3_2018_P_H2 2018P2 3WT YPDJRNaseHH⋯i505i704293 NA NA NA150010 NA
4_2018_P_H2 2018P2 4WT YPDJRNaseHH⋯i506i704304 NA NA NA149010 NA
2_2018_P_T1 2018P1 2WT YPDJTotal T⋯i507i704312 NA NA NA212010 NA
3_2018_P_T1 2018P1 3WT YPDJTotal T⋯i508i704323 NA NA NA150010 NA
4_2018_P_T1 2018P1 4WT YPDJTotal T⋯i501i705334 NA NA NA149010 NA
2_2018_P_M1 2018P1 2WT YPDJmRNA M⋯i502i705342 NA NA NA212010 NA
3_2018_P_M1 2018P1 3WT YPDJmRNA M⋯i503i705353 NA NA NA150010 NA
4_2018_P_M1 2018P1 4WT YPDJmRNA M⋯i504i705364 NA NA NA149010 NA

filter out 2018 samples

[6]:

mtdf %>%
    filter(str_detect(Label, "2019")) ->
        mtdf2019

dim(mtdf2019)
  1. 24
  2. 22

Add a label column to the counts object so that it can be merged with the metadata

[7]:

myregex <- "_S[1-9][0-9]{0,2}_L00[1-4]_ReadsPerGene.out.tab"


### dim before adding label column
dim(genecounts2019)

genecounts2019 %>%
    mutate(Label=str_replace(expid, myregex, "")) ->
        genecounts2019

### dim before adding label column

dim(genecounts2019)
  1. 96
  2. 8500
  1. 96
  2. 8501
[8]:

### Look at the first five genes across the libraries

genecounts2019[1:8, c(1:5, ncol(genecounts2019))]
dim(genecounts2019)
A tibble: 8 × 6
expidCNAG_00001CNAG_00002CNAG_00003CNAG_00004Label
<chr><int><int><int><int><chr>
1_2019_P_M1_S1_L001_ReadsPerGene.out.tab 035482231_2019_P_M1
1_2019_P_M1_S1_L002_ReadsPerGene.out.tab 043462271_2019_P_M1
1_2019_P_M1_S1_L003_ReadsPerGene.out.tab 046492321_2019_P_M1
1_2019_P_M1_S1_L004_ReadsPerGene.out.tab 034582221_2019_P_M1
10_2019_P_M1_S10_L001_ReadsPerGene.out.tab0303613010_2019_P_M1
10_2019_P_M1_S10_L002_ReadsPerGene.out.tab0373711710_2019_P_M1
10_2019_P_M1_S10_L003_ReadsPerGene.out.tab0293113510_2019_P_M1
10_2019_P_M1_S10_L004_ReadsPerGene.out.tab0232713110_2019_P_M1
  1. 96
  2. 8501

Add counts across lanes

[9]:

genecounts2019 %>%
    group_by(Label) %>%
        summarize_each(sum, -expid) ->
            cnt2019

[10]:

cnt2019[1:4, 1:5]
A tibble: 4 × 5
LabelCNAG_00001CNAG_00002CNAG_00003CNAG_00004
<chr><int><int><int><int>
1_2019_P_M1 0158201904
10_2019_P_M10119131513
11_2019_P_M10 90121573
12_2019_P_M10 81151533

Check to make sure that every label in your count file has an entry in the metadata file

[11]:

setdiff(cnt2019$Label, mtdf2019$Label)
setdiff(mtdf2019$Label, cnt2019$Label)

Save objects to image file

[13]:

curdir <- "/home/jovyan/work/scratch/analysis_output"
imgdir <- file.path(curdir, "img")

imgfile <- file.path(imgdir, "pilotcnt2019.RData")

imgfile
'/home/jovyan/work/scratch/analysis_output/img/pilotcnt2019.RData'
[14]:

save(cnt2019, mtdf2019, file = imgfile)

tools::md5sum(imgfile)
/home/jovyan/work/scratch/analysis_output/img/pilotcnt2019.RData: '0d1cdfc8bb97baf52ca4f8f33e2c0ae4'
[15]:

sessionInfo()

R version 3.6.0 (2019-04-26)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux 9 (stretch)

Matrix products: default
BLAS:   /usr/lib/openblas-base/libblas.so.3
LAPACK: /usr/lib/libopenblasp-r0.2.19.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C
 [9] LC_ADDRESS=C               LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] forcats_0.4.0   stringr_1.4.0   dplyr_0.8.1     purrr_0.3.2
[5] readr_1.3.1     tidyr_0.8.3     tibble_2.1.2    ggplot2_3.1.1
[9] tidyverse_1.2.1

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.1       cellranger_1.1.0 plyr_1.8.4       pillar_1.4.1
 [5] compiler_3.6.0   base64enc_0.1-3  tools_3.6.0      zeallot_0.1.0
 [9] digest_0.6.19    uuid_0.1-2       lubridate_1.7.4  jsonlite_1.6
[13] evaluate_0.14    nlme_3.1-139     gtable_0.3.0     lattice_0.20-38
[17] pkgconfig_2.0.2  rlang_0.3.4      cli_1.1.0        rstudioapi_0.10
[21] IRdisplay_0.7.0  IRkernel_1.0.1   haven_2.1.0      withr_2.1.2
[25] xml2_1.2.0       httr_1.4.0       repr_1.0.1       vctrs_0.1.0
[29] hms_0.4.2        generics_0.0.2   grid_3.6.0       tidyselect_0.2.5
[33] glue_1.3.1       R6_2.4.0         readxl_1.3.1     pbdZMQ_0.3-3
[37] modelr_0.1.4     magrittr_1.5     backports_1.1.4  scales_1.0.0
[41] htmltools_0.3.6  rvest_0.3.4      assertthat_0.2.1 colorspace_1.4-1
[45] stringi_1.4.3    lazyeval_0.2.2   munsell_0.5.0    broom_0.5.2
[49] crayon_1.3.4