---
title: 'Feature Counts'
output:
    html_document:
      number_sections: true
      toc: true
      theme: cosmo
      highlight: tango
---

```{r setup, include=FALSE, warning=FALSE, message=FALSE}
knitr::opts_chunk$set(
  echo = opt$echo,
  error = TRUE
)
```


# User input

```{r 'user input'}
user_input = data.frame(name = names(opt)[-1],
                        value = unlist(opt))
datatable(user_input, rownames = FALSE)
```

# Calculate feature counts

```{r 'ste[ 2'}
res = featureCounts(
  files = strsplit(opt$input_bam_paths, ',')[[1]],
  # annotation
	annot.inbuilt=opt$annot_inbuilt,
	annot.ext=opt$annot_ext,
	isGTFAnnotationFile=opt$isGTFAnnotationFile,
	GTF.featureType=opt$gtf_feature_type,
	GTF.attrType=opt$gtf_attr_type,
	chrAliases=opt$chr_aliases,
	
	# level of summarization
	useMetaFeatures=opt$use_meta_features,
	
	# overlap between reads and features
	allowMultiOverlap=opt$allow_multi_overlap,
	minOverlap=opt$min_overlap,
	largestOverlap=opt$largest_overlap,
	readExtension5=opt$read_extension_5,
	readExtension3=opt$read_extension_3,
	read2pos=opt$read_2_pos,
	
	# multi-mapping reads
	countMultiMappingReads=opt$count_multi_mapping_reads,
	fraction=opt$fraction,

	# read filtering
	minMQS=opt$min_mqs,
	splitOnly=opt$split_only,
	nonSplitOnly=opt$non_split_only,
	primaryOnly=opt$primary_only,
	ignoreDup=opt$ignore_dup,
	
	# strandness
	strandSpecific=opt$strand_specific,
	
	# exon-exon junctions
	juncCounts=opt$junc_counts,
	genome=opt$genome,
	
	# parameters specific to paired end reads
	isPairedEnd=opt$is_paired_end,
	requireBothEndsMapped=opt$require_both_ends_mapped,
	checkFragLength=opt$check_frag_length,
	minFragLength=opt$min_frag_length,
	maxFragLength=opt$max_frag_length,
	countChimericFragments=opt$count_chimeric_fragments,	
	autosort=opt$auto_sort,
	
	# miscellaneous
	nthreads=opt$n_threads,
	maxMOp=opt$max_mop,
	reportReads=opt$report_reads
)
```

# Write counts into CSV file

```{r}
colnames(res$counts) = strsplit(opt$input_bam_names, ',')[[1]]
# write count into csv file
write.table(res$counts, file = 'feature_counts.txt')
```

Display the first 100 rows.

```{r}
datatable(head(res$counts, 100))
```

# Save results into RData file

```{r}
save(opt, res, file = 'feature_counts.RData')
str(res)
```