Hi I was wondering if you could let me know as to what would be the correct way to get TMM normalized region counts.
Here is the code. Kindly let me know which is the correct way to do it :-
library(csaw) options(echo=TRUE) # if you want see commands in output file args <- commandArgs(trailingOnly = TRUE) if(length(args) < 2) { args <- c("--help") } cell_type <- as.character(unlist(args[1])) histone_name <- as.character(unlist(args[2])) bam.files <- unlist(args[3:length(args)]) TSS1000_counts <- paste("TSS1000_",histone_name,"_",cell_type,".txt",sep='') TSS_1000_human <- read.table("human_TSS1000bp.txt",sep="\t",header=F) colnames(TSS_1000_human) <- c("chr","start","end","strand","ensembl_id") TSS_1000_human_granges <- with(TSS_1000_human, GRanges(chr, IRanges(start, end), strand, id=ensembl_id)) param <- readParam(minq=10) frag.len <- 36 reg_counts_1000 <- regionCounts(bam.files,TSS_1000_human_granges,ext=frag.len,param=param) head(assay(reg_counts_1000)) binned <- windowCounts(bam.files, bin=TRUE, width=10000, param=param) normfacs <- normOffsets(binned) test_normalize_counts <- (assay(reg_counts_1000))/normfacs TSS1000_TMM <- asDGEList(reg_counts_1000, norm.factors=normfacs)