SEQUENZA error with .extract function

Karthigayini Sivaprakasam

unread,

Apr 23, 2016, 7:37:14 PM4/23/16

to Sequenza User Group

Hi

I used sequenza to run 6 different samples and it ran with no issues for all except one sample. The error comes while executing the sequenza.extract function and it is:

test <- sequenza.extract(data.file)

Processing 1: 51 variant calls; 3212 heterozygous positions; 396786 homozygous positions.

Processing 2: sed: 1: " 4e+05,670378 p": invalid command code e

gzip: error writing to output: Broken pipe

gzip: XAM6.bam_small2.seqz.gz: uncompress failed

Error in read.table(file = file, header = header, sep = sep, quote = quote, :

no lines available in input

My input file has 74MB of compressed data:

'data.frame': 6761140 obs. of 14 variables:

$ chromosome : chr "1" "1" ...

$ position : int 14011 14061 14062 14086 14119 ...

$ base.ref : chr "N" "N" ...

$ depth.normal : int 15 15 17 17 12 ...

$ depth.tumor : int 11 11 9 9 7 ...

$ depth.ratio : num 0.775 0.775 0.509 0.509 0.607 ...

$ Af : num 1 1 1 1 1 ...

$ Bf : num 0 0 0 0 0 ...

$ zygosity.normal: chr "hom" "hom" ...

$ GC.percent : num 62 62 76 76 72 ...

$ good.reads : num 49 49 25 25 3 ...

$ AB.normal : chr "N" "N" ...

$ AB.tumor : chr "." "." ...

$ tumor.strand : chr "0" "0" ...

Can anyone kindly help me to resolve this error?

Francesco Favero

unread,

Apr 25, 2016, 7:35:02 AM4/25/16

to Karthigayini Sivaprakasam, Sequenza User Group

Hi Karthigayinim,

It looks like the file have some problem…

Which OS and R version are you running?

can you try to run the command:

gc.stats <- gc.sample.stats(file = data.file)

and report back the content of:

gc.stats$file.metrics

Thanks

Francesco

--
You received this message because you are subscribed to the Google Groups "Sequenza User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sequenza-user-g...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Karthigayini Sivaprakasam

unread,

Apr 25, 2016, 1:24:51 PM4/25/16

to Sequenza User Group, ksiv...@asu.edu

Sure and thank you for replying immediately. I am using CENTOS 6.2 (Linux) and R 3.1.2

I had tried to generate the GZ file again and run in MAC OS too, and got the same error.

The result from the metrics File.

>gc.stats$file.metrics

chr n.lines start end

1 1 399998 1 399998

2 2 270379 399999 670377

3 3 197490 670378 867867

4 4 210032 867868 1077899

5 5 345355 1077900 1423254

6 6 304578 1423255 1727832

7 7 250704 1727833 1978536

8 8 204397 1978537 2182933

9 9 382349 2182934 2565282

10 10 229066 2565283 2794348

11 11 166102 2794349 2960450

12 12 198106 2960451 3158556

13 13 142367 3158557 3300923

14 14 126602 3300924 3427525

15 15 175621 3427526 3603146

16 16 154927 3603147 3758073

17 17 207168 3758074 3965241

18 18 210211 3965242 4175452

19 19 79699 4175453 4255151

20 20 300126 4255152 4555277

21 21 142457 4555278 4697734

22 22 86468 4697735 4784202

23 23 121179 4784203 4905381

24 24 178620 4905382 5084001

25 25 131277 5084002 5215278

26 26 170329 5215279 5385607

27 27 170891 5385608 5556498

28 28 143062 5556499 5699560

29 29 62759 5699561 5762319

30 30 160791 5762320 5923110

31 31 74430 5923111 5997540

32 32 75619 5997541 6073159

33 33 80798 6073160 6153957

34 34 87199 6153958 6241156

35 35 54582 6241157 6295738

36 36 86037 6295739 6381775

37 37 90747 6381776 6472522

38 38 69426 6472523 6541948

M M 548 6541949 6542496

X X 189951 6542497 6732447

JH373233.1 JH373233.1 622 6732448 6733069

JH373235.1 JH373235.1 33 6733070 6733102

JH373236.1 JH373236.1 379 6733103 6733481

JH373238.1 JH373238.1 425 6733482 6733906

JH373239.1 JH373239.1 2776 6733907 6736682

JH373241.1 JH373241.1 129 6736683 6736811

JH373242.1 JH373242.1 154 6736812 6736965

JH373243.1 JH373243.1 491 6736966 6737456

JH373245.1 JH373245.1 182 6737457 6737638

JH373247.1 JH373247.1 82 6737639 6737720

JH373248.1 JH373248.1 171 6737721 6737891

JH373249.1 JH373249.1 173 6737892 6738064

JH373250.1 JH373250.1 81 6738065 6738145

JH373252.1 JH373252.1 144 6738146 6738289

JH373253.1 JH373253.1 507 6738290 6738796

JH373254.1 JH373254.1 338 6738797 6739134

JH373255.1 JH373255.1 39 6739135 6739173

JH373256.1 JH373256.1 63 6739174 6739236

JH373258.1 JH373258.1 148 6739237 6739384

JH373261.1 JH373261.1 64 6739385 6739448

JH373264.1 JH373264.1 54 6739449 6739502

JH373265.1 JH373265.1 193 6739503 6739695

JH373267.1 JH373267.1 331 6739696 6740026

JH373268.1 JH373268.1 87 6740027 6740113

JH373272.1 JH373272.1 38 6740114 6740151

JH373273.1 JH373273.1 17 6740152 6740168

JH373278.1 JH373278.1 360 6740169 6740528

AAEX03019090.1 AAEX03019090.1 3 6740529 6740531

JH373280.1 JH373280.1 208 6740532 6740739

JH373282.1 JH373282.1 35 6740740 6740774

JH373283.1 JH373283.1 32 6740775 6740806

JH373284.1 JH373284.1 95 6740807 6740901

JH373285.1 JH373285.1 31 6740902 6740932

JH373286.1 JH373286.1 306 6740933 6741238

JH373287.1 JH373287.1 35 6741239 6741273

JH373288.1 JH373288.1 94 6741274 6741367

JH373290.1 JH373290.1 67 6741368 6741434

JH373291.1 JH373291.1 278 6741435 6741712

JH373292.1 JH373292.1 75 6741713 6741787

JH373293.1 JH373293.1 28 6741788 6741815

JH373295.1 JH373295.1 19 6741816 6741834

JH373296.1 JH373296.1 229 6741835 6742063

JH373297.1 JH373297.1 352 6742064 6742415

JH373298.1 JH373298.1 171 6742416 6742586

JH373299.1 JH373299.1 76 6742587 6742662

JH373300.1 JH373300.1 69 6742663 6742731

AAEX03019240.1 AAEX03019240.1 54 6742732 6742785

JH373301.1 JH373301.1 101 6742786 6742886

JH373304.1 JH373304.1 252 6742887 6743138

JH373306.1 JH373306.1 990 6743139 6744128

JH373307.1 JH373307.1 36 6744129 6744164

AAEX03019307.1 AAEX03019307.1 205 6744165 6744369

JH373308.1 JH373308.1 159 6744370 6744528

JH373311.1 JH373311.1 121 6744529 6744649

AAEX03019339.1 AAEX03019339.1 24 6744650 6744673

JH373315.1 JH373315.1 27 6744674 6744700

JH373316.1 JH373316.1 130 6744701 6744830

JH373318.1 JH373318.1 12 6744831 6744842

JH373319.1 JH373319.1 148 6744843 6744990

JH373321.1 JH373321.1 20 6744991 6745010

JH373322.1 JH373322.1 42 6745011 6745052

JH373324.1 JH373324.1 110 6745053 6745162

JH373327.1 JH373327.1 13 6745163 6745175

JH373329.1 JH373329.1 175 6745176 6745350

JH373330.1 JH373330.1 124 6745351 6745474

AAEX03019491.1 AAEX03019491.1 132 6745475 6745606

JH373331.1 JH373331.1 100 6745607 6745706

JH373332.1 JH373332.1 232 6745707 6745938

JH373333.1 JH373333.1 14 6745939 6745952

JH373334.1 JH373334.1 140 6745953 6746092

JH373335.1 JH373335.1 48 6746093 6746140

JH373336.1 JH373336.1 154 6746141 6746294

JH373338.1 JH373338.1 117 6746295 6746411

JH373339.1 JH373339.1 2 6746412 6746413

JH373340.1 JH373340.1 9 6746414 6746422

JH373341.1 JH373341.1 21 6746423 6746443

JH373342.1 JH373342.1 31 6746444 6746474

JH373343.1 JH373343.1 113 6746475 6746587

JH373344.1 JH373344.1 39 6746588 6746626

JH373345.1 JH373345.1 71 6746627 6746697

JH373347.1 JH373347.1 61 6746698 6746758

JH373348.1 JH373348.1 16 6746759 6746774

JH373349.1 JH373349.1 26 6746775 6746800

JH373350.1 JH373350.1 36 6746801 6746836

JH373351.1 JH373351.1 282 6746837 6747118

AAEX03019681.1 AAEX03019681.1 4 6747119 6747122

AAEX03019691.1 AAEX03019691.1 35 6747123 6747157

JH373356.1 JH373356.1 11 6747158 6747168

JH373359.1 JH373359.1 10 6747169 6747178

JH373362.1 JH373362.1 15 6747179 6747193

JH373366.1 JH373366.1 19 6747194 6747212

JH373368.1 JH373368.1 200 6747213 6747412

JH373372.1 JH373372.1 144 6747413 6747556

JH373373.1 JH373373.1 253 6747557 6747809

JH373375.1 JH373375.1 13 6747810 6747822

JH373376.1 JH373376.1 49 6747823 6747871

JH373377.1 JH373377.1 232 6747872 6748103

JH373378.1 JH373378.1 647 6748104 6748750

JH373380.1 JH373380.1 2 6748751 6748752

JH373381.1 JH373381.1 288 6748753 6749040

JH373382.1 JH373382.1 39 6749041 6749079

JH373383.1 JH373383.1 14 6749080 6749093

AAEX03019893.1 AAEX03019893.1 20 6749094 6749113

JH373385.1 JH373385.1 23 6749114 6749136

JH373386.1 JH373386.1 44 6749137 6749180

JH373387.1 JH373387.1 55 6749181 6749235

JH373389.1 JH373389.1 420 6749236 6749655

JH373390.1 JH373390.1 8 6749656 6749663

JH373392.1 JH373392.1 187 6749664 6749850

JH373394.1 JH373394.1 50 6749851 6749900

JH373395.1 JH373395.1 2 6749901 6749902

AAEX03019987.1 AAEX03019987.1 8 6749903 6749910

JH373397.1 JH373397.1 17 6749911 6749927

JH373398.1 JH373398.1 35 6749928 6749962

JH373400.1 JH373400.1 151 6749963 6750113

JH373401.1 JH373401.1 6 6750114 6750119

JH373402.1 JH373402.1 30 6750120 6750149

AAEX03020054.1 AAEX03020054.1 39 6750150 6750188

AAEX03020055.1 AAEX03020055.1 13 6750189 6750201

JH373403.1 JH373403.1 6 6750202 6750207

AAEX03020060.1 AAEX03020060.1 9 6750208 6750216

JH373404.1 JH373404.1 291 6750217 6750507

JH373405.1 JH373405.1 18 6750508 6750525

JH373407.1 JH373407.1 5 6750526 6750530

JH373408.1 JH373408.1 111 6750531 6750641

JH373409.1 JH373409.1 106 6750642 6750747

AAEX03020105.1 AAEX03020105.1 66 6750748 6750813

JH373412.1 JH373412.1 85 6750814 6750898

AAEX03020144.1 AAEX03020144.1 29 6750899 6750927

JH373414.1 JH373414.1 49 6750928 6750976

JH373415.1 JH373415.1 122 6750977 6751098

AAEX03020161.1 AAEX03020161.1 343 6751099 6751441

JH373416.1 JH373416.1 77 6751442 6751518

JH373418.1 JH373418.1 44 6751519 6751562

JH373419.1 JH373419.1 4 6751563 6751566

JH373420.1 JH373420.1 38 6751567 6751604

JH373421.1 JH373421.1 43 6751605 6751647

JH373422.1 JH373422.1 6 6751648 6751653

JH373423.1 JH373423.1 8 6751654 6751661

JH373424.1 JH373424.1 33 6751662 6751694

AAEX03020242.1 AAEX03020242.1 20 6751695 6751714

JH373425.1 JH373425.1 30 6751715 6751744

AAEX03020251.1 AAEX03020251.1 81 6751745 6751825

JH373426.1 JH373426.1 342 6751826 6752167

JH373427.1 JH373427.1 68 6752168 6752235

JH373428.1 JH373428.1 41 6752236 6752276

JH373429.1 JH373429.1 17 6752277 6752293

JH373430.1 JH373430.1 77 6752294 6752370

JH373432.1 JH373432.1 137 6752371 6752507

JH373433.1 JH373433.1 55 6752508 6752562

JH373434.1 JH373434.1 8 6752563 6752570

JH373435.1 JH373435.1 197 6752571 6752767

JH373436.1 JH373436.1 137 6752768 6752904

JH373437.1 JH373437.1 314 6752905 6753218

JH373438.1 JH373438.1 102 6753219 6753320

JH373439.1 JH373439.1 11 6753321 6753331

JH373440.1 JH373440.1 24 6753332 6753355

JH373442.1 JH373442.1 27 6753356 6753382

JH373443.1 JH373443.1 24 6753383 6753406

AAEX03020403.1 AAEX03020403.1 2 6753407 6753408

JH373445.1 JH373445.1 18 6753409 6753426

JH373446.1 JH373446.1 4 6753427 6753430

AAEX03020427.1 AAEX03020427.1 19 6753431 6753449

JH373449.1 JH373449.1 22 6753450 6753471

JH373450.1 JH373450.1 81 6753472 6753552

JH373451.1 JH373451.1 8 6753553 6753560

JH373452.1 JH373452.1 66 6753561 6753626

JH373453.1 JH373453.1 4 6753627 6753630

AAEX03020482.1 AAEX03020482.1 27 6753631 6753657

JH373457.1 JH373457.1 15 6753658 6753672

AAEX03020497.1 AAEX03020497.1 6 6753673 6753678

JH373461.1 JH373461.1 2 6753679 6753680

AAEX03020524.1 AAEX03020524.1 4 6753681 6753684

AAEX03020525.1 AAEX03020525.1 37 6753685 6753721

JH373465.1 JH373465.1 175 6753722 6753896

AAEX03020533.1 AAEX03020533.1 20 6753897 6753916

JH373467.1 JH373467.1 337 6753917 6754253

AAEX03020568.1 AAEX03020568.1 135 6754254 6754388

JH373470.1 JH373470.1 134 6754389 6754522

JH373473.1 JH373473.1 123 6754523 6754645

JH373475.1 JH373475.1 23 6754646 6754668

JH373476.1 JH373476.1 44 6754669 6754712

JH373477.1 JH373477.1 59 6754713 6754771

JH373478.1 JH373478.1 162 6754772 6754933

AAEX03020671.1 AAEX03020671.1 234 6754934 6755167

JH373482.1 JH373482.1 15 6755168 6755182

JH373485.1 JH373485.1 4568 6755183 6759750

JH373486.1 JH373486.1 26 6759751 6759776

JH373488.1 JH373488.1 79 6759777 6759855

JH373489.1 JH373489.1 7 6759856 6759862

JH373493.1 JH373493.1 82 6759863 6759944

JH373494.1 JH373494.1 59 6759945 6760003

JH373495.1 JH373495.1 66 6760004 6760069

JH373496.1 JH373496.1 19 6760070 6760088

AAEX03020761.1 AAEX03020761.1 13 6760089 6760101

JH373498.1 JH373498.1 4 6760102 6760105

JH373501.1 JH373501.1 2 6760106 6760107

JH373502.1 JH373502.1 8 6760108 6760115

JH373503.1 JH373503.1 13 6760116 6760128

JH373507.1 JH373507.1 9 6760129 6760137

JH373508.1 JH373508.1 8 6760138 6760145

JH373510.1 JH373510.1 9 6760146 6760154

JH373512.1 JH373512.1 6 6760155 6760160

JH373513.1 JH373513.1 11 6760161 6760171

JH373515.1 JH373515.1 533 6760172 6760704

AAEX03020959.1 AAEX03020959.1 132 6760705 6760836

JH373524.1 JH373524.1 14 6760837 6760850

JH373526.1 JH373526.1 22 6760851 6760872

JH373528.1 JH373528.1 6 6760873 6760878

JH373529.1 JH373529.1 11 6760879 6760889

JH373530.1 JH373530.1 251 6760890 6761140

To unsubscribe from this group and stop receiving emails from it, send an email to sequenza-user-group+unsub...@googlegroups.com.

Francesco Favero

unread,

Apr 25, 2016, 4:14:39 PM4/25/16

to Karthigayini Sivaprakasam, Sequenza User Group

Wow a lot of chromosomes there :-).

I see the problem... It uses 4e5 instead of 400000 as line number to load chr 2.
To validate that you can try:

file.lines <- gc.stats$file.metrics[ gc.stats$file.metrics$chr == "2", ]
seqz.data <- read.seqz(file, gz = gz, n.lines = c(file.lines$start, file.lines$end))

And it should give the error...

I never experienced this error. I will need to fix the R code by telling explicitly to avoid exponential notation...

As a quick fix you can try to run the command:

options("scipen"=100, "digits"=4)

(Just found this on stackoverflow, maybe you need to play with the parameters)

Before running sequenza.extract (or you can try to load only chr 2).

I hope this works for a quick fix.

Francesco

To unsubscribe from this group and stop receiving emails from it, send an email to sequenza-user-g...@googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

--
You received this message because you are subscribed to the Google Groups "Sequenza User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email to sequenza-user-g...@googlegroups.com.

Karthigayini Sivaprakasam

unread,

Apr 25, 2016, 7:20:40 PM4/25/16

to Sequenza User Group, ksiv...@asu.edu

It worked the magic. Thank you Favero.

On a general note, is sequenza human (organism) specific?

To unsubscribe from this group and stop receiving emails from it, send an email to sequenza-user-group+unsub...@googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

--
You received this message because you are subscribed to the Google Groups "Sequenza User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email to sequenza-user-group+unsub...@googlegroups.com.

Francesco Favero

unread,

Apr 26, 2016, 7:11:07 AM4/26/16

to Karthigayini Sivaprakasam, Sequenza User Group

Glad to hear it worked! :)

It shouldn’t be, but in practice the copynumber package, the dependency that we use for segmentation, uses the hg19 to get the chromosome arms (p and q). So the assembly argument in sequenza.extract takes care only of the centromere position, but still it could be troublesome fon non-human. Long story short, be aware of the centromeres…

A part from that, there is a chance of failure when using sequenza on non-heterozygous chromosomes (usually chrM and the other non canonical contigs), because no heterozygous position where found.

The gender chromosomes are “configurable” in the arguments of sequenza.fit and sequenza.results (XY = c(X = "X", Y = “Y”))

Francesco

To unsubscribe from this group and stop receiving emails from it, send an email to sequenza-user-g...@googlegroups.com.

Karthigayini Sivaprakasam

unread,

Apr 26, 2016, 1:08:33 PM4/26/16

to Sequenza User Group, ksiv...@asu.edu

So, the sequenza.extract function takes the p and q positions from humans which might differ in other organism. Since the cellularity and ploidy are dependent on CNVs, the results many be skewed. But we would not know the extent of skewness. Is my understanding correct?

I have excluded chromosome M and other contigs in the chromosome list. Hope that it mitigate some of the failures.

Can you please share your thoughts on how could we still use this tool for non-humans?

Francesco Favero

unread,

Apr 27, 2016, 4:34:04 AM4/27/16

to Karthigayini Sivaprakasam, Sequenza User Group

the critical part is only the segmentation, that hods info about the cetromeres.

It is possible to skip the segmentation by providing breaks to the sequenza.extract function.

So in practice you should be able to segment with your favourite algorithm the data, and then use sequenza to quantify tumor content, ploidy and copy numbers.

I’m in the process of using other algorithms as well for the segmentation, I have some working code using bcp and falcon (both from CRAN).

If you want I can share some quick snippets of code (quick experiments mostly). Although I think bcp is more fit to be introduced as possible alternative/replacement seg algo in sequenza, potentially bringing some improvement to the probabilistic model, I have some ‘shareable’ code only for falcon (using just the segmentation algorithm from the package).

Francesco

Karthigayini Sivaprakasam

unread,

Apr 27, 2016, 12:07:24 PM4/27/16

to Sequenza User Group, ksiv...@asu.edu

It makes sense now!

It will be great use if you could share the snippets. Thank you Dr.Favero. This is amazing.

...

Francesco Favero

unread,

Apr 28, 2016, 6:40:24 AM4/28/16

to Sequenza User Group, ksiv...@asu.edu

I've noticed that falcon was recently being updated to version 0.2, so this code works only with the new version:

First things you need to define your centromeres, with a simple dataframe:

This is an example with hg19

chromosome <- c("chr1", "chr2", "chr3", "chr4", "chr5", "chr6", "chr7", "chr8", "chr9", "chr10",
"chr11", "chr12", "chr13", "chr14", "chr15", "chr16", "chr17", "chr18", "chr19",
"chr20", "chr21", "chr22", "chrX", "chrY")
start.pos <- c(121535434, 92326171, 90504854, 49660117, 46405641, 58830166, 58054331,
43838887, 47367679, 39254935, 51644205, 34856694, 16000000, 16000000,
17000000, 35335801, 22263006, 15460898, 24681782, 26369569, 11288129,
13000000, 58632012, 10104553)
end.pos <- c(124535434, 95326171, 93504854, 52660117, 49405641, 61830166, 61054331,
46838887, 50367679, 42254935, 54644205, 37856694, 19000000, 19000000,
20000000, 38335801, 25263006, 18460898, 27681782, 29369569, 14288129,
16000000, 61632012, 13104553)
centromeres <- data.frame(chromosome = chromosome,
start.pos = start.pos,
end.pos = end.pos)

Now I'm using a mix of sequenza and falcon to make a function that returns breakpoints for a given chromosome

falcon.seg.seqz <- function(data.file, chromosome, centromeres){
require(sequenza)
seqz <- read.seqz(data.file, chr.name = chromosome)
chrom <- gsub(chromosome, pattern = "chr", replacement = "")
centromeres$chromosome <- gsub(centromeres$chromosome, pattern = "chr", replacement = "")
seqz <- seqz[seqz$zygosity.normal == "het", ]
get.tauhat <- function(seqz, ...) {
require(falcon)
at <- round(seqz$depth.tumor * seqz$Af, 0)
bt <- round(seqz$depth.tumor * seqz$Bf, 0)
an <- round(seqz$depth.normal * 0.55, 0)
bn <- round(seqz$depth.normal * 0.45, 0)
getChangepoints(data.frame(AT = at, BT = bt, AN = an, BN = bn) , ...)
}
p <- seqz$position < centromeres$start.pos[centromeres$chromosome == chrom]
q <- seqz$position > centromeres$end.pos[centromeres$chromosome == chrom]
pos.p <- seqz$position[p]
pos.q <- seqz$position[q]
l.p <- length(pos.p)
l.q <- length(pos.q)
do.breaks <- function(chrom, tauhat) {
start.pos <- tauhat
start.pos[-1] <- tauhat[-1]+1
data.frame(chrom = chrom,
start.pos = start.pos[-(length(start.pos))],
end.pos = tauhat[-1])
}
chrom <- unique(seqz$chromosome)
if (l.p > 1 & l.q > 1) {
tauhat.p <- get.tauhat(seqz[p, ], verbose = FALSE)
tauhat.p <- c(min(pos.p), pos.p[tauhat.p], max(pos.p))
tauhat.q <- get.tauhat(seqz[q, ], verbose = FALSE)
tauhat.q <- c(min(pos.q), pos.q[tauhat.q], max(pos.q))
breaks.p <- do.breaks(chrom, tauhat.p)
breaks.q <- do.breaks(chrom, tauhat.q)
rbind(breaks.p, breaks.q)
} else if (l.p < 2) {
tauhat.q <- get.tauhat(seqz[q, ], verbose = FALSE)
tauhat.q <- c(min(pos.q), pos.q[tauhat.q], max(pos.q))
do.breaks(chrom, tauhat.q)
} else if (l.q < 2 ) {
tauhat.p <- get.tauhat(seqz[p, ], verbose = FALSE)
tauhat.p <- c(min(pos.p), pos.p[tauhat.p], max(pos.p))
do.breaks(chrom, tauhat.p)
} else {
stop("Segmentation went wrong...")
}
}

that's it all we need, really. To use with sequenza test data:

data.file <- system.file("data", "example.seqz.txt.gz", package = "sequenza")
brk <- list()
chromosome.list <- c(1:22, "X")
#chromosome.list <- paste0("chr", c(1:22, "X")) ## If you have chromosome name like chr1
for (i in chromosome.list){
brk [[i]] <- falcon.seg.seqz(data.file, chromosome = i, centromeres = centromeres)
}
breaks <- do.call(rbind, brk)

Now that we have the breakpoints we can run the usual sequenza pipeline:

test <- sequenza.extract(data.file, breaks = breaks, chromosome.list= chromosome.list)
CP <- sequenza.fit(test, mc.core = 4, segment.filter= 5e6) ## Increase a little bit from the default filter, the segmentation is slightly noisier
sequenza.results(test, CP, out.dir="falcon_seqz_centromeres", sample.id = "falcon_seqz_centromeres")

The segmentation with falcon resulted a little more noisy that with the default sequenza segmentations, but this can probably be improved by tweaking the parameters.

Use this as it is, I think it's possible to do a much better job. I wrote this as a proof of concept, it is not intended for "production" but just for testing purpose.

In a future version of sequenza I hope to introduce new segmentation options that would enable running non-human organism doable out of the box.

Francesco

Karthigayini Sivaprakasam

unread,

Apr 28, 2016, 2:11:06 PM4/28/16

to Sequenza User Group, ksiv...@asu.edu

Thank you. I would try this and get back to you. Between , Is the start and stop positions the starting and ending loci of the centromeres?

Francesco Favero

unread,

Apr 29, 2016, 4:20:13 AM4/29/16

to Karthigayini Sivaprakasam, Sequenza User Group

On 28 Apr 2016, at 20:11, Karthigayini Sivaprakasam <ksiv...@asu.edu> wrote:

Thank you. I would try this and get back to you. Between , Is the start and stop positions the starting and ending loci of the centromeres?

Yes :)

Karthigayini Sivaprakasam

unread,

May 4, 2016, 9:53:59 PM5/4/16

to Sequenza User Group, ksiv...@asu.edu

Here comes the another issue. Do you have a work around (for this work around, I am sorry) if we do not have the centromere intervals?

...

Francesco Favero

unread,

May 6, 2016, 4:32:04 AM5/6/16

to Karthigayini Sivaprakasam, Sequenza User Group

Well, it’s not necessary to have the centromeres coordinates, it is only needed to avoid a segment to span between the two arms of a chromosome.

If you don’t have the centromeres then, I would run the normal pipeline without modification. It will not make a huge difference after all.

Best

Francesco

Karthigayini Sivaprakasam

unread,

May 6, 2016, 1:29:49 PM5/6/16

to Sequenza User Group, ksiv...@asu.edu

Thank you.

To unsubscribe from this group and stop receiving emails from it, send an email to sequenza-user-group+unsub...@googlegroups.com.

Karthigayini Sivaprakasam

unread,

May 10, 2016, 3:01:21 PM5/10/16

to Sequenza User Group, ksiv...@asu.edu

Hi

I tried this code with the centromere intervals for both whole genome and exome and got exactly the same result as without the centromeres (the usual sequenza procedure). Is there any other factor that is playing a role here?

Francesco Favero

unread,

May 11, 2016, 4:36:58 AM5/11/16

to Karthigayini Sivaprakasam, Sequenza User Group

Hi,

What is exactly wrong with the results?

In theory nothing completely wrong happens if the centromeres are not correct, just you will have a breaks where the centromere is suppose to be in the hg19 genome, and you might not have a breakpoint where the actual cent. is in the specific organism.

So the difference might lay only in the breakpoints coordinates.

Francesco

Karthigayini Sivaprakasam

unread,

May 11, 2016, 1:37:00 PM5/11/16

to Sequenza User Group, ksiv...@asu.edu

Hi

Nothing is wrong with the results.

As you mentioned, one should except loss of data in the regions where centromere is not used, if the samples are from non-humans. I assumption is that it will affect the ploidy calling and thereby the cellularity would be little skewed (different) when the centromeric regions are mentioned.

To unsubscribe from this group and stop receiving emails from it, send an email to sequenza-user-group+unsub...@googlegroups.com.

Reply all

Reply to author

Forward