Help to identify this error.

[carolinehey]

Hello, I'm new to Rscripts using terminal and shallowHRD. I tried running the old_hg38 script but get these error messages. Any help on how to solve them would be appreciated.

(shallowHRD) root@srvodeappgsq02:/projekt/G150-2023-Functional-Genomics/shallowHRD/shallowHRD# Rscript old_versions/shallowHRD_hg 38_old.R QDNAseq/105025949544_filtered.bam_ratio.txt test cytoband_adapted_hg38.csv Loading required package: ggplot2 Loading required package: stats4 Loading required package: BiocGenerics

Attaching package: ‘BiocGenerics’

The following object is masked from ‘package:gridExtra’:

combine

The following objects are masked from ‘package:stats’:

IQR, mad, sd, var, xtabs

The following objects are masked from ‘package:base’:

anyDuplicated, append, as.data.frame, basename, cbind, colnames,
dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
union, unique, unsplit, which.max, which.min

Loading required package: S4Vectors

Attaching package: ‘S4Vectors’

The following objects are masked from ‘package:base’:

expand.grid, I, unname

Loading required package: IRanges Loading required package: GenomeInfoDb

Sample : 105025949544_filtered Warning message: NaNs produced Warning message: NaNs produced on going... Error in X[i, 5] : subscript out of bounds Error in h(simpleError(msg, call)) : error in evaluating the argument 'x' in selecting a method for function 'diff': argument 'x' must be numeric Calls: localMinima -> diff -> density -> density -> density.default Error in density.default(test) : argument 'x' must be numeric Calls: density -> density -> density.default Don't know how to automatically pick scale for object of type . Defaulting to continuous. Error in geom_vline(): ! Problem while computing aesthetics. ℹ Error occurred in the 2nd layer. Caused by error in FUN(): ! object 'Threshold' not found Backtrace: ▆

1. ├─ggplot2::ggsave(...)
2. │ ├─grid::grid.draw(plot)
3. │ └─ggplot2:::grid.draw.ggplot(plot)
4. │ ├─base::print(x)
5. │ └─ggplot2:::print.ggplot(x)
6. │ ├─ggplot2::ggplot_build(x)
7. │ └─ggplot2:::ggplot_build.ggplot(x)
8. │ └─ggplot2:::by_layer(...)
9. │ ├─rlang::try_fetch(...)
   10. │ │ ├─base::tryCatch(...)
   11. │ │ │ └─base (local) tryCatchList(expr, classes, parentenv, handlers)
   12. │ │ │ └─base (local) tryCatchOne(expr, names, parentenv, handlers[[1L]])
   13. │ │ │ └─base (local) doTryCatch(return(expr), name, parentenv, handler)
   14. │ │ └─base::withCallingHandlers(...)
   15. │ └─ggplot2 (local) f(l = layers[[i]], d = data[[i]])
   16. │ └─l$compute_aesthetics(d, plot)
   17. │ └─ggplot2 (local) compute_aesthetics(..., self = self)
   18. │ └─base::lapply(aesthetics, eval_tidy, data = data, env = env)
   19. │ └─rlang (local) FUN(X[[i]], ...)
   20. └─base::.handleSimpleError(...)
   21. └─rlang (local) h(simpleError(msg, call))
   22. └─handlers[1L]
   23. └─cli::cli_abort(...)
   24. └─rlang::abort(...) Warning message: In cbind(...) : number of rows of result is not a multiple of vector length (arg 1) Error: object 'Threshold' not found Error: object 'THR' not found Error: object 'THR' not found Error in [<-(*tmp*, 1, c_conf, value = 1) : subscript out of bounds Calls: breakSmoothToLGA -> getSegmentID on going... Error in $<-.data.frame(*tmp*, num_line, value = 1:0) : replacement has 2 rows, data has 0 Calls: $<- -> $<-.data.frame Error in geom_point(): ! Problem while computing aesthetics. ℹ Error occurred in the 1st layer. Caused by error in FUN(): ! object 'num_line' not found Backtrace: ▆
10. ├─base::suppressWarnings(...)
11. │ └─base::withCallingHandlers(...)
12. ├─ggplot2::ggsave(...)
13. │ ├─grid::grid.draw(plot)
14. │ └─ggplot2:::grid.draw.ggplot(plot)
15. │ ├─base::print(x)
16. │ └─ggplot2:::print.ggplot(x)
17. │ ├─ggplot2::ggplot_build(x)
18. │ └─ggplot2:::ggplot_build.ggplot(x)
    10. │ └─ggplot2:::by_layer(...)
    11. │ ├─rlang::try_fetch(...)
    12. │ │ ├─base::tryCatch(...)
    13. │ │ │ └─base (local) tryCatchList(expr, classes, parentenv, handlers)
    14. │ │ │ └─base (local) tryCatchOne(expr, names, parentenv, handlers[[1L]])
    15. │ │ │ └─base (local) doTryCatch(return(expr), name, parentenv, handler)
    16. │ │ └─base::withCallingHandlers(...)
    17. │ └─ggplot2 (local) f(l = layers[[i]], d = data[[i]])
    18. │ └─l$compute_aesthetics(d, plot)
    19. │ └─ggplot2 (local) compute_aesthetics(..., self = self)
    20. │ └─ggplot2:::scales_add_defaults(...)
    21. │ └─base::lapply(aesthetics[new_aesthetics], eval_tidy, data = data)
    22. │ └─rlang (local) FUN(X[[i]], ...)
    23. └─base::.handleSimpleError(...)
    24. └─rlang (local) h(simpleError(msg, call))
    25. └─handlers[1L]
    26. └─cli::cli_abort(...)
    27. └─rlang::abort(...) Error in tmp[k, c_chr + 1] : subscript out of bounds Calls: LGA_control Error in tmp[k, c_chr + 1] : subscript out of bounds Calls: LGA_control Error in eval(quote(list(...)), env) : object 'WC' not found Calls: cbind -> standardGeneric -> eval -> eval -> eval Error in names(x) <- value : 'names' attribute [9] must be the same length as the vector [0] Calls: colnames<- -> colnames<- Error in if (test[l, 9] == 0) { : argument is of length zero Warning message: Using size aesthetic for lines was deprecated in ggplot2 3.4.0. ℹ Please use linewidth instead. Error in if (higlight_CCNE1[1, 5] >= 0.5) { : missing value where TRUE/FALSE needed Error: object 'WC' not found

[aeeckhou]

Hello,

First of all, don't use the old script ! The performance is not as good as the latest version and to some extent could introduce wrong calls.

Secondly, I am unsure of how it behaves with the way you wrote the relative paths. Try : "QDNAseq/105025949544_filtered.bam_ratio.txt" -> "./QDNAseq/105025949544_filtered.bam_ratio.txt" "test" -> "./test" "cytoband_adapted_hg38.csv" -> "./cytoband_adapted_hg38.csv"

Thridly, how did you generate you input, with what software ? Please provide an example !

I insist that the latest scripts should really be your primary focus. If you can use : "QDNAseq_from_bam_no_chrX.R" to generate you input "shallowHRD_hg19_1.13_QDNAseq_no_chrX.R" to run the analysis.

Best regards, Alexandre Eeckhoutte

[carolinehey]

Okay, i will try align to hg19 then. The change in relative paths did not solve the problem. I generated the input using your QDNAseq_from_bam_chrX.R script. As far as i see the input looks like your example.

Section of my input: "chromosome" "start" "ratio" "ratio_median" 1 850001 0.174 -0.002 1 900001 -0.144 -0.002 1 950001 0.025 -0.002 1 1000001 0.011 -0.002 1 1050001 -0.135 -0.002 1 1100001 -0.037 -0.002 1 1200001 0.01 -0.002 1 1250001 -0.126 -0.002 1 1300001 0.114 -0.002 1 1350001 0.057 -0.002 1 1450001 -0.148 -0.002 1 1550001 0.006 -0.002 1 1750001 0.068 -0.002 1 1800001 0.059 -0.002 1 1850001 0.037 -0.002

[aeeckhou]

The old script doesn't work with the chrX, it was added after in the lastest versions. "QDNAseq_from_bam_chrX.R" work for hg19 (default of QDNAseq) and a small adaptation should be made for hg38 !

Going with hg19 and "QDNAseq_from_bam_no_chrX.R" followed by "shallowHRD_hg19_1.13_QDNAseq_no_chrX.R" should make the run smooth.

Keep me posted!

Alexandre

[carolinehey]

Using hg19 did solve my problems. Thank you wery much. Another question: I do not align using only chr 1-22 because it will introduce some sort of bias in my opinion. The scripts do work anyway, but do you think it is possible to remove everything but chr 1-22 after alignment instead. Then i do not see it as a problem.

[aeeckhou]

Hello yes it is possible! I'm pretty sure it won't change much regarding the profiles but you can remove everything but the chr 1-22 after alignement.

Best, Alexandre

[aeeckhou]

Hello,

Do you have other questions or can I close this issue ?

Best, Alexandre Eeckhoutte

[carolinehey]

Yes, you can close it. Thank you for the help.

Best, Caroline