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 <data.frame>. 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]](cnd)
    

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: $&lt;- -&gt; $&lt;-.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:
▆

1. ├─base::suppressWarnings(...)
2. │ └─base::withCallingHandlers(...)
3. ├─ggplot2::ggsave(...)
4. │ ├─grid::grid.draw(plot)
5. │ └─ggplot2:::grid.draw.ggplot(plot)
6. │ ├─base::print(x)
7. │ └─ggplot2:::print.ggplot(x)
8. │ ├─ggplot2::ggplot_build(x)
9. │ └─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]](cnd)
    

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