Body: <p>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.</p>
  <p>My bedgraph files look like this:</p>
  <pre><code>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
  </code></pre>
  <p>I must load a complete set of files called RDS which basically bears a list of the individual bed graph files.</p>
  <p>The read lengths for each of the files look like this:</p>
  <pre><code>seqnames start  end normcount
  73694699     opt2     0   52 0.0000000
  73694700     opt2    52   56 0.0243838
  73694701     opt2    56   83 0.0731515
  </code></pre>
  <p>The Tn5 R normalization scripts which are relevant</p>
  <pre class="lang-r prettyprint-override"><code>library(dplyr, lib.loc = &quot;/scg/apps/software/r/4.1.2/lib64/R/library&quot;)
  library(gridExtra, lib.loc = &quot;/scg/apps/software/r/4.1.2/lib64/R/library&quot;)
  library(Gviz, lib.loc = &quot;/home/ag111/R/x86_64-pc-linux-gnu-library/4.1&quot;) #needed to be installed
  library(ggplot2)
  library(RcppArmadillo, lib.loc = &quot;/scg/apps/software/r/4.1.2/lib64/R/library&quot;)
  #library(stringi, lib.loc = &quot;/scg/apps/software/r/4.0.3/lib64/R/library&quot;)
  library(tidyr, lib.loc = &quot;/scg/apps/software/r/4.1.2/lib64/R/library&quot;)
  library(RColorBrewer, lib.loc = &quot;/scg/apps/software/r/4.1.2/lib64/R/library&quot;)
  library(ggpubr, lib.loc = &quot;/scg/apps/software/r/4.1.2/lib64/R/library&quot;) #needed to be installed
  library(GenomicFeatures, lib.loc = &quot;/scg/apps/software/r/4.1.2/lib64/R/library&quot;)
  library(GenomicRanges, lib.loc = &quot;/scg/apps/software/r/4.1.2/lib64/R/library&quot;)
  
  source('AG_CnRT_functions.r')
  dir.create(file.path(&quot;Tn5exITR_InVivo_norm_full&quot;))
  
  intermediates_dirs &lt;- c(SeqDepth='intermediates/')
  #intermediates_dirs &lt;- c(SeqDepth='/labs/markay/Aranyak/Tn5_NHP_rename/Data_analysis/intermediates/')
  
  normalization_method &lt;- 'SeqDepth'
  
  AAV_gsize &lt;- 4410
  AAV_anno &lt;- read.delim('AAV_opt2_anno.txt')
  
  #figures_dir &lt;- c(SeqDepth='Tn5exITR_InVivo_norm_full2/')
  
  figures_dir &lt;- c(SeqDepth='Tn5exITR_InVivo_norm_full2/')
  
  ############## Tn5 ###########
  AAV_normbedgraphs &lt;- readRDS(&quot;/labs/markay/Aranyak/Tn5_NHP_rename/Data_analysis/intermediates/AAV_rmdup_SeqDepthNormbedgraphs.rds&quot;)
  
  naming_order &lt;- c('replicate','species',&quot;Msample&quot;,'serotype','protein') 
  
  AAV_normbedgraphs_exITR &lt;- AAV_normbedgraphs %% lapply(function(x){
    df &lt;- x %% subset(start = 185 &amp; end &lt;=  4167)
    return(df)
  })
  
  AAV_normcov_regions_meta_exITR &lt;- list()
  AAV_normcov_regions_meta_exITR[[normalization_method]] &lt;- 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]]<span class="math-container">$serotype[AAV_normcov_regions_meta_exITR[[normalization_method]]$</span>serotype == &quot;LK&quot;] &lt;- &quot;LK03&quot;
  AAV_normcov_regions_meta_exITR[[normalization_method]]<span class="math-container">$serotype[AAV_normcov_regions_meta_exITR[[normalization_method]]$</span>serotype == &quot;iG&quot;] &lt;- &quot;AM&quot;
  AAV_normcov_regions_meta_exITR[naming_order] &lt;- lapply(AAV_normcov_regions_meta_exITR[naming_order], factor)
  AAV_normcov_regions_meta_exITR[[normalization_method]]<span class="math-container">$serotype &lt;- ordered(AAV_normcov_regions_meta_exITR[[normalization_method]]$</span>serotype, levels = c(&quot;AM&quot;, &quot;LK03&quot;, &quot;DJ&quot;))
  
  AAV_normcov_regions_meta_exITR[[normalization_method]]<span class="math-container">$serotype[AAV_normcov_regions_meta_exITR[[normalization_method]]$</span>serotype == &quot;LK&quot;] &lt;- &quot;LK03&quot;
  AAV_normcov_regions_meta_exITR[[normalization_method]]<span class="math-container">$serotype[AAV_normcov_regions_meta_exITR[[normalization_method]]$</span>serotype == &quot;iG&quot;] &lt;- &quot;KP1&quot;
  AAV_normcov_regions_meta_exITR[naming_order] &lt;- lapply(AAV_normcov_regions_meta_exITR[naming_order], factor)
  AAV_normcov_regions_meta_exITR[[normalization_method]]<span class="math-container">$serotype &lt;- ordered(AAV_normcov_regions_meta_exITR[[normalization_method]]$</span>serotype, levels = c(&quot;KP1&quot;, &quot;LK03&quot;))
  
  
  metadata &lt;- get_metadata_AG3(AAV_normbedgraphs=AAV_normbedgraphs_exITR,
                               #naming_order=naming_order)
  metadata &lt;- get_metadata_KP(AAV_normbedgraphs=AAV_normbedgraphs_exITR,
                               naming_order=naming_order)
  
  head(metadata)
  
  #coverage separated by AAV regions and species
  #to_plot &lt;- c()
  #to_plot &lt;- AAV_normcov_regions_meta_exITR[[normalization_method]] %% gather('region','normcov',2:6) 
  #to_plot<span class="math-container">$region &lt;- factor(to_plot$</span>region, levels=c(as.character(AAV_anno_exITR<span class="math-container">$region),'AAV_normcov'))
  #max_y &lt;- ceiling(max( log2(to_plot$</span>normcov1) ) )
  #min_y &lt;- floor( min( log2(to_plot$normcov1) ) ).
  </code></pre>
  <p>After loading of the rds file the error I am getting while I am executing the code block</p>
  <pre class="lang-r prettyprint-override"><code> AAV_normcov_regions_meta_exITR &lt;- list()
   #AAV_normcov_regions_meta_exITR[[normalization_method]] &lt;- get_normcov_regions_meta_AG3(
   AAV_normcov_regions_meta_exITR[[normalization_method]] &lt;- 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<span class="math-container">$sample == metadata$</span>sample)) { : 
    missing value where TRUE/FALSE needed.
  </code></pre>
  <p>The code block which I think for the function file that is valid  is</p>
  <pre class="lang-r prettyprint-override"><code>AAV_normcov_regions[is.na(AAV_normcov_regions)] &lt;- 0
    AAV_normcov_regions[is.null(AAV_normcov_regions)] &lt;- 0
    AAV_normcov_regions[is.na(AAV_normcov_regions)] = 0
    AAV_normcov_regions &lt;- 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).
  </code></pre>
  <p>This function is working for another set of samples to generate a plot list.
  The plot list I am having is to_plot</p>
  <pre class="lang-r prettyprint-override"><code>[1] sample      AAV_normcov replicate   species     Msample     serotype    protein     condition   group      
  &lt;0 rows&gt; 
  </code></pre>
  <p>(or 0-length row.names) where as what is expected is something like this</p>
  <pre class="lang-r prettyprint-override"><code>  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
  </code></pre>
  <p>I am happy to share other code blocks if it is required to answer the question.</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-03-16T10:23:15.167Z
Comments.Text: Please see our help page on [formatting](https://bioinformatics.stackexchange.com/help/formatting). I tried to format your code as code, but please check and make sure I didn't miss something.
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-03-16T20:37:52.957Z
Comments.Text: Is there a typo in your code? You have a `x %% subset` which should probably be `x %>% subset`. Did the code run or did you edit it after copy-pasting and make a mistake there?
Comments.UserDisplayName: Ram RS
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-03-16T21:07:26.917Z
Comments.Text: It will be indeed be x %>% subset. There was a mistake during copy pasting.
Comments.UserDisplayName: Aranyak Goswami
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-03-17T17:42:11.723Z
Comments.Text: Check your code for more copy-paste errors and make sure everything lines up. Edit your post and fix all the errors you find, including the `%% subset` error. Please use proper code formatting for inline code and code blocks.
Comments.UserDisplayName: Ram RS
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-03-21T20:43:12.707Z
Comments.Text: The file that you include as a bedgraph example doesn't look like [bedgraphs](http://genome.ucsc.edu/goldenPath/help/bedgraph.html).
Comments.UserDisplayName: Maximilian Press
Tags: <r><genomics>
Title: Cut and Tag analysis
ThreadText: # Cut and Tag analysis
  
  (Asked by: username_0 on 2023-03-15)
  
  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
  
  ```r
  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/username_0/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/username_0/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
  
  ```r
   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
  
  ```r
  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
  
  ```r
  [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
  
  ```r
    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.
  
  **Tags:** r, genomics
  
  > **Comments:**
  > **username_1**: Please see our help page on [formatting](https://bioinformatics.stackexchange.com/help/formatting). I tried to format your code as code, but please check and make sure I didn't miss something.  
  > **username_2**: Is there a typo in your code? You have a `x %% subset` which should probably be `x %>% subset`. Did the code run or did you edit it after copy-pasting and make a mistake there?  
  > **username_0**: It will be indeed be x %>% subset. There was a mistake during copy pasting.  
  > **username_2**: Check your code for more copy-paste errors and make sure everything lines up. Edit your post and fix all the errors you find, including the `%% subset` error. Please use proper code formatting for inline code and code blocks.  
  > **username_3**: The file that you include as a bedgraph example doesn't look like [bedgraphs](http://genome.ucsc.edu/goldenPath/help/bedgraph.html).  
  
  ---
  
  ## 0 answers
  
  