|
0
|
1 ---
|
|
11
|
2 title: 'Feature Counts'
|
|
0
|
3 output:
|
|
|
4 html_document:
|
|
|
5 number_sections: true
|
|
|
6 toc: true
|
|
|
7 theme: cosmo
|
|
|
8 highlight: tango
|
|
|
9 ---
|
|
|
10
|
|
|
11 ```{r setup, include=FALSE, warning=FALSE, message=FALSE}
|
|
|
12 knitr::opts_chunk$set(
|
|
7
|
13 echo = opt$echo,
|
|
0
|
14 error = TRUE
|
|
|
15 )
|
|
|
16 ```
|
|
|
17
|
|
|
18
|
|
11
|
19 # User input
|
|
0
|
20
|
|
11
|
21 ```{r 'user input'}
|
|
14
|
22 user_input = data.frame(name = names(opt)[-1],
|
|
|
23 value = unlist(opt))
|
|
|
24 datatable(user_input, rownames = FALSE)
|
|
0
|
25 ```
|
|
|
26
|
|
11
|
27 # Calculate feature counts
|
|
|
28
|
|
0
|
29 ```{r 'ste[ 2'}
|
|
10
|
30 res = featureCounts(
|
|
11
|
31 files = strsplit(opt$input_bam_paths, ',')[[1]],
|
|
10
|
32 # annotation
|
|
|
33 annot.inbuilt=opt$annot_inbuilt,
|
|
|
34 annot.ext=opt$annot_ext,
|
|
|
35 isGTFAnnotationFile=opt$isGTFAnnotationFile,
|
|
|
36 GTF.featureType=opt$gtf_feature_type,
|
|
|
37 GTF.attrType=opt$gtf_attr_type,
|
|
|
38 chrAliases=opt$chr_aliases,
|
|
|
39
|
|
|
40 # level of summarization
|
|
|
41 useMetaFeatures=opt$use_meta_features,
|
|
|
42
|
|
|
43 # overlap between reads and features
|
|
|
44 allowMultiOverlap=opt$allow_multi_overlap,
|
|
|
45 minOverlap=opt$min_overlap,
|
|
|
46 largestOverlap=opt$largest_overlap,
|
|
|
47 readExtension5=opt$read_extension_5,
|
|
|
48 readExtension3=opt$read_extension_3,
|
|
|
49 read2pos=opt$read_2_pos,
|
|
|
50
|
|
|
51 # multi-mapping reads
|
|
|
52 countMultiMappingReads=opt$count_multi_mapping_reads,
|
|
|
53 fraction=opt$fraction,
|
|
0
|
54
|
|
10
|
55 # read filtering
|
|
|
56 minMQS=opt$min_mqs,
|
|
|
57 splitOnly=opt$split_only,
|
|
|
58 nonSplitOnly=opt$non_split_only,
|
|
|
59 primaryOnly=opt$primary_only,
|
|
|
60 ignoreDup=opt$ignore_dup,
|
|
|
61
|
|
|
62 # strandness
|
|
|
63 strandSpecific=opt$strand_specific,
|
|
|
64
|
|
|
65 # exon-exon junctions
|
|
|
66 juncCounts=opt$junc_counts,
|
|
|
67 genome=opt$genome,
|
|
|
68
|
|
|
69 # parameters specific to paired end reads
|
|
|
70 isPairedEnd=opt$is_paired_end,
|
|
|
71 requireBothEndsMapped=opt$require_both_ends_mapped,
|
|
|
72 checkFragLength=opt$check_frag_length,
|
|
|
73 minFragLength=opt$min_frag_length,
|
|
|
74 maxFragLength=opt$max_frag_length,
|
|
|
75 countChimericFragments=opt$count_chimeric_fragments,
|
|
|
76 autosort=opt$auto_sort,
|
|
|
77
|
|
|
78 # miscellaneous
|
|
|
79 nthreads=opt$n_threads,
|
|
|
80 maxMOp=opt$max_mop,
|
|
|
81 reportReads=opt$report_reads
|
|
|
82 )
|
|
11
|
83 ```
|
|
|
84
|
|
|
85 # Write counts into CSV file
|
|
|
86
|
|
|
87 ```{r}
|
|
|
88 colnames(res$counts) = strsplit(opt$input_bam_names, ',')[[1]]
|
|
|
89 # write count into csv file
|
|
13
|
90 write.table(res$counts, file = 'feature_counts.txt')
|
|
|
91 ```
|
|
|
92
|
|
|
93 Display the first 100 rows.
|
|
|
94
|
|
|
95 ```{r}
|
|
|
96 datatable(head(res$counts, 100))
|
|
11
|
97 ```
|
|
|
98
|
|
|
99 # Save results into RData file
|
|
|
100
|
|
|
101 ```{r}
|
|
14
|
102 save(opt, res, file = 'feature_counts.RData')
|
|
10
|
103 str(res)
|
|
0
|
104 ```
|
|
|
105
|
|
|
106
|
