2
$\begingroup$

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.

$\endgroup$
5
  • $\begingroup$ Please see our help page on formatting. I tried to format your code as code, but please check and make sure I didn't miss something. $\endgroup$
    – terdon
    Mar 16, 2023 at 10:23
  • $\begingroup$ 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? $\endgroup$
    – Ram RS
    Mar 16, 2023 at 20:37
  • $\begingroup$ It will be indeed be x %>% subset. There was a mistake during copy pasting. $\endgroup$ Mar 16, 2023 at 21:07
  • $\begingroup$ 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. $\endgroup$
    – Ram RS
    Mar 17, 2023 at 17:42
  • $\begingroup$ The file that you include as a bedgraph example doesn't look like bedgraphs. $\endgroup$ Mar 21, 2023 at 20:43

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.