<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>

CC BY-SA 4.0
2023-03-16T10:23:15.167Z
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.
terdon
CC BY-SA 4.0
2023-03-16T20:37:52.957Z
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?
Ram RS
CC BY-SA 4.0
2023-03-16T21:07:26.917Z
It will be indeed be x %>% subset. There was a mistake during copy pasting.
Aranyak Goswami
CC BY-SA 4.0
2023-03-17T17:42:11.723Z
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.
Ram RS
CC BY-SA 4.0
2023-03-21T20:43:12.707Z
The file that you include as a bedgraph example doesn't look like [bedgraphs](http://genome.ucsc.edu/goldenPath/help/bedgraph.html).
Maximilian Press
<r><genomics>
Cut and Tag analysis
# 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
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.

**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

