add fpkm export to featcounts_deseq

dge_workflow/featcounts_deseq.R

@@ -358,6 +358,50 @@ geneInfo <- quote({
         biomaRt::getBM(mart=mart)
 }) %>% cache_it("geneInfo")
 
+## save as reference for downstream analysis
+write_tsv(geneInfo, path="geneInfo.txt")
+
+
+#' Export counts per replicate as double-normalized FPKMs=(read_count * 10^9)/(gene_length * library_size)
+fpkms <- countMatrix %>%
+    rownames2column("ensembl_gene_id") %>%
+    gather(replicate, num_tags, -ensembl_gene_id) %>% tbl_df %>%
+    left_join(transmute(geneInfo, ensembl_gene_id, gene_length=end_position-start_position), by="ensembl_gene_id") %>%
+#    left_join(mutate(geneInfo, ensembl_gene_id, gene_length=end_position-start_position), by="ensembl_gene_id") %>%
+## calculate normalized counts
+    group_by(replicate) %>%
+    mutate(lib_size=sum(num_tags), fpkm=((1E9/lib_size)*(num_tags/gene_length)) %>% round(digits=3)) %>%
+    ungroup()
+
+#fpkms %>% ggplot(aes(fpkm)) + geom_histogram() + scale_x_log10(labels=comma)
+
+write_tsv(fpkms, "fpkms_by_replicate.long.txt")
+fpkms %>%
+    transmute(ensembl_gene_id, replicate, fpkm) %>%
+    spread(replicate, fpkm) %>% write_tsv("fpkms_by_replicate.wide.txt")
+
+#' Also export averaged counts per sample as double-normalized FPKMs=(read_count * 10^9)/(gene_length * library_size)
+
+sampleFpkms <- countMatrix %>%
+    rownames2column("ensembl_gene_id") %>%
+    gather(replicate, num_tags, -ensembl_gene_id) %>% tbl_df %>%
+    mutate(sample=get_sample_from_replicate(replicate)) %>%
+    group_by(sample, ensembl_gene_id) %>%
+    summarize(median_tags=median(num_tags)) %>%
+    left_join(transmute(geneInfo, ensembl_gene_id, gene_length=end_position-start_position), by="ensembl_gene_id") %>%
+    group_by(sample) %>%
+    mutate(lib_size=sum(median_tags), fpkm=((1E9/lib_size)*(median_tags/gene_length)) %>% round(digits=3)) %>%
+    ungroup()
+
+#fpkms %>% ggplot(aes(fpkm)) + geom_histogram() + scale_x_log10(labels=comma)
+
+write_tsv(sampleFpkms, "fpkms_by_sample.long.txt")
+sampleFpkms %>%
+    transmute(ensembl_gene_id, sample, fpkm) %>%
+    spread(sample, fpkm) %>% write_tsv("fpkms_by_sample.wide.txt")
+
+#'  [FPKMs as table](fpkms.long.txt) [FPKMs as Matrix](fpkms.wide.txt)
+
 
 # Export the complete dataset for later analysis
 deAnnot <- deResults %>%