File I/O

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

Strucrure of count data from STAR

• column 1: gene ID

• column 2: counts for unstranded RNA-seq

• column 3: counts for the 1st read strand aligned with RNA (htseq-count option -s yes)

• column 4: counts for the 2nd read strand aligned with RNA (htseq-count option -s reverse)

Get list of files

count_dir <- "/data/hts_2019_data/hts2019_pilot_counts"

files <- list.files(path=count_dir, pattern="*.tab", full.names = TRUE)

files[1]

'/data/hts_2019_data/hts2019_pilot_counts/1_2019_P_M1_S1_L001_ReadsPerGene.out.tab'

Read in first 2 files

d1 <- read_tsv(files[1], col_names=FALSE, col_type=cols())

d2 <- read_tsv(files[2], col_names=FALSE, col_type=cols())

d1 %>% head(5)

A tibble: 5 × 4

X1	X2	X3	X4
<chr>	<dbl>	<dbl>	<dbl>
N_unmapped	21162	21162	21162
N_multimapping	95540	95540	95540
N_noFeature	17075	4260222	26923
N_ambiguous	418493	1416	606
CNAG_04548	0	0	0

d2 %>% head(5)

A tibble: 5 × 4

X1	X2	X3	X4
<chr>	<dbl>	<dbl>	<dbl>
N_unmapped	19644	19644	19644
N_multimapping	94070	94070	94070
N_noFeature	16407	4195892	26053
N_ambiguous	412656	1410	638
CNAG_04548	0	0	0

Combining files

bind_cols(d1, d2) %>% head(5)

A tibble: 5 × 8

X1	X2	X3	X4	X11	X21	X31	X41
<chr>	<dbl>	<dbl>	<dbl>	<chr>	<dbl>	<dbl>	<dbl>
N_unmapped	21162	21162	21162	N_unmapped	19644	19644	19644
N_multimapping	95540	95540	95540	N_multimapping	94070	94070	94070
N_noFeature	17075	4260222	26923	N_noFeature	16407	4195892	26053
N_ambiguous	418493	1416	606	N_ambiguous	412656	1410	638
CNAG_04548	0	0	0	CNAG_04548	0	0	0

bind_rows(d1, d2) %>% head(5)

A tibble: 5 × 4

X1	X2	X3	X4
<chr>	<dbl>	<dbl>	<dbl>
N_unmapped	21162	21162	21162
N_multimapping	95540	95540	95540
N_noFeature	17075	4260222	26923
N_ambiguous	418493	1416	606
CNAG_04548	0	0	0

inner_join(x=d1, y=d2, by='X1') %>% head(5)

A tibble: 5 × 7

X1	X2.x	X3.x	X4.x	X2.y	X3.y	X4.y
<chr>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
N_unmapped	21162	21162	21162	19644	19644	19644
N_multimapping	95540	95540	95540	94070	94070	94070
N_noFeature	17075	4260222	26923	16407	4195892	26053
N_ambiguous	418493	1416	606	412656	1410	638
CNAG_04548	0	0	0	0	0	0

Custom function to read files

read_file <- function(file) {
    read_tsv(file, col_names = FALSE, skip=4, col_types = cols()) %>%
    mutate(source=basename(file)) %>%
    select(source=source, gene=X1, count=X4) %>%
    separate(source, sep='_', into=c("sample", "year", "p", "method", "s", "lane", "junk")) %>%
    select(sample, year, p, method, s, lane, gene, count)
}

read_file(files[1]) %>% head(5)

A tibble: 5 × 8

sample	year	p	method	s	lane	gene	count
<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<dbl>
1	2019	P	M1	S1	L001	CNAG_04548	0
1	2019	P	M1	S1	L001	CNAG_07303	0
1	2019	P	M1	S1	L001	CNAG_07304	14
1	2019	P	M1	S1	L001	CNAG_00001	0
1	2019	P	M1	S1	L001	CNAG_07305	1