Get Star Count of the Pilot 2019 Data

Load paths and libraries

library(foreach)
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() ──
✖ purrr::accumulate() masks foreach::accumulate()
✖ dplyr::filter()     masks stats::filter()
✖ dplyr::lag()        masks stats::lag()
✖ purrr::when()       masks foreach::when()

data_root_dir <- "/data/hts_2019_data/"
raw_fastq_dir <- paste(data_root_dir, "hts2019_pilot_rawdata", sep = "")
metadata_file <- file.path(raw_fastq_dir, "2019_pilot_metadata.tsv")
count_dir <- paste(data_root_dir, "hts2019_pilot_counts", sep = "")

Create directories

Below is an illustration of our folder structure.

scratch
└── bioinf_intro
└── analysis_output
    ├── out -> the folder to store all our output data files
    └── img -> the folder to store all our images

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

system(paste("mkdir -p", outdir), intern = TRUE)
system(paste("mkdir -p", imgdir), intern = TRUE)

count_suffix <- "_ReadsPerGene.out.tab"

list.files(count_dir, pattern = paste0(count_suffix,"$"), full.names = FALSE) -> countfiles

Load STAR Count Data

list.files(count_dir, pattern = paste0(count_suffix,"$"), full.names = FALSE) -> countfiles
countfiles

1. '1_2019_P_M1_S1_L001_ReadsPerGene.out.tab'
2. '1_2019_P_M1_S1_L002_ReadsPerGene.out.tab'
3. '1_2019_P_M1_S1_L003_ReadsPerGene.out.tab'
4. '1_2019_P_M1_S1_L004_ReadsPerGene.out.tab'
5. '10_2019_P_M1_S10_L001_ReadsPerGene.out.tab'
6. '10_2019_P_M1_S10_L002_ReadsPerGene.out.tab'
7. '10_2019_P_M1_S10_L003_ReadsPerGene.out.tab'
8. '10_2019_P_M1_S10_L004_ReadsPerGene.out.tab'
9. '11_2019_P_M1_S11_L001_ReadsPerGene.out.tab'
10. '11_2019_P_M1_S11_L002_ReadsPerGene.out.tab'
11. '11_2019_P_M1_S11_L003_ReadsPerGene.out.tab'
12. '11_2019_P_M1_S11_L004_ReadsPerGene.out.tab'
13. '12_2019_P_M1_S12_L001_ReadsPerGene.out.tab'
14. '12_2019_P_M1_S12_L002_ReadsPerGene.out.tab'
15. '12_2019_P_M1_S12_L003_ReadsPerGene.out.tab'
16. '12_2019_P_M1_S12_L004_ReadsPerGene.out.tab'
17. '13_2019_P_M1_S13_L001_ReadsPerGene.out.tab'
18. '13_2019_P_M1_S13_L002_ReadsPerGene.out.tab'
19. '13_2019_P_M1_S13_L003_ReadsPerGene.out.tab'
20. '13_2019_P_M1_S13_L004_ReadsPerGene.out.tab'
21. '14_2019_P_M1_S14_L001_ReadsPerGene.out.tab'
22. '14_2019_P_M1_S14_L002_ReadsPerGene.out.tab'
23. '14_2019_P_M1_S14_L003_ReadsPerGene.out.tab'
24. '14_2019_P_M1_S14_L004_ReadsPerGene.out.tab'
25. '15_2019_P_M1_S15_L001_ReadsPerGene.out.tab'
26. '15_2019_P_M1_S15_L002_ReadsPerGene.out.tab'
27. '15_2019_P_M1_S15_L003_ReadsPerGene.out.tab'
28. '15_2019_P_M1_S15_L004_ReadsPerGene.out.tab'
29. '16_2019_P_M1_S16_L001_ReadsPerGene.out.tab'
30. '16_2019_P_M1_S16_L002_ReadsPerGene.out.tab'
31. '16_2019_P_M1_S16_L003_ReadsPerGene.out.tab'
32. '16_2019_P_M1_S16_L004_ReadsPerGene.out.tab'
33. '17_2019_P_M1_S17_L001_ReadsPerGene.out.tab'
34. '17_2019_P_M1_S17_L002_ReadsPerGene.out.tab'
35. '17_2019_P_M1_S17_L003_ReadsPerGene.out.tab'
36. '17_2019_P_M1_S17_L004_ReadsPerGene.out.tab'
37. '18_2019_P_M1_S18_L001_ReadsPerGene.out.tab'
38. '18_2019_P_M1_S18_L002_ReadsPerGene.out.tab'
39. '18_2019_P_M1_S18_L003_ReadsPerGene.out.tab'
40. '18_2019_P_M1_S18_L004_ReadsPerGene.out.tab'
41. '19_2019_P_M1_S19_L001_ReadsPerGene.out.tab'
42. '19_2019_P_M1_S19_L002_ReadsPerGene.out.tab'
43. '19_2019_P_M1_S19_L003_ReadsPerGene.out.tab'
44. '19_2019_P_M1_S19_L004_ReadsPerGene.out.tab'
45. '2_2018_P_H1_S25_L001_ReadsPerGene.out.tab'
46. '2_2018_P_H1_S25_L002_ReadsPerGene.out.tab'
47. '2_2018_P_H1_S25_L003_ReadsPerGene.out.tab'
48. '2_2018_P_H1_S25_L004_ReadsPerGene.out.tab'
49. '2_2018_P_H2_S28_L001_ReadsPerGene.out.tab'
50. '2_2018_P_H2_S28_L002_ReadsPerGene.out.tab'
51. '2_2018_P_H2_S28_L003_ReadsPerGene.out.tab'
52. '2_2018_P_H2_S28_L004_ReadsPerGene.out.tab'
53. '2_2018_P_M1_S34_L001_ReadsPerGene.out.tab'
54. '2_2018_P_M1_S34_L002_ReadsPerGene.out.tab'
55. '2_2018_P_M1_S34_L003_ReadsPerGene.out.tab'
56. '2_2018_P_M1_S34_L004_ReadsPerGene.out.tab'
57. '2_2018_P_T1_S31_L001_ReadsPerGene.out.tab'
58. '2_2018_P_T1_S31_L002_ReadsPerGene.out.tab'
59. '2_2018_P_T1_S31_L003_ReadsPerGene.out.tab'
60. '2_2018_P_T1_S31_L004_ReadsPerGene.out.tab'
61. '2_2019_P_M1_S2_L001_ReadsPerGene.out.tab'
62. '2_2019_P_M1_S2_L002_ReadsPerGene.out.tab'
63. '2_2019_P_M1_S2_L003_ReadsPerGene.out.tab'
64. '2_2019_P_M1_S2_L004_ReadsPerGene.out.tab'
65. '20_2019_P_M1_S20_L001_ReadsPerGene.out.tab'
66. '20_2019_P_M1_S20_L002_ReadsPerGene.out.tab'
67. '20_2019_P_M1_S20_L003_ReadsPerGene.out.tab'
68. '20_2019_P_M1_S20_L004_ReadsPerGene.out.tab'
69. '21_2019_P_M1_S21_L001_ReadsPerGene.out.tab'
70. '21_2019_P_M1_S21_L002_ReadsPerGene.out.tab'
71. '21_2019_P_M1_S21_L003_ReadsPerGene.out.tab'
72. '21_2019_P_M1_S21_L004_ReadsPerGene.out.tab'
73. '22_2019_P_M1_S22_L001_ReadsPerGene.out.tab'
74. '22_2019_P_M1_S22_L002_ReadsPerGene.out.tab'
75. '22_2019_P_M1_S22_L003_ReadsPerGene.out.tab'
76. '22_2019_P_M1_S22_L004_ReadsPerGene.out.tab'
77. '23_2019_P_M1_S23_L001_ReadsPerGene.out.tab'
78. '23_2019_P_M1_S23_L002_ReadsPerGene.out.tab'
79. '23_2019_P_M1_S23_L003_ReadsPerGene.out.tab'
80. '23_2019_P_M1_S23_L004_ReadsPerGene.out.tab'
81. '24_2019_P_M1_S24_L001_ReadsPerGene.out.tab'
82. '24_2019_P_M1_S24_L002_ReadsPerGene.out.tab'
83. '24_2019_P_M1_S24_L003_ReadsPerGene.out.tab'
84. '24_2019_P_M1_S24_L004_ReadsPerGene.out.tab'
85. '3_2018_P_H1_S26_L001_ReadsPerGene.out.tab'
86. '3_2018_P_H1_S26_L002_ReadsPerGene.out.tab'
87. '3_2018_P_H1_S26_L003_ReadsPerGene.out.tab'
88. '3_2018_P_H1_S26_L004_ReadsPerGene.out.tab'
89. '3_2018_P_H2_S29_L001_ReadsPerGene.out.tab'
90. '3_2018_P_H2_S29_L002_ReadsPerGene.out.tab'
91. '3_2018_P_H2_S29_L003_ReadsPerGene.out.tab'
92. '3_2018_P_H2_S29_L004_ReadsPerGene.out.tab'
93. '3_2018_P_M1_S35_L001_ReadsPerGene.out.tab'
94. '3_2018_P_M1_S35_L002_ReadsPerGene.out.tab'
95. '3_2018_P_M1_S35_L003_ReadsPerGene.out.tab'
96. '3_2018_P_M1_S35_L004_ReadsPerGene.out.tab'
97. '3_2018_P_T1_S32_L001_ReadsPerGene.out.tab'
98. '3_2018_P_T1_S32_L002_ReadsPerGene.out.tab'
99. '3_2018_P_T1_S32_L003_ReadsPerGene.out.tab'
100. '3_2018_P_T1_S32_L004_ReadsPerGene.out.tab'
101. '3_2019_P_M1_S3_L001_ReadsPerGene.out.tab'
102. '3_2019_P_M1_S3_L002_ReadsPerGene.out.tab'
103. '3_2019_P_M1_S3_L003_ReadsPerGene.out.tab'
104. '3_2019_P_M1_S3_L004_ReadsPerGene.out.tab'
105. '4_2018_P_H1_S27_L001_ReadsPerGene.out.tab'
106. '4_2018_P_H1_S27_L002_ReadsPerGene.out.tab'
107. '4_2018_P_H1_S27_L003_ReadsPerGene.out.tab'
108. '4_2018_P_H1_S27_L004_ReadsPerGene.out.tab'
109. '4_2018_P_H2_S30_L001_ReadsPerGene.out.tab'
110. '4_2018_P_H2_S30_L002_ReadsPerGene.out.tab'
111. '4_2018_P_H2_S30_L003_ReadsPerGene.out.tab'
112. '4_2018_P_H2_S30_L004_ReadsPerGene.out.tab'
113. '4_2018_P_M1_S36_L001_ReadsPerGene.out.tab'
114. '4_2018_P_M1_S36_L002_ReadsPerGene.out.tab'
115. '4_2018_P_M1_S36_L003_ReadsPerGene.out.tab'
116. '4_2018_P_M1_S36_L004_ReadsPerGene.out.tab'
117. '4_2018_P_T1_S33_L001_ReadsPerGene.out.tab'
118. '4_2018_P_T1_S33_L002_ReadsPerGene.out.tab'
119. '4_2018_P_T1_S33_L003_ReadsPerGene.out.tab'
120. '4_2018_P_T1_S33_L004_ReadsPerGene.out.tab'
121. '4_2019_P_M1_S4_L001_ReadsPerGene.out.tab'
122. '4_2019_P_M1_S4_L002_ReadsPerGene.out.tab'
123. '4_2019_P_M1_S4_L003_ReadsPerGene.out.tab'
124. '4_2019_P_M1_S4_L004_ReadsPerGene.out.tab'
125. '5_2019_P_M1_S5_L001_ReadsPerGene.out.tab'
126. '5_2019_P_M1_S5_L002_ReadsPerGene.out.tab'
127. '5_2019_P_M1_S5_L003_ReadsPerGene.out.tab'
128. '5_2019_P_M1_S5_L004_ReadsPerGene.out.tab'
129. '6_2019_P_M1_S6_L001_ReadsPerGene.out.tab'
130. '6_2019_P_M1_S6_L002_ReadsPerGene.out.tab'
131. '6_2019_P_M1_S6_L003_ReadsPerGene.out.tab'
132. '6_2019_P_M1_S6_L004_ReadsPerGene.out.tab'
133. '7_2019_P_M1_S7_L001_ReadsPerGene.out.tab'
134. '7_2019_P_M1_S7_L002_ReadsPerGene.out.tab'
135. '7_2019_P_M1_S7_L003_ReadsPerGene.out.tab'
136. '7_2019_P_M1_S7_L004_ReadsPerGene.out.tab'
137. '8_2019_P_M1_S8_L001_ReadsPerGene.out.tab'
138. '8_2019_P_M1_S8_L002_ReadsPerGene.out.tab'
139. '8_2019_P_M1_S8_L003_ReadsPerGene.out.tab'
140. '8_2019_P_M1_S8_L004_ReadsPerGene.out.tab'
141. '9_2019_P_M1_S9_L001_ReadsPerGene.out.tab'
142. '9_2019_P_M1_S9_L002_ReadsPerGene.out.tab'
143. '9_2019_P_M1_S9_L003_ReadsPerGene.out.tab'
144. '9_2019_P_M1_S9_L004_ReadsPerGene.out.tab'

length(countfiles)

144

So there are 96 files about 2019 pilot data under this folder. Each file is generated from a fastq file using STAR.

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")
}


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

# Data type for each column
coltypes <- list(col_character(), col_integer(), col_integer(), col_integer())

read the count files and combine them

out <- foreach(file = countfiles, .combine = mycombine) %do% {
    cntfile <- mystarfile(count_dir, file)
    readr::read_tsv(cntfile, col_names = FALSE, col_types = coltypes ) %>%
        dplyr::select(X1, X4) %>%
            dplyr::rename_(.dots=setNames(names(.), c("gene", file)))
}

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.”

dim(out)

1. 8503
2. 145

out[1:6, 1:6]

A tibble: 6 × 6

gene	1_2019_P_M1_S1_L001_ReadsPerGene.out.tab	1_2019_P_M1_S1_L002_ReadsPerGene.out.tab	1_2019_P_M1_S1_L003_ReadsPerGene.out.tab	1_2019_P_M1_S1_L004_ReadsPerGene.out.tab	10_2019_P_M1_S10_L001_ReadsPerGene.out.tab
<chr>	<int>	<int>	<int>	<int>	<int>
N_unmapped	21162	19644	19307	17356	6456
N_multimapping	95540	94070	97787	98723	75280
N_noFeature	26923	26053	27257	26775	26470
N_ambiguous	606	638	615	607	464
CNAG_04548	0	0	0	0	0
CNAG_07303	0	0	0	0	0

Gather and spread the genes to get a count matrix (genecounts)

out %>%
    slice(-(1:4)) %>%
    gather(expid, value, -gene) %>%
    spread(gene, value) -> genecounts

genecounts[1:6, 1:6]

A tibble: 6 × 6

expid	CNAG_00001	CNAG_00002	CNAG_00003	CNAG_00004	CNAG_00005
<chr>	<int>	<int>	<int>	<int>	<int>
1_2019_P_M1_S1_L001_ReadsPerGene.out.tab	0	35	48	223	5
1_2019_P_M1_S1_L002_ReadsPerGene.out.tab	0	43	46	227	7
1_2019_P_M1_S1_L003_ReadsPerGene.out.tab	0	46	49	232	8
1_2019_P_M1_S1_L004_ReadsPerGene.out.tab	0	34	58	222	2
10_2019_P_M1_S10_L001_ReadsPerGene.out.tab	0	30	36	130	5
10_2019_P_M1_S10_L002_ReadsPerGene.out.tab	0	37	37	117	7

Gather and spread the first four rows to nmisc

For nmisc, we will take the first 4 rows of out since those are the summarizing features. Next, we want to transform the data frame so that it is in statistical format (the samples are the rows and the feature types are the columns). Using a combination of gather and spread, we can transpose our matrix into the desired format.

out %>%
    slice(1:4) %>%
    gather(expid, value, -gene) %>%
    spread(gene, value) %>%
    rename_(.dots=setNames(names(.), c("expid", "namb", "nmulti", "nnofeat","nunmap"))) -> nmisc

nmisc[1:6, ]

A tibble: 6 × 5

expid	namb	nmulti	nnofeat	nunmap
<chr>	<int>	<int>	<int>	<int>
1_2019_P_M1_S1_L001_ReadsPerGene.out.tab	606	95540	26923	21162
1_2019_P_M1_S1_L002_ReadsPerGene.out.tab	638	94070	26053	19644
1_2019_P_M1_S1_L003_ReadsPerGene.out.tab	615	97787	27257	19307
1_2019_P_M1_S1_L004_ReadsPerGene.out.tab	607	98723	26775	17356
10_2019_P_M1_S10_L001_ReadsPerGene.out.tab	464	75280	26470	6456
10_2019_P_M1_S10_L002_ReadsPerGene.out.tab	445	73190	25919	6097

Gather and spread the gene rows

For each samples, we want to sum up all the counts, so we can create a variable denoting the number of total genes mapped for each sample by summing across the rows.

out %>%
    slice(-(1:4)) %>%
    gather(expid, value, -gene) %>%
    spread(gene, value) %>%
    mutate(ngenemap=rowSums(.[-1])) %>%
    select(expid, ngenemap) -> ngene

ngene[1:6, ]

A tibble: 6 × 2

expid	ngenemap
<chr>	<dbl>
1_2019_P_M1_S1_L001_ReadsPerGene.out.tab	4660549
1_2019_P_M1_S1_L002_ReadsPerGene.out.tab	4591006
1_2019_P_M1_S1_L003_ReadsPerGene.out.tab	4715846
1_2019_P_M1_S1_L004_ReadsPerGene.out.tab	4681095
10_2019_P_M1_S10_L001_ReadsPerGene.out.tab	3261459
10_2019_P_M1_S10_L002_ReadsPerGene.out.tab	3208423

merge in the 4 misc counts and add summaries

So far, we can create a comprehensive data frame mapresults which will combine ngene with nmisc. This data frame will have summarizing mapping features in addition to proportion features.

ngene %>%
    full_join(nmisc, by="expid") %>%
    mutate(depth = as.integer(ngenemap + namb + nmulti + nnofeat + nunmap)) %>%
    mutate(prob.gene = ngenemap / depth) %>%
    mutate(prob.nofeat = nnofeat / depth) %>%
    mutate(prob.unique = (ngenemap+nnofeat) / depth) -> mapresults

mapresults[1:6, ]

A tibble: 6 × 10

expid	ngenemap	namb	nmulti	nnofeat	nunmap	depth	prob.gene	prob.nofeat	prob.unique
<chr>	<dbl>	<int>	<int>	<int>	<int>	<int>	<dbl>	<dbl>	<dbl>
1_2019_P_M1_S1_L001_ReadsPerGene.out.tab	4660549	606	95540	26923	21162	4804780	0.9699818	0.005603378	0.9755851
1_2019_P_M1_S1_L002_ReadsPerGene.out.tab	4591006	638	94070	26053	19644	4731411	0.9703249	0.005506391	0.9758313
1_2019_P_M1_S1_L003_ReadsPerGene.out.tab	4715846	615	97787	27257	19307	4860812	0.9701766	0.005607499	0.9757841
1_2019_P_M1_S1_L004_ReadsPerGene.out.tab	4681095	607	98723	26775	17356	4824556	0.9702644	0.005549733	0.9758141
10_2019_P_M1_S10_L001_ReadsPerGene.out.tab	3261459	464	75280	26470	6456	3370129	0.9677549	0.007854299	0.9756092
10_2019_P_M1_S10_L002_ReadsPerGene.out.tab	3208423	445	73190	25919	6097	3314074	0.9681205	0.007820888	0.9759414

outfile <- file.path(outdir, "hts-pilot-2019.RData")
save(mapresults, genecounts, file=outfile)
tools::md5sum(outfile)

/home/jovyan/work/scratch/analysis_output/out/hts-pilot-2019.RData: 'f617195665ef950ffb36a89f2d18a9ca'