--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/rmarkdown_feature_counts.Rmd	Fri Dec 29 15:03:18 2017 -0500
@@ -0,0 +1,104 @@
+---
+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'}
+opt
+```
+
+# 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(res, file = 'feature_counts.RData')
+str(res)
+```
+
+