Cut and Tag analysis

Question

I am trying to do Cut and Tag analysis for a set of files and I want to perform a procedure called Tn5 normalization of samples. The basic goal is to generate bedgraph files and does normalization with the weights of a sequencing depth. I have a set of function files that I have defined and then a set of Tn5 normalization scripts. The scripts are both quite long and hence I will be sharing parts of the scripts to get hold of my error and also provide as much information as possible.

My bedgraph files look like this:

replicate species Msample serotype protein            sample      condition    group
1        R1     NHP      02     LK03      NA R1-NHP-02-LK03-NA R1-NHP-LK03-NA LK03-NHP
2        R1     NHP      04     LK03      NA R1-NHP-04-LK03-NA R1-NHP-LK03-NA LK03-NHP
3        R1     NHP      13      KP1      NA R1-NHP-13-KP1-NA R1-NHP-KP1-NA KP1-NHP
4        R1     NHP      42      KP1      NA R1-NHP-42-KP1-NA R1-NHP-KP1-NA KP1-NHP
5        R1     NHP      55     LK03      NA R1-NHP-55-LK03-NA R1-NHP-LK03-NA LK03-NHP
6        R1     NHP      57      KP1      NA R1-NHP-57-KP1-NA R1-NHP-KP1-NA KP1-NHP

I must load a complete set of files called RDS which basically bears a list of the individual bed graph files.

The read lengths for each of the files look like this:

seqnames start  end normcount
73694699     opt2     0   52 0.0000000
73694700     opt2    52   56 0.0243838
73694701     opt2    56   83 0.0731515

The Tn5 R normalization scripts which are relevant

library(dplyr, lib.loc = "/scg/apps/software/r/4.1.2/lib64/R/library")
library(gridExtra, lib.loc = "/scg/apps/software/r/4.1.2/lib64/R/library")
library(Gviz, lib.loc = "/home/ag111/R/x86_64-pc-linux-gnu-library/4.1") #needed to be installed
library(ggplot2)
library(RcppArmadillo, lib.loc = "/scg/apps/software/r/4.1.2/lib64/R/library")
#library(stringi, lib.loc = "/scg/apps/software/r/4.0.3/lib64/R/library")
library(tidyr, lib.loc = "/scg/apps/software/r/4.1.2/lib64/R/library")
library(RColorBrewer, lib.loc = "/scg/apps/software/r/4.1.2/lib64/R/library")
library(ggpubr, lib.loc = "/scg/apps/software/r/4.1.2/lib64/R/library") #needed to be installed
library(GenomicFeatures, lib.loc = "/scg/apps/software/r/4.1.2/lib64/R/library")
library(GenomicRanges, lib.loc = "/scg/apps/software/r/4.1.2/lib64/R/library")

source('AG_CnRT_functions.r')
dir.create(file.path("Tn5exITR_InVivo_norm_full"))

intermediates_dirs <- c(SeqDepth='intermediates/')
#intermediates_dirs <- c(SeqDepth='/labs/markay/Aranyak/Tn5_NHP_rename/Data_analysis/intermediates/')

normalization_method <- 'SeqDepth'

AAV_gsize <- 4410
AAV_anno <- read.delim('AAV_opt2_anno.txt')

#figures_dir <- c(SeqDepth='Tn5exITR_InVivo_norm_full2/')

figures_dir <- c(SeqDepth='Tn5exITR_InVivo_norm_full2/')

############## Tn5 ###########
AAV_normbedgraphs <- readRDS("/labs/markay/Aranyak/Tn5_NHP_rename/Data_analysis/intermediates/AAV_rmdup_SeqDepthNormbedgraphs.rds")

naming_order <- c('replicate','species',"Msample",'serotype','protein') 

AAV_normbedgraphs_exITR <- AAV_normbedgraphs %% lapply(function(x){
  df <- x %% subset(start = 185 & end <=  4167)
  return(df)
})

AAV_normcov_regions_meta_exITR <- list()
AAV_normcov_regions_meta_exITR[[normalization_method]] <- get_normcov_regions_meta_AG( 
  AAV_normbedgraphs=AAV_normbedgraphs_exITR,
  naming_order=naming_order,
  AAV_gsize=AAV_gsize,
  AAV_anno=AAV_anno)


Merging data
AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype[AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype == "LK"] <- "LK03"
AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype[AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype == "iG"] <- "AM"
AAV_normcov_regions_meta_exITR[naming_order] <- lapply(AAV_normcov_regions_meta_exITR[naming_order], factor)
AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype <- ordered(AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype, levels = c("AM", "LK03", "DJ"))

AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype[AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype == "LK"] <- "LK03"
AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype[AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype == "iG"] <- "KP1"
AAV_normcov_regions_meta_exITR[naming_order] <- lapply(AAV_normcov_regions_meta_exITR[naming_order], factor)
AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype <- ordered(AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype, levels = c("KP1", "LK03"))


metadata <- get_metadata_AG3(AAV_normbedgraphs=AAV_normbedgraphs_exITR,
                             #naming_order=naming_order)
metadata <- get_metadata_KP(AAV_normbedgraphs=AAV_normbedgraphs_exITR,
                             naming_order=naming_order)

head(metadata)

#coverage separated by AAV regions and species
#to_plot <- c()
#to_plot <- AAV_normcov_regions_meta_exITR[[normalization_method]] %% gather('region','normcov',2:6) 
#to_plot$region <- factor(to_plot$region, levels=c(as.character(AAV_anno_exITR$region),'AAV_normcov'))
#max_y <- ceiling(max( log2(to_plot$normcov1) ) )
#min_y <- floor( min( log2(to_plot$normcov1) ) ).

After loading of the rds file the error I am getting while I am executing the code block

 AAV_normcov_regions_meta_exITR <- list()
 #AAV_normcov_regions_meta_exITR[[normalization_method]] <- get_normcov_regions_meta_AG3(
 AAV_normcov_regions_meta_exITR[[normalization_method]] <- get_normcov_regions_meta_AG( 
   AAV_normbedgraphs=AAV_normbedgraphs_exITR,
   naming_order=naming_order,
   AAV_gsize=AAV_gsize,
   AAV_anno=AAV_anno)

Error in if (!all(AAV_normcov_regions$sample == metadata$sample)) { : 
  missing value where TRUE/FALSE needed.

The code block which I think for the function file that is valid is

AAV_normcov_regions[is.na(AAV_normcov_regions)] <- 0
  AAV_normcov_regions[is.null(AAV_normcov_regions)] <- 0
  AAV_normcov_regions[is.na(AAV_normcov_regions)] = 0
  AAV_normcov_regions <- AAV_normcov_regions[ match(metadata$sample,AAV_normcov_regions$sample) , ]
  if ( !all(AAV_normcov_regions$sample == metadata$sample) ){ 
    stop('Unsuccessful matching to metadata') 
  }    
  return(AAV_normcov_regions).

This function is working for another set of samples to generate a plot list. The plot list I am having is to_plot

[1] sample      AAV_normcov replicate   species     Msample     serotype    protein     condition   group      
<0 rows> 

(or 0-length row.names) where as what is expected is something like this

  sample   Enhancer        FLuc      polyA   Promoter       WPRE AAV_normcov replicate batch species drug serotype condition group
1  R1-LiD-H-N-DJ 0.42244801 0.290929537 0.43640203 0.22379898 0.34125606  0.26650388        R1   LiD       H    N       DJ      H-DJ  H-DJ
2  R1-LiD-H-N-iG 1.42930414 0.259549496 0.57248120 0.37425687 0.42189340  0.36092640        R1   LiD       H    N       AM      H-iG  H-iG
3  R1-LiD-H-N-LK 0.25754367 0.317786560 0.27373645 0.19470267 0.29818285  0.24753280        R1   LiD       H    N     LK03      H-LK  H-LK
4  R1-LiD-M-N-DJ 0.24578258 0.097169668 0.14489638 0.06470209 0.11091813  0.09308374        R1   LiD       M    N       DJ      M-DJ  M-DJ

I am happy to share other code blocks if it is required to answer the question.