Introduction

This page will walk through the steps to go from the raw bulk RNA sequencing fastq files to a count matrix. This involves running a Snakemake pipeline located in the code directory and some code in R. The end product is a count matrix corresponding to genes x samples containing information from all 18 bulk samples from the study. We will also load in sample meta information from a .csv file in this repository.

Run the Snakemake Pipeline after downloading necessary files and creating the conda environment (see README.txt)

1. Move the fastq files from all samples (these are output files from bulk RNA sequencing runs) into the folder code/bulkRNA_preprocessing/fastq/. Undetermined data files are not required.

2. Modify the code/bulkRNA_preprocessing/Pipeline/cluster_solo.json file to correspond to the computing cluster you are working with.

3. Create a STAR genome index for hg38 and point to its location in the config_bulk_hg38.yaml file.

4. Download a GTF file for hg38 and point to its location in the config_bulk_hg38.yaml file.

5. Install the conda working by running “conda env create –file environment.yaml”

6. Run “source activate bulkRNAseq” to activate the bulkRNAseq conda environment.

7. Run the file “submit.sh”.

Load files into R

After running the Snakemake pipeline, you will have output files (located in code/bulkRNA_preprocessing/out) corresponding to the count matrix and quality control of the raw fastq files. In this repository is also included sample information in a csv. We will now load the count matrix and sample information into R and save them as R objects for subsequent analysis. # Load in bulk RNA sequencing data from pilot experiment and visualize correlations between samples

library("gplots")

Attaching package: 'gplots'

The following object is masked from 'package:stats':

    lowess

library("ggplot2")
library("RColorBrewer")
library("scales")
library("edgeR")

Loading required package: limma

#These sample names correspond to the individual_replicate_treatmentCondition.
samplenames <- c("18855_3_S","19160_3_S", "18856_3_U",
                 "18856_1_U","18855_2_S", "18856_2_S",
                 "19160_3_U","18855_2_U", "19160_2_S",
                 "18855_1_S","18856_1_S", "19160_1_S",
                 "19160_2_U","19160_1_U", "18855_1_U",
                 "18856_3_S","18856_2_U", "18855_3_U"
                 )

# Load colors 
colors <- colorRampPalette(c(brewer.pal(9, "Blues")[1],brewer.pal(9, "Blues")[9]))(100)
pal <- c(brewer.pal(9, "Set1"), brewer.pal(8, "Set2"), brewer.pal(12, "Set3"))

# load in counts
raw_counts <- read.table("code/bulkRNA_preprocessing/out/counts/counts.txt", header = T)

# assign row.names based on the Geneid column
row.names(raw_counts) <- raw_counts$Geneid

# remove extra data columns from count matrix
raw_counts <- raw_counts[, -c(1:6)] 

# assign col.names 
names(raw_counts) <- samplenames

#Sample information
sampleinfo <- read.csv("data/Sample.info.RNAseq.csv")
head(sampleinfo)

  Sample_ID Individual Sex Replicate treatment  RIN LibraryPrepBatch
1 18855_1_S   NA18855    F         1    Strain  9.9                1
2 18856_1_S   NA18856    M         1    Strain  8.5                1
3 19160_1_S   NA19160    M         1    Strain  9.9                1
4 18855_2_S   NA18855    F         2    Strain 10.0                1
5 18856_2_S   NA18856    M         2    Strain  9.6                1
6 19160_2_S   NA19160    M         2    Strain 10.0                1
   LibSize
1 19712744
2 13571439
3 20979958
4 18256163
5 16782764
6 16622407

#Re-order sample information to match the columns of count matrix
sampleinfo <- sampleinfo[match(samplenames, sampleinfo$Sample_ID),]

Based on the QC metrics, one of our samples did not pass with high quality and is therefore removed from subsequent analyses. We will, however, keep a copy of the full data to demonstrate that the sample does cluster as expected with its corresponding technical replicates (NA19160 Unstrain, first Replicate).

# save relabeled raw count data matrix
saveRDS(raw_counts, "data/raw_counts_relabeled_includeNA19160.rds")
# save reorganized sample information
saveRDS(sampleinfo, "data/Sample.info.RNAseq.reordered_includeNA19160.rds")

# remove 02 sample that does not pass QC. More information in manuscript.
raw_counts <- raw_counts[, -14] 
#remove 02 sample that does not pass QC
sampleinfo <- sampleinfo[-14,]

Save the raw count and sample information data

# save relabeled raw count data matrix
saveRDS(raw_counts, "data/raw_counts_relabeled.rds")
# save reorganized sample information
saveRDS(sampleinfo, "data/Sample.info.RNAseq.reordered.rds")

Session information

sessionInfo()

R version 3.6.1 (2019-07-05)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Scientific Linux 7.4 (Nitrogen)

Matrix products: default
BLAS/LAPACK: /software/openblas-0.2.19-el7-x86_64/lib/libopenblas_haswellp-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] edgeR_3.26.5       limma_3.40.6       scales_1.1.1      
[4] RColorBrewer_1.1-2 ggplot2_3.3.3      gplots_3.0.1.1    

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.5         compiler_3.6.1     pillar_1.4.7      
 [4] later_1.1.0.1      git2r_0.26.1       workflowr_1.6.2   
 [7] bitops_1.0-6       tools_3.6.1        digest_0.6.27     
[10] lattice_0.20-41    evaluate_0.14      lifecycle_0.2.0   
[13] tibble_3.0.4       gtable_0.3.0       pkgconfig_2.0.3   
[16] rlang_0.4.10       yaml_2.2.1         xfun_0.8          
[19] withr_2.3.0        stringr_1.4.0      dplyr_1.0.2       
[22] knitr_1.23         generics_0.0.2     fs_1.3.1          
[25] vctrs_0.3.6        gtools_3.8.1       caTools_1.17.1.2  
[28] locfit_1.5-9.1     tidyselect_1.1.0   rprojroot_2.0.2   
[31] grid_3.6.1         glue_1.4.2         R6_2.5.0          
[34] rmarkdown_1.13     gdata_2.18.0       purrr_0.3.4       
[37] magrittr_2.0.1     whisker_0.3-2      promises_1.1.1    
[40] htmltools_0.5.0    ellipsis_0.3.1     colorspace_2.0-0  
[43] httpuv_1.5.1       KernSmooth_2.23-15 stringi_1.4.6     
[46] munsell_0.5.0      crayon_1.3.4