Merge pull request #84 from bcbio/fix_singleend_lp

Fix single cell read data metrics.

Author: lpantano

inst/templates/rnaseq/00_libs/load_data.R

@@ -1,7 +1,10 @@
 library(tidyverse)
 library(SummarizedExperiment)
 library(janitor)
-load_metrics <- function(se_object, multiqc_data_dir, gtf_fn, counts){
+load_metrics <- function(se=se_object, multiqc=multiqc_data_dir,
+                         gtf=gtf_fn,
+                         counts=counts,
+                         single_end=FALSE){
 
   # bcbio input
   if (!is.na(se_object)){
@@ -32,14 +35,21 @@ load_metrics <- function(se_object, multiqc_data_dir, gtf_fn, counts){
       summarize(total_reads = sum(fastqc_raw_total_sequences))
 
     # This renames to user-friendly names the metrics columns
+    if (single_end){
+      metrics <- metrics %>%
+        dplyr::filter(!is.na(fastqc_raw_total_sequences))
+    }else{
+      metrics <- metrics %>%
+        dplyr::filter(is.na(fastqc_raw_total_sequences))
+    }
     metrics <- metrics %>%
-      dplyr::filter(is.na(fastqc_raw_total_sequences)) %>%
       remove_empty(which = 'cols') %>%
       full_join(total_reads) %>%
       mutate(mapped_reads = samtools_reads_mapped) %>%
-      mutate(exonic_rate = exonic/(star_uniquely_mapped * 2)) %>%
-      mutate(intronic_rate = intronic/(star_uniquely_mapped * 2)) %>%
-      mutate(intergenic_rate = intergenic/(star_uniquely_mapped * 2)) %>%
+      rowwise() %>%
+      mutate(exonic_rate = exonic/(exonic + intronic + intergenic)) %>%
+      mutate(intronic_rate = intronic/(exonic + intronic + intergenic)) %>%
+      mutate(intergenic_rate = intergenic/(exonic + intronic + intergenic)) %>%
       mutate(x5_3_bias = qualimap_5_3_bias)
 
     # Sometimes we don't have rRNA due to mismatch annotation, We skip this if is the case

inst/templates/rnaseq/01_quality_assesment/QC.Rmd

@@ -17,12 +17,13 @@ output:
 params:
   # Fill this file with the right paths to nfcore output
   # Put hg38, mm10, mm39, or other
-  # params_file: params_qc_nf-core-example.R # example data
-  params_file: params_qc_nf-core-example.R
+  # params_file: params_qc-example.R # example data
+  params_file: params_qc-example.R
   genome: hg38
+  single_end: false
+  factor_of_interest: sample_type
   project_file: ../information.R
   functions_file: ../00_libs/load_data.R
-  factor_of_interest: sample_type
 ---
 
 Template developed with materials from https://hbctraining.github.io/main/.
@@ -45,6 +46,11 @@ This code is in this ![](https://img.shields.io/badge/status-stable-green) revis
 #    this is used to color plots, it needs to be part of the metadata
 factor_of_interest=params$factor_of_interest
 genome=params$genome
+<<<<<<< HEAD
+singleend=params$singleend
+=======
+single_end=params$single_end
+>>>>>>> c37fcd42e9e4ba276293a9bc591e6298affd14e5
 # 2. Set input files in this file
 source(params$params_file)
 # 3. If you set up this file, project information will be printed below and
@@ -120,7 +126,8 @@ coldata$sample=row.names(coldata)
 counts <- load_counts(counts_fn)
 counts <- counts[,colnames(counts) %in% coldata$sample]
 
-metrics <- load_metrics(se_object, multiqc_data_dir, gtf_fn, counts) %>% 
+metrics <- load_metrics(se_object, multiqc_data_dir,
+                        gtf_fn, counts, single_end) %>% 
   left_join(coldata, by = c('sample')) %>% 
   as.data.frame()
 metrics <- subset(metrics, metrics$sample %in% coldata$sample)

inst/templates/rnaseq/01_quality_assesment/params_qc-example.R

@@ -4,7 +4,6 @@
 # Example data
 coldata_fn='https://raw.githubusercontent.com/bcbio/bcbioR-test-data/devel/rnaseq/nf-core/coldata.csv'
 counts_fn=url('https://raw.githubusercontent.com/bcbio/bcbioR-test-data/devel/rnaseq/nf-core/star_salmon/salmon.merged.gene_counts.tsv')
-se_object=url('https://raw.githubusercontent.com/bcbio/bcbioR-test-data/devel/rnaseq/nf-core/star_salmon/salmon.merged.gene_counts.rds')
 # This folder is in the output directory inside multiqc folder
 multiqc_data_dir='https://raw.githubusercontent.com/bcbio/bcbioR-test-data/devel/rnaseq/nf-core/multiqc/star_salmon/multiqc-report-data/'
 # This file is inside the genome folder in the output directory