text:
# Find two integers from an unordered list to sum to the input

Find 2 integers from an unordered list that sum to a given integer.

1. Given an unordered list of integers, find 2 integers that sum to a given value, print these 2 integers, and indicate success (exit 0). They don't need to be any particular numbers (i.e. the first 2 integers summing to the right number), any pair that sums to the value will work.
2. an integer is positive and greater than zero.
3. a list of integers can be in any data structure including a file of integers - one integer per line.
4. if no integers can be found, indicate a failure (exit 1).
5. two integers at different positions in list must be returned. (i.e. you can't return the same number from the same position twice)

(Note: in the video, these are not exactly the requirements. The 'interviewer' changed his multiple times.)

eg.

sum2 8 <<EOF
1
7
4
6
5
3
8
2
EOF

Prints 3 and 5 and exit status is 0. Note that in this 1,7 and 2,6 would also be allowed results.

sum2 8 <<EOF
1
2
3
4

Returns exit status 1 since no possible combo. 4,4 isn't allowed, per rule 5.

• This could have been a great question if it had had a chance to shake out some of the loose ends in the Sandbox first. For instance, for something like this I would expect to write a function that returned either a falsy value or a pair of numbers. – Neil Mar 20 '17 at 11:55
• In the example, why is the returned pair is (3,5) and not (1,7)? – Rod Mar 20 '17 at 12:09
• How can there be a "first" pair in an unordered list? That's inherently self-contradictory. – Peter Taylor Mar 20 '17 at 12:39
• I don't really think the exit 0/exit 1 thing is a good idea. Many languages can't exist easily like that, and it's generally allowed to exit with an error (i.e. ignore STDERR) Many golfing languages don't even have an easy way to exit by exit code I think – Rɪᴋᴇʀ Mar 20 '17 at 14:22
• On second thought, some answers have gone through some effort to produce exit code 1, so it may be better not to change the requirements now – Luis Mendo Mar 20 '17 at 18:06

# Bash, 84 bytes

My implementation of (roughly) Google's engineer's solution but using bash and an input stream - not my solution, so this does not count.

while read V;do((V<$1))&&{ ((T=R[V]))&&echo$T $V&&exit;((R[$1-V]=V));};done;exit 1

Method

while we can read integer V from input stream if less than target $1 then if already seen$1-V then print $1-V and V and exit 0 (else) save candidate for input$1-V exit 1

# Brachylog, 9 bytes

h⊇Ċ.+~t?∧

Try it online!

Assuming I understood the challenge correctly…

### Explanation

h⊇Ċ          Ċ ('couple') has two elements, and is a subset of the head of the input
Ċ.         Output = Ċ
.+~t?     The sum of the elements of the Output is the tail of the Input
∧    (disable implicit unification)

# Perl 6, 59 bytes

$_=get;put lines().combinations(2).first(*.sum==$_)//exit 1

## Expanded:

$_ = get; # get one line (the value to sum to) put # print with trailing newline lines() # get the rest of the lines of input .combinations(2) # get the possible combinations .first( # find the first one *.sum ==$_    # that sums to the input
)
//                 # if there is no value (｢Nil｣)
exit 1           # exit with a non-zero value (｢put｣ is not executed)

# JavaScript ES6, 58 70 68 64 bytes

a=>b=>{for(i in a)if(a.includes(b-a[i],i+1))return[a[i],b-a[i]]}

Returns a pair of numbers in the form of an array if found, otherwise returns undefined, a falsy value.

f=a=>b=>{for(i in a)if(a.includes(b-a[i],i+1))return[a[i],b-a[i]]}

console.log(f([1,7,4,6,5,3,8,2])(8));
console.log(f([1,2,3,4,5,6,7,8])(8));
console.log(f([1,2,3,4])(8));
console.log(f([2,2])(4));

• The example was 3, 5 but this outputs 1, 7... – Neil Mar 20 '17 at 12:06
• @Neil, sorry, i have amended rules because I messed up. 1,7 is ok. – philcolbourn Mar 20 '17 at 12:57
• It won't work for f([2,2] 4) ? – cliffroot Mar 20 '17 at 14:32
• @cliffroot should work for that case now – Tom Mar 20 '17 at 15:51
• Nice includes trick. – Neil Mar 20 '17 at 18:08

## JavaScript (ES6), 6157 56 bytes

Takes the array of integers a and the expected sum s in currying syntax (a)(s). Returns a pair of matching integers as an array, or undefined if no such pair exists.

a=>s=>(r=a.find((b,i)=>a.some(c=>i--&&b+c==s)))&&[r,s-r]

### Formatted and commented

a =>                      // given an array of integers (a)
s => (                  // and an expected sum (s)
r = a.find((b, i) =>  // look for b at position i in a such that:
a.some(c =>         //   there exists another c in a:
i-- &&            //     - at a different position
b + c == s        //     - satisfying b + c == s
)                   //   end of some()
)                     // end of find(): assign the result to r
) &&                    // if it's not falsy:
[r, s - r]              // return the pair of integers

### Test

let f =

a=>s=>(r=a.find((b,i)=>a.some(c=>i--&&b+c==s)))&&[r,s-r]

console.log(f([1,7,4,6,5,3,8,2])(8))
console.log(f([1,2,3,4])(8))

# Jelly, 14 bytes

ŒcS=⁹$$ÐfḢṄo⁶H  Try it online! This is a function (not a full program) that outputs to standard output. (The TIO link has a wrapper that runs a function and disregards its return value.) This program could be 4 bytes shorter if not for the exit code requirement; returning an exit code of 1 in Jelly is fairly hard. (It's possible that there's a terser way to do this that I've missed.) ## Explanation ŒcS=⁹$$ÐfḢṄo⁶H
Œc                All pairs of values from {the first argument}
Ðf         Take only those which
S=⁹               sum to {the second argument}
           Parse the preceding three builtins as a group

Ṅ       Output this result (plus a newline) on standard output
o⁶     If this value is falsey, replace it with a space character
H    Halve every element of the value

We can halve every integer in a pair just fine, so the o⁶H will do nothing if we found a result, other than returning a useless return value that isn't relevant anyway (the Ṅ serves as a convenient single-byte method to determine the function's return value early, under PPCG rules). However, if we didn't find a result, we end up trying to halve a space character, an operation so meaningless it causes the Jelly interpreter to crash. Fortunately, this crash produces an exit code of 1.

# Perl 5, 51 bytes

46 bytes of code + for 5 bytes for -pli flags.

$\="$_ $v"if$h{$v=$^I-$_};$h{$_}=1}{$\||exit 1

Try it online!

The idea is to iterate on the input list: on a number x ($_), if we previously saw n-x ($^I-$_) then we found what we were looking for, and set $\ to these two values ("$_$v"). At the end, if $\ isn't set, then we exit 1, else it will be implicitly printed. • Does a literal tab work in place of the two characters ^I? – user62131 Mar 20 '17 at 21:58 • @ais523 It looks like I can't. Maybe it was possible on older versions of Perl though. – Dada Mar 21 '17 at 7:57 # Röda, 60 56 bytes f s,a{seq 1,s|{|x|[[x,s-x]]if[x in a,s-x in a-x]}_|pull}  Try it online! This code throws an error if there is no answer. It generates all possible pairs that can form the sum s, ie. 1, s-1, 2, s-2, 3, s-3, ... Then it checks if both numbers are in the array a and if so, pushes them to the stream. pull reads one value from the stream and returns it. If there are no values in the stream, it throws an error. a-x returns the array a with x removed. ## Python 2, 60 bytes This short, until rules with exiting with code 1 are clarified. Now exits with error if nothing is found. -5 bytes thanks to @Peilonrayz -4 bytes thanks to @Rod Try it online a,s=input() while a: x=a.pop() if s-x in a:r=s-x,x print r  • @Peilonrayz Didn't realize that, thank you! – Dead Possum Mar 20 '17 at 15:54 • @Peilonrayz This would violate fith rule: two integers at different positions in list must be returned. (i.e. you can't return the same number from the same position twice) – Dead Possum Mar 20 '17 at 15:59 • You can use spaces+tabs for mixed indentation to reduce 2 bytes or switch to input() to reduce 4 bytes – Rod Mar 20 '17 at 18:54 • @Rod Thanks! Input seems nicer – Dead Possum Mar 20 '17 at 19:50 • @Eric Duminil Yeah. It's equivalent to eval(raw_input()) (I think). – Yytsi Mar 20 '17 at 22:58 # C++ 133 bytes (compiled with clang 4 and gcc 5.3 -std=c++14) #include <set> auto f=[](auto s,int v,int&a,int&b){std::set<int>p;for(auto i:s)if(p.find(i)==end(p))p.insert(v-i);else{a=v-i;b=i;}};  # C 108 bytes void f(int*s,int*e,int v,int*a,int*b){do{int*n=s+1;do if(v-*s==*n){*a=*s;*b=*n;}while(++n<e);}while(++s<e);}  • Welcome to the site! Unfortunately, I think you need to add 15 bytes for #include <set> and a few more for std::set. Although you can also save some bytes if you remove the braces around p.insert(v-i); – DJMcMayhem Mar 20 '17 at 17:39 • @DJMcMayhem oh, thank you. So should i include main()? – em2er Mar 20 '17 at 17:56 • @em2er No, you don't need to include main. We consider (unless stated otherwise in the challenge) that a function is a valid submission. (welcome on the site btw!) – Dada Mar 20 '17 at 17:56 • No, a function submission is fine. (And much shorter because you can take the input as arguments) You just need to count any includes that your submission requires. – DJMcMayhem Mar 20 '17 at 18:00 • @DJMcMayhem @Dada thanks a lot! i am not sure with end too, but it compiles on gcc without std:: (and set if course not) – em2er Mar 20 '17 at 18:05 # Haskell, 34 bytes (n:v)#s|elem(s-n)v=(n,s-n)|1<2=v#s  Try it online! For each element of the list, this function checks if (sum-element) is in the following part of the list. Returns the first couple it finds. If the function reaches the end of the list it throws a "non -exhaustive patterns" error and exits with code 1. • I'm afraid this approach does not work for inputs like [2,2]#4. – Laikoni Mar 20 '17 at 22:38 • @Laikoni Thank you, I hadn't read the challenge well enough. This new version should be correct (and shorter^^) – Leo Mar 21 '17 at 0:12 ## PowerShell, 109 97 bytes param($i,$a)($c=0..($a.count-1))|%{$c-ne($f=$_)|%{if($a[$f]+$a[$_]-eq$i){$a[$f,$_];exit}}};exit 1

Took a 12 byte deal that AdmBorkBork offered

Explanation

# Get the parameter passed where $i is the addition target from the array of numbers in$a
param($i,$a)

($c=0..($a.count-1))|%{
# We are going to have two loops to process the array elements.
# The first loop element will be held by $f$f=$_ # Create a second loop that will be the same as the first except for the position of$f to
# prevent counting the same number twice.
$c|?{$_-ne$f}|%{ # Check if the number at the current array indexes add to the target value. If so print and exit. if($a[$f]+$a[$_]-eq$i){$a[$f],$a[$_];exit}
}

}
# If nothing was found in the loop then we just exit with error.
exit 1

The current rules look for exit code which this does. Those could be removed and just check for numbers being returned and a falsy.

Sample Usage

If the code above was saved as function s

s 8 @(1,2,3,4)
s 8 @(1,7,4,6,5,3,8,2)

• You can save a few more bytes by eliminating $c and looping downward -- ($a.count-1)..1|%{$f=$_;--$_..0|%{if... – AdmBorkBork Mar 21 '17 at 20:12 ## R, 49 bytes function(x,y){r=combn(x,2);r[,colSums(r)==y][,1]}  This finds all 2-combinations of x and returns a matrix. Then, sums by column and finds all the sums that equal to y (so without the [,1] part at the end it will print all the combinations that their sums equal to y) # Japt, 9 bytes Saved many bytes thanks to @ETHproductions à2 æ_x ¥V  Try it online! ### Explanation à2 æ_x ¥V à2 // Creates all combinations of the input, length 2 æ // Returns the first item where: _x // The sum of the two items in each set ¥V // == Second input  ### Example Input: [1,2,3], 4 à2 // [[1,2],[1,3],[2,3]] æ_x // [3, 4, 5 ] ¥V // 3!=4, 4==4 ✓ Output: // 1,3  ## Javascript, 1149686 84 bytes a=>b=>{c=b.length;for(x=0;x<c;x++)for( y=x;++y<c;)if(b[x]+b[y]==a)return[b[x],b[y]]}  Saved 1 byte thanks to @Cyoce and another 8 bytes thanks to @ETHProductions This returns a tuple with the first combination of list-elements that sum up to the given input, or nothing for no match. I've removed the vars in the function; REPL.it crashes without them, but the Chrome Dev Console handles this just fine... Try it online! • Doesn't exit code 1, as the challenge specifically requests for invalid input. It's an invalid answer for now, but I've asked OP about this requirement for languages that can't do that easily. – Rɪᴋᴇʀ Mar 20 '17 at 14:26 • @Matt Yes, that rule is observed: y=x+1 takes care of that. – steenbergh Mar 20 '17 at 15:29 • You can use a=>b=>... to save a byte – Cyoce Mar 20 '17 at 19:18 • You can save another three bytes with for(y=x;++y<b.length;){. Also, you can remove all sets of braces except the outermost one, and you can remove the space after return – ETHproductions Mar 23 '17 at 12:39 ## Clojure, 77 bytes #(first(mapcat(fn[i a](for[b(drop(inc i)%):when(=(+ a b)%2)][a b]))(range)%))  Returns the first such pair or nil. ## Haskell, 62 bytes r=return;s#[]=r 1;s#(a:b)|elem(s-a)b=print(a,s-a)>>r 0|1<2=s#b  I still don't know what's allowed by the challenge and what not. I'm going for a function that prints a pair of numbers and returns 0 if there's a solution and prints nothing and returns 1 if there's no solution. As printing is I/O, I have to lift the return values into the IO-Monad (via return) and the actual type of the function is Num a => IO a. Usage example (with return value printed by the repl): *Main> 4 # [2,2] (2,2) 0  If raising exceptions is allowed, fail will save some bytes (total 51): s#[]=fail"";s#(a:b)|elem(s-a)b=print(a,s-a)|1<2=s#b  # Jelly, 9 bytes ŒcS=¥ÐfḢZ  Jelly has no way of setting the exit code to arbitrary values, so this produces a TypeError for input without a valid solution that will cause the parent interpreter to exit with exit code 1. Try it online! ### How it works ŒcS=¥ÐfḢZ Main link. Argument: A (array of integers), n (integer) Œc Yield all 2-combinations of different elements of A. Ðf Filter by the link to the left. ¥ Combine the two links to the left into a dyadic chain. S Take the sum of the pair. = Compare the result with n. Ḣ Head; extract the first pair of the resulting array. This yields 0 if the array is empty. Z Zip/transpose the result. This doesn't (visibly) alter pairs, but it raise a TypeError for 0.  # Nova, 101 bytes q(Int[] a,Int x)=>a{if(Int y=a.firstWhere({a.contains(x-a.remove(0))}))return [y,x-y];System.exit(1)}  One nice thing about code golf is that it helps me find bugs in my language. e.g. the space required between return and [y,x-y]. Once I add push/pop functions to Array.nova and fix return, would be 96 bytes: q(Int[] a,Int x)=>a{if(Int y=a.firstWhere({a.contains(x-a.pop())}))return[y,x-y];System.exit(1)}  Usage: class Test { static q(Int[] a,Int x)=>a{if(Int y=a.firstWhere({a.contains(x-a.remove(0))}))return [y,x-y];System.exit(1)} public static main(String[] args) { Console.log(q([1, 2, 3, 4, 5], 8)) // [5, 3] Console.log(q([1, 2, 3, 4, 5], 5)) // [1, 4] Console.log(q([1, 2, 3, 4], 8)) // exit code 1 } }  Edit: Also, there's this way at 73 bytes (69 using pop), too: q(Int[] a,Int x)=>[Int y=a.firstOrThrow({a.contains(x-a.remove(0))}),x-y]  firstOrThrow will throw an Exception, which will be uncaught and therefore ultimately exiting the program with exit code 1. ;) This way seems more readable too. # Pyth, 12 bytes hfqsThQ.ceQ2  ### Explanation  .ceQ2 Get all pairs from the second input fqsThQ Find the ones whose sum is the first input h Take the first (exits with error code 1 if there aren't any)  # PHP, 88 bytes for($i=1;$a=$argv[$k=++$i];)for(;$b=$argv[++$k];)if($a+$b==$argv[1])die("$a$b");die(1);

takes input from command line arguments, sum first. Run with -nr.

Fortunately, die/ exit exits with 0 when you give it a string as parameter.

I tried to merge the loops to one; but it requires a longer initialization and test this time.

## Mathematica, 76 bytes

f::e="1";If[(l=Cases[#~Subsets~{2},x_/;Tr@x==#2])=={},Message@f::e,First@l]&

Fairly straightforward: #~Subsets~{2} gets all 2-element subsets of the list, then Cases[...,x_/;Tr@x==#2] picks only the ones whose sum is the number we want. If there are none of these, If[l=={}, Message@f::e,First@l] prints the error message f::e : 1 that we defined earlier (since I have no idea what else "exit status 1" might mean for Mathematica); otherwise, it returns the first entry in the list of pairs that sum to the correct thing.

If we're allowed to return a falsey value instead of doing that weird exit status thing, the following code has 58 bytes:

If[(l=Cases[#~Subsets~{2},x_/;Tr@x==#2])=={},1<0,First@l]&

# Scala, 55 41 bytes

(l,n)=>l combinations 2 find(_.sum==n)get

Returns a list of the two numbers if they exist and throws an error otherwise. Uncaught, this error will result in an exit status of 1.

# Ruby, 53 48 bytes

->a,s{p(a.combination(2).find{|x,y|x+y==s})?0:1}

Input: a is the list, s is the expected sum.

If the 2 numbers are found, print them and return 0, otherwise return 1, as in the specification.

# TI-Basic, 59 bytes

Prompt L1
Prompt X
While 1
L1(1→B
seq(L1(C),C,2,dim(L1→L1
If sum(not(X-L1-B
Then
Disp B,X-B
Return
End
End

Explanation:

Prompt L1               # 4 bytes, input array like "{1, 2, 3}"
Prompt X                # 3 bytes, Input target sum
While 1                 # 3 bytes, until the list is empty
L1(1→B                  # 7 bytes, try the first element (now B)
seq(L1(C),C,2,dim(L1→L1  # 18 bytes, remove first element from list
If sum(not(X-L1-B       # 10 bytes, if any element in the list plus B is the target
Then                    # 2 bytes, then...
Disp B,X-B              # 7 bytes, print it and it's "complement"
Return                  # 2 bytes, and exit gracefully
End                     # 2 bytes
End                     # 1 byte

If the program did not exit gracefully, it will cause an error when there are not enough elements in the list for it to continue.

## CJam, 23 bytes

l~_,1>{e!2f<::+#)}{;;}?

Input is sum numbers. For example: 6 [3 2 3]. Leaves a positive number for truthy and an empty string or 0 for falsey.

Explanation:

l~    e# Read input and evaluate:  | 7 [3 2 3]
_     e# Duplicate:                | 7 [3 2 3] [3 2 3]
,     e# Take the length:          | 7 [3 2 3] 3
1>{   e# If more than 1:           | 7 [3 2 3]
e!  e#   Unique permutations:    | 7 [[2 3 3] [3 2 3] [3 3 2]]
2f< e#   Slice each to length 2: | 7 [[2 3] [3 2] [3 3]]
::+ e#   Some each:              | 7 [5 5 6]
#   e#   Index:                  | -1
)   e#   Increment:              | 0
}{    e# Else:                     | 7 [3 2 3]
;   e#   Pop                     | 7
;   e#   pop                     |
}?    e# Endif
e# Implicit output: 0


text:
How to globally determine (Seasonal) ARIMA Parameters

I have 15-minutely data (96 values per day) over several years for around 340 entities (i.e. 340 data sets or long ts).

Now my task is to forecast a 4-hour window (i.e. 16 observations) for each day in the test set (the last 30% of each entity), for which I use ARIMA models.

In my opinion, it does not make sense to take the whole ts for each entity into account for forecasting the 4-hour window. Because of that, my approach is to take a much shorter period on which I create the Model, for example the last week before the 4-hour window (i.e. 96*7 = 672 Last values).

However, because auto.arima in R takes around 2 minutes to estimate the parameters of the model, I want to estimate the parameters for each entity globally, Meaning: If I have a very long time series (e.g. 105216 data points for an entity), which is stationary and seasonal with a frequency of 96, how would be the approach to determine the Parameters for an ARIMA models which fits best to each 672-data-point-sub-ts (or, more intuitive, each week)?

my initial idea was to loop through each week of the test data, create one model for each week, store the results and compare them for each entity. however, as stated above this takes a lot of time. Is there a better approach?

• Due to the high frequency of your seasonality and the multiple seasonality you might be better of with "tbats" than with "(S)ARIMA" – Ferdi Oct 31 '18 at 8:50
• I want to benchmark (S)ARIMA compared to other models, thats why I'm interested in the parameters. Once I have solid parameters, a model with 672 data points is created in less then a second. So I'm looking for a scientific approach to determine those, as opposed to just "guessing" – bk_ Oct 31 '18 at 9:01
• You might want to look at stats.stackexchange.com/questions/254466/… as it deals with exogenous factors like day-of-the-week , holidays etc while incorporating daily aggregates and arima memory – IrishStat Oct 31 '18 at 9:58
• thanks for the hint, unfortunately it does no help me identify (global) ARIMA parameters... – bk_ Oct 31 '18 at 10:40
• Rob Hyndman has once answered a question here on how to fit a model for multiple time series where the parameter values are identical across the time series. His solution might not work for your dataset that is so large (over $10^5$ observations), but you might try looking the thread up nevertheless. – Richard Hardy Oct 31 '18 at 14:40

I found an answer to the problem.

for each of the 340 entities, I estimated an initial model using auto.arima, based on the last few values. This model is used for the first forecast of the 16-observation forecast. For the next forecasting period (i.e. on the next day), I use my pre trained model, and pass it to the Arima function, in that way the parameters for each entity do not need to be re-estimated. It seems like the models perform pretty well. Furthermore, the Arima update is quite fast.

some pseudo code:

# day 1 fit <- auto.arima(train, seasonal = TRUE)

train consists of the last 4 days, since data is 15-minutely, 96 values per day. Hence, train is a ts object with frequency = 96.

# all following days fit <- Arima(train2, model = fit)

Where train2 are the data points in between.

format: plain


text:
Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles and JavaScript.

# A comprehensive assessment of somatic mutation detection in cancer using whole-genome sequencing

## Abstract

As whole-genome sequencing for cancer genome analysis becomes a clinical tool, a full understanding of the variables affecting sequencing analysis output is required. Here using tumour-normal sample pairs from two different types of cancer, chronic lymphocytic leukaemia and medulloblastoma, we conduct a benchmarking exercise within the context of the International Cancer Genome Consortium. We compare sequencing methods, analysis pipelines and validation methods. We show that using PCR-free methods and increasing sequencing depth to 100 × shows benefits, as long as the tumour:control coverage ratio remains balanced. We observe widely varying mutation call rates and low concordance among analysis pipelines, reflecting the artefact-prone nature of the raw data and lack of standards for dealing with the artefacts. However, we show that, using the benchmark mutation set we have created, many issues are in fact easy to remedy and have an immediate positive impact on mutation detection accuracy.

## Introduction

The International Cancer Genome Consortium (ICGC) is characterizing over 25,000 cancer cases from many forms of cancer1. Currently, there are 74 projects supported by different national and international funding agencies. As innovation and development of sequencing technologies have driven prices down and throughput up, projects have been transitioning from exome to whole-genome sequencing (WGS) of tumour and matched germline samples, supplemented by transcript and methylation analyses when possible, facilitating the discovery of new biology for many different forms of cancer2,3,4,5,6,7,8,9,10. However, as data from the different projects began to be collected and centralized (https://dcc.icgc.org/), it became apparent that there are marked differences in how teams generate WGS data and analyse it. On the basis of cost, capacity and analytical experience, it was initially determined that comprehensive identification of tumour-specific somatic mutations requires WGS with a minimum of 30 × sequence coverage of each the tumour and normal genomes11 with paired reads on the order of 100–250 bp in length, depending on the platform. However, from project to project the sample preparation, coverage of tumour and normal samples and read lengths vary. Even more variability exists in the approaches to identify differences between tumour and normal genomes, evidenced by the many strategies developed to identify somatic single-base mutations (SSM)12, somatic insertion/deletion mutations (SIM) and larger structural changes (rearrangements and chromosome segment copy number changes)5.

This variation makes comparison of mutation calls across cancers challenging because of the unknown contributions of individual pipeline components and parameters on the accuracy of the calls. Benchmark data sets and analytical tools13,14,15,16,17 have been developed for variant calling on normal genomes, while those for cancer have largely focused on SSM detection from exome sequencing12,18. Benchmarking of mutation calling from exome data from The Cancer Genome Atlas has raised concerns about biased inferences and highlights the need for benchmark data sets19. In our study we set out to investigate the factors that need to be considered to generate high-quality whole-genome sequence data and high-confidence variant calls from tumour-normal pairs, including new sources of bias and pitfalls not encountered in exome data.

We explored several benchmarking strategies. First, we evaluated somatic mutation calling pipelines using a common set of 40 × WGS reads of average quality corresponding to a case of chronic lymphocytic leukaemia (CLL). In a second benchmark, we evaluated both sequencing methods and somatic mutation calling pipelines using matched samples from a case of medulloblastoma (MB, a malignant pediatric brain tumour arising in the cerebellum20,21) from the ICGC PedBrain Tumor project. Both cancers exhibit a high degree of tumour purity (95–98%). For each case, we made available unaligned sequence reads of a tumour (40 × genome coverage) and its corresponding normal genome (30 × coverage) to members of the ICGC consortium, who then returned somatic mutation calls. In contrast to the approach taken in a recent benchmark of SSM calling using three simulated tumour genomes22, we have used the sequence from a real tumour-normal pair and made a concerted effort to manually curate both SSMs and SIMs detectable at a sequencing depth 8–10 times in excess of the standard amount (300 ×). We argue that real, not simulated, mutations are more useful for dissecting performance of mutation callers with respect to real genome-wide mutational signatures, and methods for detecting insertion–deletion mutations, an even bigger challenge to somatic mutation callers, must also be benchmarked. Our study has two main results: one, we identify outstanding issues in somatic mutation analysis from WGS data and begin to formulate a set of best practices to be adopted more widely by genome researchers, and, two, we provide two benchmark data sets for testing or developing new somatic mutation calling pipelines.

## Results

### WGS data generation

We conducted a first benchmark exercise using WGS data generated from a CLL tumour-normal pair and then a second using a case of MB. Both tumour types were expected to have relatively low mutational load and not very pervasive structural changes: CLL has a few known translocations and large copy number variants, and MB exhibits a large degree of tetraploidy but is otherwise typically free of complex rearrangements. The quality of the CLL data was below today’s standards but typical for the time it was produced, while the quality of the MB library preparation and produced sequence were of high quality. The validation strategies also differed. For CLL we chose to validate submitted mutations by target capture and sequencing with two independent platforms (MiSeq and IonTorrent). This approach was limited by technical issues inherent to target capture and sequencing on these other platforms, which led to a low rate of independently validated mutations. Moreover, the real false-negative (FN) rates were underestimated because of the limited coverage provided to participants. For the MB data set, 300 × in sequence reads were generated by five different sequencing centres, which we joined and used to create a curated set of mutations. The MB tumour presented with a tetraploid background combined with other changes of ploidy in chromosomes 1, 4, 8, 15 and 17 (Supplementary Fig. 1), giving us the opportunity to benchmark performance at lower mutant allele frequencies. Moreover, the high depth and relatively unbiased coverage of the genome enabled higher sensitivity in mutation detection leading to a more inclusive set of curated mutations. For these reasons, we present here only the results for the MB benchmark.

### Evaluation of sequencing library construction methods

Several different protocols were used for generating sequencing libraries at the five contributing sequencing centres, which varied in their reagent supplier, methods for selecting the fragment size of library inserts and use of amplifying PCR steps (Table 1 and Supplementary Table 1). Interestingly, these differences resulted in marked variation in the evenness of coverage genome-wide as well as in key regions of interest such as exons. PCR-free libraries were found to give the most even coverage, with very little effect of GC content on coverage levels (Fig. 1a). Evenness is directly correlated with coverage of the genome: when we downsampled each data set to an average 30 × tumour and control coverage (libraries sequenced to less than 28 ×, L.G and L.H, were excluded from further analysis), we see that in the best-performing libraries (L.A and L.B controls), 73–74% of the genome was covered at or above 25 ×, while the worst-performing library (L.F tumour) had only 46% of the genome covered at this level (Fig. 1b). In general, the coverage distribution was more even and the percentage of well-covered regions was higher in the control libraries compared with the tumour libraries, reflecting the different copy number states of the tumour. An unusual pattern of GC content distribution in control library E, however, meant that this was slightly worse than its tumour counterpart. The percentage of exonic regions covered at ≤10 × (that is, likely insufficient to accurately call mutations) also varied, with a range from less than 1% ‘missing’ in the best-performing libraries to more than 10% in the worst (Fig. 1c), demonstrating that sequencing library preparation performance can have a significant impact on the ability to identify coding variants in downstream analyses. Performance in other regions of interest, such as enhancers and untranslated repeats, was similarly variable (Fig. 1c and Supplementary Fig. 2).

### Evaluation of sequencing depth

Combining the sequencing data generated from each participating centre gave us the opportunity to investigate a tumour-normal pair with very deep coverage. After merging each of the individual pairs, the combined tumour coverage was 314 ×, and the control 272 ×. To remove already identified artefacts (Supplementary Note 1), we excluded the tumour library from centre E and the slightly contaminated control library from centre B. For comparison of mutation-calling metrics at a range of coverage levels, the combined tumour and normal sets were randomly serially downsampled to 250, 200, 150, 100, 50, 30 and 20 × coverage and then analysed using the standard DKFZ pipeline (MB.I, Supplementary Methods). The total number of mutations increases when going from 30 × to 50 × and further to 100 × coverage; however, no striking increase is seen above this level (at 100 ×, 95% of the maximum mutation number are detected, in contrast to only 77% at 30 × ; Fig. 2a and Supplementary Table 2). While the majority of mutations were called at the 30 × level, there were some notable differences in the number and type of mutations detected as the coverage increased. The sensitivity for detecting mutations with lower mutant allele frequencies (that is, subclonal alterations and/or events happening after polysomic changes but also major somatic mutations in samples with low tumour cell content) was much greater with higher coverage, as seen from density plots of mutations versus allele frequency (AF, Fig. 2b). This effect was even more striking when considering mutation calls per chromosome, which clearly shows the difference between low and high coverage when looking for late-occurring mutations after whole-chromosome copy number changes (Supplementary Fig. 3).

### Effect of tumour purity on mutation calling

Since MBs tend to show a very high tumour cell content (usually above 95%, and for this sample 98%, because of their nature as masses of small, round, tightly packed tumour cells), the high coverage data set also provided an opportunity to model the dynamics of mutation calling with increasing coverage and with increasing proportions of ‘contaminating’ normal tissue. We found that the mutation calls with increasing coverage were accurately modelled with a Michaelis–Menten equation, reaching ‘saturation’ (no or minimal additional mutations called as coverage increases) at around 100 × (Fig. 2c). The impact of normal cells on SSM detection could be thought of as a ‘mixed-type inhibition’ of mutation detection sensitivity, which we examined by mixing increasing proportions of normal sequence reads (17, 33 and 50%) into the tumour data set and re-calling mutations. Each curve displayed the same plateau after 100 × as the pure tumour sample; however, the addition of any normal content meant that the maximum mutation count from the pure tumour could not be reached, even at 250 × total coverage. At 100 ×, the detected proportions of mutation calls from the pure sample were 95%, 90% and 85%, respectively, for 17%, 33% and 50% ‘contamination’ (Fig. 2c). At lower coverage, the normal cell content had a proportionally larger impact. At 30 ×, only 92%, 83% or 68% of the calls from the 30 × pure sample were called when adding 17%, 33% or 50% normal reads, respectively (Supplementary Table 2). For SIMs called using the DKFZ pipeline, a different picture was observed. SIM calling at present likely suffers at least as much from low specificity as from low sensitivity, as indicated by the fact that increasing coverage actually reduces the number of called variants (that is, the FP rate decreases; Supplementary Fig. 4).

### Effect of tumour to normal sequencing depth ratio

We investigated the effect of tumour-normal coverage ratios on variant calling to assess whether increasing coverage of the tumour alone is sufficient to increase mutation detection sensitivity. The 250 × tumour genome was therefore compared with control genomes at 200, 150, 100, 50 and 30 × coverages. Down to the 150 × level, few differences are seen in the mutations called when compared with the 250 × /250 × standard (Supplementary Fig. 5a,b). At lower control coverage levels, a notable increase is observed in the overall number of mutations reported because of a sharp rise in those called with a low allele fraction. Since these mutations are not called in the 250 × versus 250 × set, it is almost certain that they are sequencing artefacts arising in a very small proportion of calls, which appear to be somatic when the control coverage is insufficient to show the same phenomenon. These new calls are dominated by a T>G base change arising in a particular sequence context (GpTpG, Supplementary Fig. 5c). Indeed, performing a motif analysis on the wider context of these changes revealed that the majority occur at a thymine base within a homopolymer run of guanines (Supplementary Fig. 5d). Keeping the ratio of tumour:normal coverage closer to one therefore appears to play a role in maintaining the accuracy of mutation calling with standard pipelines, since any systematic artefacts are then balanced out in both the tumour and control data sets. While it may be possible to apply additional filters to account for the new false positives (FPs) seen in unbalanced comparisons, this would potentially come at the cost of a reduced sensitivity for detecting true mutations with low allele frequencies (that is, tumour subpopulations), which are of particular interest when increasing sequencing coverage depth.

### Curation of a Gold Set of somatic mutations

We used the high-coverage (314 × :272 ×) data set from the sequencing benchmark to curate a Gold Set of verified somatic mutations (Methods). Gold Set mutations were classified (Table 2) according to the potential issues that may lead to an incorrect call: Tier 1 mutations have a mutant AF≥10%, Tier 2 mutations have AF≥5%, Tier 3 includes all verified mutations supported by unambiguous alignments, while Tier 4 includes additional mutations with more complicated or ambiguous local alignments and Tier 5 includes those with unusually high or low read depth.

The MB Gold Set had a total of 1,620 bona fide mutation calls across all tiers (Table 2), with 962, 1,101, 1,255 and 1,263 SSMs in Tiers 1, 2, 3 and 4, respectively, and 337 and 347 SIMs in Tiers 1 and 4, respectively. The mutational load of this tumour was 0.5 mutations per Mbp. Of these, there were eleven exonic SSMs (seven missense, three synonymous and one early stop) and one splice site mutation. We found that 32% of SSMs are in RepeatMasked sequence, 9% in tandem repeats (4% in homopolymer tracts) and 4.4% adjacent to tandem repeats. About a quarter of SSMs (27%) exhibits a mutant AF in the tumour of less than 10%, with 6% being very close to the alternate AF in the normal sample (Supplementary Table 3). For curated SIMs, 83% fall in tandem repeats (71% in homopolymers; Supplementary Table 4).

### Evaluation of somatic mutation calling pipelines

A submission and revision process was set up for the MB benchmark with guidelines for the mutation call format. Participating centres were provided with the best-quality sequence data set from the sequencing benchmark (L.A). Using these FASTQs, they produced SSM and SIM calls and submitted them for evaluation. We received 18 SSM and 16 SIM submissions.

Submissions were compared among themselves and to the Gold Set. Figure 3 shows the overlap of mutation call sets (private calls are shown in Supplementary Figs 6 and 7). We found that only 205 SSMs and one SIM were agreed upon by all submitters (Fig. 3 and Supplementary Figs 6 and 7). Agreement among SSM sets was much greater than agreement among SIMs in general. In Fig. 4 we show the precision versus recall of the submitted mutation calls. Each letter corresponds to a submission compared with the Gold Set of Tiers 1, 2 or 3, with comparison with Tier 1 (AF>10% mutant AF) having the highest value for recall in the plot, Tier 2 (AF>5%) the second highest and Tier 3 (AF>2%) the lowest. Precision is always calculated against Tier 4, which also includes mutations that are complex or have ambiguous positions.

We observed a cluster of well-performing SSM submissions with high values for both precision and recall. Those with the highest F1 scores (Table 3) were MB.Q and MB.J, pipelines that combine two different somatic mutation callers: qSNP23 with GATK24 and SGA25 with FreeBayes26 (Supplementary Methods). Submissions with a high number of calls did not necessarily achieve higher recall; about two-thirds of all mutations (or >80% of Tier 1 mutations) can be detected without making many false-positive (FP) calls, after which increases in recall are accompanied by precipitous declines in precision. This is because of the fact that at 40 × depth, a fraction of the curated mutations is impossible to detect. Likewise, the one submission whose precision was the highest (MB.L1) was not much more precise than MB.B or MB.Q, which both found over twice as many true mutations (Table 3). For SIMs (Fig. 4b), some submissions achieved precisions greater than 0.9; however, their sensitivities were still low. The highest F1 score (0.65 for MB.I) is noticeably lower than that obtained for SSMs (0.79). Overall, SIM detection appears to be more challenging, with performance lagging behind than that of SSM detection.

### Correlation of pipeline components with shared mutations

Could those submissions that cluster together in terms of precision and recall be calling the same mutations or have similarities in terms of their pipelines? Using a measure of pairwise overlap, the Jaccard index, we clustered the submissions and display the results as both a heatmap and hierarchical clustering in Supplementary Fig. 8. Correspondence analysis gave similar results. We also broke them down by true positives (TPs) or FPs (Supplementary Figs 9–12) and clustered the pipelines based on shared components or parameters (Supplementary Fig. 13, input data in Supplementary Data 1). We find that when submissions agree, they tend to agree on true mutations, and when they disagree, these are more likely to be FPs or true negatives. Some notable exceptions can be observed among the FP SIMs, where MB.L1 and MB.L2 cluster as do MB.F and MB.N. These concordant FPs may indicate incompleteness of the Gold Set and/or similarities in the pipelines. In this case, MB.L1 is a filtered subset of MB.L2, explaining the high degree of overlap. MB.F and MB.N both use Strelka27 to call SIMs, possibly explaining the overlap in FPs. Indeed, some overlap of MB.F and MB.N FP SIMs is seen with MB.L2, which also uses Strelka. For TP SIMs, pipelines that share components tend to have higher overlap, for example, among Strelka calls or among GATK SomaticIndelDetector calls. Sometimes, we observe higher Jaccard index values for pipelines using different software, for example Platypus28 and Atlas-indel29, which are two of the most sensitive mutation callers. There is much more concordance among SSM calls; therefore, trends are harder to see among the FPs. Logically, SSM submissions with the highest F1 scores have the highest Jaccard indices.

### Genomic or alignment features affecting accuracy

We asked what genomic features or sample, sequence or alignment features might correlate with better or worse performance. We used ‘rainfall’ plots that show density of features along a genome. By plotting the distance from the previous mutation (Fig. 5 and Supplementary Fig. 14), we can observe clustering or hotspots. SSM and SIM calls are coloured according to TP, FP and FN status. The Gold Set exhibits no mutational hotspots; therefore, any deviation is likely to be caused by a feature of the pipeline. Indeed, we detect quite different patterns: MB.Q (Fig. 5a) and MB.B (Fig. 5b) do not display any notable hotspots, while MB.C (Fig. 5c), for example, has many FP mutation calls in centromeric regions. MB.D and other call sets display varying degrees of this problem, which may arise if alignment of reads is performed to the GRCH37 reference without the d5 decoy sequence and/or no ultrahigh-signal blacklist is employed. MB.K overcalls (Fig. 5d) but a more subtle pattern is also apparent: higher ploidy chromosomes (for example, 17) display a greater density of calls and lower ploidy chromosomes (8, 15, X and Y) demonstrate a lower density of calls, presumably because of coverage.

Other genomic features such as tandem or interspersed repeats, as well as some key sample/sequence/alignment features, also create problems for mutation callers but are not detectable at the chromosomal scale. We annotated the Gold Set and all submitted SSM and SIM calls for each feature, indicating membership (Boolean flags) or a score (for continuous characteristics; Supplementary Tables 3 and 4). The frequencies or mean scores, respectively, were computed for three subsets of each submission (TPs, FPs or FNs) and for the Gold Sets. To highlight problematic features for each submission, the differences with respect to the Gold Set were computed and multiplied by the FP or FN rate, accordingly. The problematic features of FPs in the MB SSM data set are shown as a heat map (Fig. 6). While nearly all sets of FPs are enriched in low-frequency mutations, which are harder to discriminate from background noise (also reflected by the ‘same AF’ metric), some call sets (MB.K, H, C, D and M) do less well. MB.H seems to also have a problem with segmental duplications and multimappable regions, and MB.D with duplications only. MB.K and MB.M, to a lesser extent, are enriched in SSMs located in tandem repeats, simple repeats and homopolymers. MB.C has issues with FPs falling in blacklisted regions, specifically centromeric and simple repeats. The three submissions with fewer FPs immediately adjacent to tandem repeats than in the Gold Set (MB.H, MB.C and MB.D) do not use the Burrows-Wheeler Aligner (BWA) for the primary alignment step—instead, the mappers Novoalign or GEM are used, or the detection method is not based on mapping (SMUFIN30). The corresponding heatmaps for MB SSM FNs, and MB FP and FN SIMs are shown in Supplementary Figs 15–17. In general, both tandem repeats and segmental duplications and interspersed repeats cause sensitivity issues, with some pipelines more affected than others. The results for SIMs show that SIMs in tandem repeats (the majority being homopolymers) are undercalled, being under-represented in FPs and over-represented in FNs. Interestingly, nested repeats and duplications show the opposite trend, indicating that many FPs likely arise from low mapping quality.

Correspondence analysis confirms some of the above findings for MB SSM FPs (Supplementary Fig. 18). MB.C clusters with EncMap, dukeMap and centr, suggesting that MB.C FPs occur in some blacklisted regions. MB.H (and MB.G and MB.O to lesser extent) FPs are associated with segmental duplications. MB.K FPs are associated with tandem repeats (and microsatellites and simple repeats).

### Effect of mapper on mutation calling

The differences between sets of mutations submitted by the participating groups raised questions about the impact of individual pipeline components on the results. The extent of observed pipeline customization (Supplementary Methods and Supplementary Data 1) did not allow for exhaustive testing of all potentially important analysis steps; however, three pipeline components were selected for closer inspection because of their expected high impact: mapper, reference genome build and mutation caller. Four mappers (Novoalign2, BWA, BWA-mem and GEM), two SSM callers (MuTect31 and Strelka) and three versions of the human reference genome (b37, b37+decoy and ‘hg19r’—a reduced version of hg19, with unplaced contigs and haplotypes removed) were selected for testing, based on their usage by the benchmarking groups (Supplementary Methods for software versions and settings). To limit the effect of non-tested software on the produced mutation sets, a simple SSM-calling pipeline was established. First, we compared the effect of the mapper with each of the SSM callers. With a single SSM caller employed, a considerable fraction of unfiltered SSM calls for a given mapper (0.22–0.69, depending on the mapper–caller combination) is not reproducible by that caller with any other mapper (Supplementary Fig. 19). When compared with the Gold Set (Tier 3 SSMs), calls supported by a single mapper are almost exclusively FPs (precision <0.02). On the other hand, a large majority of calls supported by all four mappers are TPs (with precision ranging from 0.87 for MuTect to 0.99 for Strelka).

### Effect of primary mutation caller on mutation calling

Similar trends are observed when SSM callers are compared while holding the mapper constant (Supplementary Table 5). A sizable fraction (0.22–0.87, depending on the mapper) of unfiltered SSM calls for any given mapper–caller combination is not reproducible by the other caller on the same alignment file. Remarkably, in case of Novoalign2, the same alignment file leads to the most somatic calls and the lowest overall precision when used with MuTect, but the fewest somatic calls and highest overall precision when used with Strelka. When compared with the Gold Set, calls private to a single caller appear to be mostly FPs, with precision ranging from 0.01 to 0.05. Calls supported by both callers prove to be mostly correct (with precision between 0.89 and 0.93; Supplementary Table 6). The consensus sets seem to be robust—considerably improving the precision rates while only minimally lowering the sensitivity. The results of reference genome choice and a detailed examination of the alignment characteristics of the different aligners are presented in Supplementary Note 1.

### Improvement of pipelines using the benchmark data set

As a demonstration of the utility of the benchmark data set to improve pipelines, we set out to improve the MB.F pipeline further (already the MB.E pipeline, which uses SomaticSniper, was replaced with the MB.F pipeline, which uses Strelka, based on the analysis of the CLL benchmark results). Using the Gold Set to devise and tune a set of filters for various metrics (Supplementary Table 7), including mapping quality and distance from the end of the read alignment block (Supplementary Fig. 20), we are able to outperform (in terms of F1) all other MB SSM submissions (Fig. 7a and Supplementary Table 8). Despite choosing reasonably conservative thresholds, we were still worried about the possibility of overfitting; thus, we tested the adjusted pipeline on the CLL benchmark data set. We achieved similar results (Fig. 7b), demonstrating that the filter settings work well on at least one other cancer type. Removal of the repeat copy filter in Strelka also improve both MB and CLL SIM sensitivity without greatly affecting precision (Supplementary Fig. 21 and Supplementary Table 8).

Additional information regarding the CLL benchmark is provided in Supplementary Note 1, Supplementary Figs 22–36 and Supplementary Tables 9–13, which present the analogous information presented here for MB. Additional sequencing analyses are described in Supplementary Note 1, Supplementary Figs 37–40 and Supplementary Tables 14 and 15. Controls for genome reference builds and effect of mapper choice are presented in Supplementary Note 1, Supplementary Figs 41–47 and Supplementary Table 16. All pipeline details (as given by each submitter) are presented in Supplementary Methods, Supplementary Fig. 48 and Supplementary Tables 17 and 18.

## Discussion

This benchmarking exercise has highlighted the importance of carefully considering all stages of the laboratory and analysis pipelines required to generate consistent and high-quality whole-genome data for cancer analysis. In this study we have isolated and tested individual library construction/sequencing methods and complete analysis pipelines. Analysis pipelines themselves are also multicomponent; therefore, we have also evaluated mappers and two popular mutation callers in isolation.

By preparing libraries and generating sequence from the same MB tumour-normal pair at five different sequencing centres, we obtained results that suggest that PCR-free library preparation protocols should be the method of choice to ensure evenness of coverage, and that a sequencing depth of close to 100 × for both tumour and normal ought to be aimed for (particularly in situations where subclonal mutations or noncoding alterations are suspected to be playing a role). With platforms such as the Illumina HiSeq X now coming online in more centres, such an increase in coverage may be feasible without dramatically increasing costs.

This exercise also afforded a unique opportunity to compare validation schemes for benchmark creation. We found that high-depth (300 ×) WGS, in contrast to targeted resequencing, allowed us to more accurately assess FN rates in addition to FP rates, as well as better enable us to determine sweet spots of pipelines.

We have found that, contrary to common perception, identifying somatic mutations, be they SSMs or SIMs, from WGS data is still a major challenge. Calling mutations with different pipelines on differently prepared sequence read sets resulted in a low level of consensus. Using a standard pipeline had the potential of improving on this but still suffered from inadequate controls for library preparation and sequencing artefacts. Using a common high-quality sequence data set yielded higher concordance, but still resulted in substantial discrepancies in somatic mutation call rates and the calls themselves in the hands of different analysis groups. In some cases we were able to identify the Achilles’ heels of pipelines for failing to identify or for overcalling mutations as a function of genomic features or properties of sequencing depth, quality or alignment problems. We found that dominating features of both FP and FN SSMs were low coverage in the tumour, and aberrant coverage of tumour and normal. Underlying these artefacts are features such as segmental duplications and centromeric repeats. Use of an appropriate reference sequence (including decoy sequences) and/or good use of blacklists (problematic high-coverage regions or low mappability) can reduce FPs, while only an increase in overall sequencing depth or more sensitive algorithms (MuTect, for example) can address the FN rate. In contrast, we found that the vast majority of curated SIMs fall in simple/tandem repeats, and yet they are often filtered out because of a concern that they may be artefacts. We found little basis for this concern, at least in our data set that came from a no-PCR library. We found that adjustment of filters related to the number of copies of a repeat unit can increase sensitivity for this type of mutation.

Data analysis pipelines are constructed by integrating software from diverse sources. We found that the particular choice of a pipeline software is not as critical as how each piece of software is applied and how the results are filtered. In many instances, the approaches used in the different software and the assumptions made therein are not evident to the user and many parts are black boxes. We found that certain combinations show much higher compatibility than others, for example, with Novoalign alignments as input, MuTect produces an SSM call set with the lowest overall precision, while the SSM call set produced by Strelka on the same alignment file has the highest overall precision. Using combinations of tools for the same process step, assuming that a result shared by two different approaches has a higher likelihood to be correct, results in higher accuracy32. Indeed, we found that some of the higher accuracy pipelines utilize consensus of more than one mutation caller. Our controlled experiment intersecting Strelka and MuTect calls bore this out as well.

Recommended checklist for WGS cancer studies:

• PCR-free library preparation

• Tumour coverage >100 ×

• Control coverage close to tumour coverage (±10%)

• Reference genome hs37d5 (with decoy sequences) or GRCh38 (untested)

• Optimize aligner/variant caller combination

• Combine several mutation callers

• Allow mutations in or near repeats (regions of the genome likely to be more prone to mutation)

• Filter by mapping quality, strand bias, positional bias, presence of soft-clipping to minimize mapping artefacts

To account for many unknowns, variant calling pipeline developers often resort to calibrating their pipelines against known results, for example from a genotyping experiment performed on similar samples. This approach might only have limited validity genome-wide, as genotyping assays are typically biased towards less complex areas of the genome. We show that for cancer WGS experiments our benchmark set has the potential to be an unbiased and powerful calibration tool. The sequencing reads and curated Gold Set mutations described here are available to the research community through the ICGC DACO and the EGA to benchmark and calibrate their pipelines.

The issues that we have addressed in this study must be resolved (and we think they can, with the use of our benchmark data set) before WGS for cancer analysis can be wholly adopted for clinical use. However, we also suggest that further benchmarks be established to resolve even more difficult mutational features that tumour samples and genomes can present. These include, but are not limited to, low sample purity (contamination of tumour cells by normal cells and also the normal contaminated by tumour cells or viral components), subclonality, structural rearrangements and karyotype aberrations. Real cancers are complex and this complexity continues to challenge somatic mutation calling pipelines. In summary, this valuable resource can serve as a useful tool for the comparative assessment of sequencing pipelines, and gives important new insights into sequencing and analysis strategies as we move into the next big expansion phase of the high-throughput sequencing era.

## Methods

### Patient material

An Institutional Review Board ethical vote (Medical Faculty of the University of Heidelberg) as well as informed consent were obtained according to ICGC the guidelines (www.icgc.org). A limited amount of the original MB DNA can be made available on request.

### Library preparation and sequencing

The libraries were prepared at the different sequencing centres: the National Center for Genome Analysis (CNAG), Barcelona, Spain; the German Cancer Research Center (DKFZ), Heidelberg, Germany; the RIKEN Institute, Tokyo, Japan; the Ontario Institute for Cancer Research (OICR), Toronto, Canada, and the Wellcome Trust Sanger Institute, Hinxton, UK. Some libraries actually comprise a mixture of different libraries (as per the centre’s standard protocols); others comprise one library only. An overview of the composition of the different libraries and differences in the library preparation protocols is given in Table 1 and Supplementary Table 1. All samples were sequenced using Illumina technology and chemistry. The majority of reads are of 2 × 100 bp length and are derived from HiSeq2000 or HiSeq2500 sequencers, however, in one read set (L.A), a low number of 2 × 250-bp MiSeq reads are also included.

### Comparison of SSM calls

Each of the participating centres performed mutation calling using the respective in house pipelines (Supplementary Methods). The raw SSM calls were provided in the form of customized Variant Call Format (VCF) files. To provide a fair comparison, only single base point mutations were considered. A call was considered to be equal when both the position and the exact substitution reported were identical. The calls were then sorted according to the number of centres that made this particular call using a custom Perl script. The resulting file was plotted using a custom R-script (both available on request).

### Merging of the bam files to get the 300 × files

To create the high coverage 300 × bam files, the raw fastq files were aligned using bwa 0.6.2-r126-tpx aln -t 12 -q 20, followed by bwa-0.6.2-tpx sample -P -T -t 8 -a 1000 -r. The bam files for each centre/library were merged, and duplicates were marked using Picard tools MarkDuplicates Version 1.61. Finally, all merged per-centre bam files were merged using picard-1.95 MergeSamFiles and the header was adjusted using samtools-0.1.19 reheader. Since only reads from different libraries were merged at this step, duplicates were not marked. The coverage was calculated using an in-house tool, taking into account only non-N bases.

### Downsampling of the 300 × files

The 300 × bam files were serially downsampled to different coverage levels (250 ×, 200 ×, 150 ×, 100 ×, 50 ×, 30 ×, 20) using picard-1.95 DownsampleSam, and the coverage was determined after each step.

### Determination of library GC bias

To determine the GC bias of the libraries, we first created 10 kb windows over the whole genome using bedtools (v2.16.2) makewindows. Then, the GC content for each window was calculated using bedtools (v2.16.2) nuc. Windows containing more than 100 ‘N’ bases were excluded (awk-3.1.6 ‘BEGIN{FS=’\t’}{if ($10 <=100 &&$11 <=100) print $1"\t"$2"\t"$3"\t"$5}’). Finally, the coverage for each of the remaining windows was calculated using bedtools (v2.16.2) multicov. Since the total coverage of the different libraries was not the same, the coverage was normalized by dividing the coverage for each window by the mean coverage across all windows for each of the samples. To visualize the GC bias, we then plotted the normalized coverage against the GC content.

### Calculation of low coverage in special regions of interest

The regions of interest were defined as previously described33. To determine the percentage of bases covered with fewer than 10 reads, we first determined the coverage over the whole genome in per-base resolution using genomeCoverageBed (v2.16.2) -bga. The resulting coverage file was compressed using bgzip, and an index was produced with tabix-0.2.5 -p bed. We then extracted the coverage for our regions of interest using tabix-0.2.5. From the resulting extracted coverage files, we computed the number of bases covered by a certain number of reads using intersectBed and a custom perl script. This table was then used to determine the percentage of bases covered by ≤10 reads.

### Extraction of mutation signatures

Mutational catalogues were generated based on the somatic mutations detected in the tumours. The 3′ and 5′ sequence context of all selected mutations was extracted, and the resulting trinucleotides were converted to the pyrimidine context of the substituted base. Considering six basic substitution types with surrounding sequence context, this results in a mutation-type vector of length 96. The mutational catalogue was set up by counting the occurrence of each of these 96 mutation types per sample.

The proportions of the signatures published by Alexandrov et al.34,35 contributing to the mutational profile of each sample were estimated based on the probabilities of point mutations with their trinucleotide context in the signatures. The respective exposures were extracted sample-wise by quadratic programming. Exposures were plotted if they accounted for at least 5% of the SSMs in a sample.

### Somatic mutation calling benchmark data set

The sequencing reads provided to pipeline benchmark participants were produced at the CNAG using a no-PCR library preparation procedure that was adapted from the KAPA Library Preparation Kit protocol used together with Illumina TruSeq adaptors and omitting PCR amplification, each for the MB tumour and the corresponding normal DNA sample (L.A). For each sample two libraries were prepared with smaller (roughly 300 bases) and larger (roughly 450 bases) insert size. Sizing was performed using agarose gel separation and excision of corresponding size bands. The two tumour libraries were sequenced to 40.5 × and the two normal libraries to 29.6 × using a combination of Illumina HiSeq2000 (2 × 100 bp) and Illumina MiSeq (2 × 250 bp). MiSeq reads contributed about 2 × to each tumour and normal data. Reads in the FASTQ format were generated using the RTA software provided by Illumina.

### Evaluation of submissions

Automatic validation was performed on the submission server to minimize formatting problems. In addition, the submitted VCF files were sorted, filtered to restrict calls to chromosomes 1–22, X and Y and SIMs were left-aligned. Submissions of both CLL and MB SSMs and SIMs were evaluated against their respective Gold Sets, whose derivation is described above. For calculation of recall, the curated mutations were classified into three tiers according to alternate (mutation) AF. Only positions were considered, not the genotypes reported. For calculation of precision, all Tier 3 mutations plus ambiguously aligned mutations (class 4) were included so as to not penalize difficult to align but otherwise convincing differences between tumour and normal samples. For SIMs, no stratification into tiers 2 or 3 was performed; for recall, Tier 1 SIMs were used while, for precision, Tier 4 SIMs were used. To compare overall performance, we used a balanced measure of accuracy, the F1 score, defined as 2 × (Precision × Recall)/(Precision+Recall).

Overlap calculations for the purpose of clustering and heat map generation were performed using Tier 3 for SSMs and Tier 1 for SIMs. The Jaccard index is defined as the intersection divided by the union.

Correspondence analysis was performed using the ‘ca’ package in R on a table where each row corresponds to a genomic position at which at least one submission calls a somatic mutation in the Gold Set. The columns comprise the presence or absence in each call set, and Boolean values indicating whether certain genomic features such as repeats or presence in a blacklisted region apply, as well as sequence data such as AF or depth.

Rainfall plots represent the distance for each SSM call from its immediately prior SSM on the reference genome. For each submission, the SSM set used was made of SSM called classified as TP or FP and SSM from the Gold Set Tier 4 that were absent from the submission, classified as FN positions.

Feature analysis was conducted as follows. Tandem repeats were annotated with two programmes. Tandem repeats finder39 was run on 201-bp windows around each SSM and SIM calls. Any repeats greater than or equal to six repeat units overlapping or immediately adjacent to the mutation position was annotated accordingly. SGA25 was also used to annotate homopolymers specifically, giving a richer annotation of the repeat context and change induced by the mutation. Mappability was calculated using gem mappability40 with one mismatch at both 100mer and 150mer lengths. The average 100mer and 150mer mappabilities for each mutation were calculated for a window of −90 to −10 or −140 to −10 with respect to its position, respectively. Mult100 and mult150 are defined as 1—mappability 100 and 150. Same AF is defined as 1−(2 × (SCOREsnape-cmp-counts−0.5)) for SCORE≥0.5 else 0.

### Control of pipeline components

The protocol consisted of choosing a genome reference, mapping, alignment processing, variant calling and then analysis of alignments or variant calls.

Genome references tested:

• ‘b37d’ (‘human_g1k_v37_decoy’ from GATK bundle 2.8)

• ‘b37’ (‘human_g1k_v37’ from GATK bundle 2.8)

• ‘hg19r’ (‘ucsc.hg19’ from GATK bundle 2.3, with its unplaced contigs and haplotype chromosomes removed)

Mappers tested:

• Novoalign2 (v2.08.03) options: -i PE 360,60 -r All 10

• BWA (0.6.2-r126-tpx)

• BWA-mem (0.7.7-r441) options: -t 8 -M -B 3

• GEM (1.828)

Mutation callers tested:

• MuTect31 (v. 1.1.4; dbSNP v. 138; COSMIC v. 64)

• Strelka27 (1.0.13)

Alignment post-processing ensuring format compatibility with the downstream tools was performed as follows:

• Merge multiple SAM/BAM output files, coordinate-sorting of alignments with Picard tools (v. 1.84)

• Discard secondary alignments (alignments with SAM FLAG 0 × 100 set) with samtools (v. 0.1.18) (‘samtools view -h -b -F 256’)

• Mark duplicates with Picard tools’ MarkDuplicates.

• Realignment (RealignerTargetCreator, IndelRealigner) around indels with GATK (v. 2.3-9-ge5ebf34). The tumour and control were processed together.

• Apply Picard tools’ FixMateInformation

No base-quality recalibration or mutation filtration was applied. Two hundred fifty-one-nucleotide-long reads mapped by Novoalign were truncated to 150 nucleotides; this affected 4.4% of tumour reads (causing possible FNs because of missing mutation support in the tumour) and 4.5% of control reads (causing possible FPs because of missing mutation evidence in the control).

The programme qProfiler (http://sourceforge.net/p/adamajava/wiki/qProfiler/) was run on each mapper’s alignment files to investigate systematic mapping differences that potentially influenced subsequent mutation calling. Specifically, distributions of values in SAM fields ‘RNAME’ (alignment reference sequence ID), ‘MAPQ’ (alignment mapping quality score), ‘TLEN’ (observed insert size), ‘CIGAR’ (alignment-level indel details together with ‘soft clippings’—mapper-induced read trimming) and ‘MD’ (alignment mismatch details) were of interest. The same alignment files were used for mutation calling and qprofiler analysis. However, since mutation calling results were limited to chromosomes 1–22, X and Y, alignment files serving as qprofiler input were first filtered so as to contain only alignments to chromosomes 1–22, X and Y in order for the statistics being relevant (contig coverage statistics being the only exception, since mapping to the decoy contig appears to be important).

### Availability of data

Sequence data for this study have been deposited in the European Genome-phenome Archive (EGA) under the accession number EGAS00001001539.

How to cite this article: Alioto, T. S. et al. A comprehensive assessment of somatic mutation detection in cancer using whole-genome sequencing. Nat. Commun. 6:10001 doi: 10.1038/ncomms10001 (2015).

## References

1. 1

Hudson, T. J. et al. International network of cancer genome projects. Nature 464, 993–998 (2010).

2. 2

Mardis, E. R. & Wilson, R. K. Cancer genome sequencing: a review. Hum. Mol. Genet. 18, R163–R168 (2009).

3. 3

Ley, T. J. et al. DNMT3A mutations in acute myeloid leukemia. N. Engl. J. Med. 363, 2424–2433 (2010).

4. 4

Puente, X. S. et al. Whole-genome sequencing identifies recurrent mutations in chronic lymphocytic leukaemia. Nature 475, 101–105 (2011).

5. 5

Alkodsi, A., Louhimo, R. & Hautaniemi, S. Comparative analysis of methods for identifying somatic copy number alterations from deep sequencing data. Brief Bioinform. 16, 242–254 (2014).

6. 6

Dewey, F. E. et al. Clinical interpretation and implications of whole-genome sequencing. JAMA 311, 1035–1045 (2014).

7. 7

Kandoth, C. et al. Mutational landscape and significance across 12 major cancer types. Nature 502, 333–339 (2013).

8. 8

Jones, D. T. et al. Dissecting the genomic complexity underlying medulloblastoma. Nature 488, 100–105 (2012).

9. 9

Cancer Genome Atlas Research, N. Genomic and epigenomic landscapes of adult de novo acute myeloid leukemia. N. Engl. J. Med. 368, 2059–2074 (2013).

10. 10

Li, H. Toward better understanding of artifacts in variant calling from high-coverage samples. Bioinformatics 30, 2843–2851 (2014).

11. 11

McGinn, S. & Gut, I. G. DNA sequencing—spanning the generations. N. Biotechnol. 30, 366–372 (2013).

12. 12

Xu, H., DiCarlo, J., Satya, R. V., Peng, Q. & Wang, Y. Comparison of somatic mutation calling methods in amplicon and whole exome sequence data. BMC Genomics 15, 244 (2014).

13. 13

Highnam, G. et al. An analytical framework for optimizing variant discovery from personal genomes. Nat. Commun. 6, 6275 (2015).

14. 14

Zook, J. M. et al. Integrating human sequence data sets provides a resource of benchmark SNP and indel genotype calls. Nat. Biotechnol. 32, 246–251 (2014).

15. 15

Pabinger, S. et al. A survey of tools for variant analysis of next-generation genome sequencing data. Brief Bioinform. 15, 256–278 (2014).

16. 16

Fang, H. et al. Reducing INDEL calling errors in whole genome and exome sequencing data. Genome Med. 6, 89 (2014).

17. 17

O'Rawe, J. et al. Low concordance of multiple variant-calling pipelines: practical implications for exome and genome sequencing. Genome Med. 5, 28 (2013).

18. 18

Wang, Q. et al. Detecting somatic point mutations in cancer genome sequencing data: a comparison of mutation callers. Genome Med. 5, 91 (2013).

19. 19

Kim, S. Y. & Speed, T. P. Comparing somatic mutation-callers: beyond Venn diagrams. BMC Bioinformatics 14, 189 (2013).

20. 20

Louis, D. N. et al. The 2007 WHO classification of tumours of the central nervous system. Acta Neuropathol. 114, 97–109 (2007).

21. 21

Taylor, M. D. et al. Molecular subgroups of medulloblastoma: the current consensus. Acta Neuropathol. 123, 465–472 (2012).

22. 22

Ewing, A. D. et al. Combining tumor genome simulation with crowdsourcing to benchmark somatic single-nucleotide-variant detection. Nat. Methods 12, 623–630 (2015).

23. 23

Kassahn, K. S. et al. Somatic point mutation calling in low cellularity tumors. PLoS ONE 8, e74380 (2013).

24. 24

McKenna, A. et al. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 20, 1297–1303 (2010).

25. 25

Simpson, J. T. & Durbin, R. Efficient construction of an assembly string graph using the FM-index. Bioinformatics 26, i367–i373 (2010).

26. 26

Garrison, E. & Marth, G. Haplotype-based variant detection from short-read sequencing. Preprint at arXiv:1207.3907 (2012).

27. 27

Saunders, C. T. et al. Strelka: accurate somatic small-variant calling from sequenced tumor-normal sample pairs. Bioinformatics 28, 1811–1817 (2012).

28. 28

Rimmer, A. et al. Integrating mapping-, assembly- and haplotype-based approaches for calling variants in clinical sequencing applications. Nat. Genet. 46, 912–918 (2014).

29. 29

Challis, D. et al. An integrative variant analysis suite for whole exome next-generation sequencing data. BMC Bioinformatics 13, 8 (2012).

30. 30

Moncunill, V. et al. Comprehensive characterization of complex structural variations in cancer by directly comparing genome sequence reads. Nat. Biotechnol. 32, 1106–1112 (2014).

31. 31

Cibulskis, K. et al. Sensitive detection of somatic point mutations in impure and heterogeneous cancer samples. Nat. Biotechnol. 31, 213–219 (2013).

32. 32

Goode, D. L. et al. A simple consensus approach improves somatic mutation prediction accuracy. Genome Med. 5, 90 (2013).

33. 33

Rieber, N. et al. Coverage bias and sensitivity of variant calling for four whole-genome sequencing technologies. PLoS ONE 8, e66621 (2013).

34. 34

Alexandrov, L. B. et al. Signatures of mutational processes in human cancer. Nature 500, 415–421 (2013).

35. 35

Alexandrov, L. B., Nik-Zainal, S., Wedge, D. C., Campbell, P. J. & Stratton, M. R. Deciphering signatures of mutational processes operative in human cancer. Cell Rep. 3, 246–259 (2013).

36. 36

Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics 25, 1754–1760 (2009).

37. 37

Marco-Sola, S., Sammeth, M., Guigo, R. & Ribeca, P. The GEM mapper: fast, accurate and versatile alignment by filtration. Nat. Methods 9, 1185–1188 (2012).

38. 38

Raineri, E., Dabad, M. & Heath, S. A note on exact differences between beta distributions in genomic (Methylation) studies. PLoS ONE 9, e97349 (2014).

39. 39

Benson, G. Tandem repeats finder: a program to analyze DNA sequences. Nucleic Acids Res. 27, 573–580 (1999).

40. 40

Derrien, T. et al. Fast computation and applications of genome mappability. PLoS ONE 7, e30377 (2012).

## Acknowledgements

We thank the DKFZ Genomics and Proteomics Core Facility and the OICR Genome Technologies Platform for provision of sequencing services. Financial support was provided by the consortium projects READNA under grant agreement FP7 Health-F4-2008-201418, ESGI under grant agreement 262055, GEUVADIS under grant agreement 261123 of the European Commission Framework Programme 7, ICGC-CLL through the Spanish Ministry of Science and Innovation (MICINN), the Instituto de Salud Carlos III (ISCIII) and the Generalitat de Catalunya. Additional financial support was provided by the PedBrain Tumor Project contributing to the ICGC, funded by German Cancer Aid (109252) and by the German Federal Ministry of Education and Research (BMBF, grants #01KU1201A, MedSys #0315416C and NGFNplus #01GS0883), the BMBF-funded ICGC projects on early-onset prostate cancer and malignant lymphoma (#01KU1001A, #01KU1002B) and via the BMBF-funded de.NBI HD-HuB network (#031A537A, #031A537C); the Ontario Institute for Cancer Research to PCB and JDM through funding provided by the Government of Ontario, Ministry of Research and Innovation; Genome Canada; the Canada Foundation for Innovation and Prostate Cancer Canada with funding from the Movember Foundation (PCB). P.C.B. was also supported by a Terry Fox Research Institute New Investigator Award, a CIHR New Investigator Award and a Genome Canada Large-Scale Applied Project Contract. The Synergie Lyon Cancer platform has received support from the French National Institute of Cancer (INCa) and from the ABS4NGS ANR project (ANR-11-BINF-0001-06). The ICGC RIKEN study was supported partially by RIKEN President’s Fund 2011, and the supercomputing resource for the RIKEN study was provided by the Human Genome Center, University of Tokyo. M.D.E., L.B., A.G.L. and C.L.A. were supported by Cancer Research UK, the University of Cambridge and Hutchison-Whampoa Limited. S.D. is supported by the Torres Quevedo subprogramme (MICINN) under grant agreement PTQ-12-05391. E.H. is supported by the Research Council of Norway under grant agreements 221580 and 218241 and by the Norwegian Cancer Society under grant agreement 71220—PR-2006-0433. We specially thank Jennifer Jennings for administrating the activity of the ICGC Verification Working Group and Anna Borrell for administrative support.

## Author information

Authors

### Contributions

T.S.A., I.B., D.T.W.J. and I.G.G. planned the study and wrote the paper. T.S.A. coordinated mutation call submissions, led analysis team and performed primary analysis. I.B. and B.H. coordinated the sequencing benchmark and analysis. S.D. coordinated metadata collection and analysis. T.S.A., S.D., M.D.E., E.H., I.B., B.H., P.G., D.T.W.J., L.E.H., T.A.B., J.T.S., L.T., A.-S.S., A.-M.P., P.R., V.Q., R.V.-M., S.N., D.V., A.G.L., R.E.D., E.R., M.D., S.C.H., P.S.T., P.J.C., P.C.B., X.S.P., J.D.M. and I.G.G. contributed to analysis. S.D., M.D.E., E.H., I.B., B.H., B.B., R.D., R.K., S.G., A.K., D.T.W.J., L.E.H., T.A.B., J.T.S., L.T., A.-S.S., P.S.T., D.J., L.S., L.F., K.R., J.H., J.W.T., A.M., R.S., A.P.B., A.-M.P., P.R., V.Q., R.V.-M., S.N., D.V., L.B., A.G.L., C.L.A., N.J.H., T.N.Y., N.W., J.V.P., S.M.G., F.C.G., S.B., N.J., N.P., M.H., M.S., R.D., N.P., M.S., M.P., P.S., A.F., H.N., M.H., C.K., S.L., J.Z., L.L., S.M., S.S., D.T., L.X., D.A.W., C.L.-O., P.J.C., P.C.B. and X.S.P. contributed to pipeline development and formatting of submissions. E.C. and D.T.W.J. provided samples. S.C., S.S., N.D., C.P., H.N., A.F., X.S.P., R.E.D., J.D.M. and M.G. contributed to sequencing effort. P.B., J.D.M., S.M.P., R.E., P.L., D.G. and T.H. provided organization and additional input on the manuscript.

### Corresponding author

Correspondence to Ivo G. Gut.

## Ethics declarations

### Competing interests

The authors declare no competing financial interests.

## Supplementary information

### Supplementary Information

Supplementary Figures 1-48, Supplementary Tables 1-18, Supplementary Note 1, Supplementary Methods and Supplementary References (PDF 9288 kb)

### Supplementary Data 1

Gold Set mutations, accuracy of pipelines, pipeline parameter tables, and low coverage regions. (XLSX 233 kb)

## Rights and permissions

Reprints and Permissions

Alioto, T., Buchhalter, I., Derdak, S. et al. A comprehensive assessment of somatic mutation detection in cancer using whole-genome sequencing. Nat Commun 6, 10001 (2015). https://doi.org/10.1038/ncomms10001

• Accepted:

• Published:

• ### Uniform genomic data analysis in the NCI Genomic Data Commons

• Zhenyu Zhang
• , Kyle Hernandez
• , Jeremiah Savage
• , Shenglai Li
• , Dan Miller
• , Stuti Agrawal
• , Francisco Ortuno
• , Louis M. Staudt
• , Allison Heath
•  & Robert L. Grossman

Nature Communications (2021)

• ### A practical guide to cancer subclonal reconstruction from DNA sequencing

• Maxime Tarabichi
• , Amit G. Deshwar
• , Máire Ni Leathlobhair
• , Jeff Wintersinger
• , David C. Wedge
• , Peter Van Loo
• , Quaid D. Morris
•  & Paul C. Boutros

Nature Methods (2021)

• ### A unified haplotype-based method for accurate and comprehensive variant calling

• Daniel P. Cooke
• , David C. Wedge
•  & Gerton Lunter

Nature Biotechnology (2021)

• ### Short and long-read genome sequencing methodologies for somatic variant detection; genomic analysis of a patient with diffuse large B-cell lymphoma

• Hannah E. Roberts
• , Maria Lopopolo
• , Alistair T. Pagnamenta
• , Eshita Sharma
• , Duncan Parkes
• , Lorne Lonie
• , Colin Freeman
• , Samantha J. L. Knight
• , Gerton Lunter
• , Helene Dreau
• , Helen Lockstone
• , Jenny C. Taylor
• , Anna Schuh
• , Rory Bowden
•  & David Buck

Scientific Reports (2021)

• ### A Systematic Evaluation of High-Throughput Sequencing Approaches to Identify Low-Frequency Single Nucleotide Variants in Viral Populations

• David J. King
• , Graham Freimanis
• , Lidia Lasecka-Dykes
• , Amin Asfor
• , Paolo Ribeca
• , Ryan Waters
• , Donald P. King
•  & Emma Laing

Viruses (2020)


text:
# Lateral surface

Straight circular cylinder with unrolled lateral surface

The lateral surface of an object is the area of all the sides of the object, excluding the area of its base and top.

For a cube, the lateral surface area would be the area of four sides. Consider the edge of the cube as ${\displaystyle a}$. The area of one square face Aface = a ⋅ a = a2. Thus the lateral surface of a cube will be the area of four faces: a ⋅ a ⋅ 4 = 4a2. The lateral surface can also be calculated by multiplying the perimeter of the base by the height of the prism.[1]

For a right circular cylinder of radius r and height h, the lateral area is the area of the side surface of the cylinder: A = 2πrh.

For a pyramid, the lateral surface area is the sum of the areas of all of the triangular faces but excluding the area of the base.

For a cone, the Lateral Surface Area would be ${\displaystyle \pi rl}$ where ${\displaystyle r}$ is the radius of the circle at the bottom of the cone and ${\displaystyle l}$ is the lateral height (the length of a line segment from the apex of the cone along its side to its base) of the cone (given by the Pythagorean theorem ${\displaystyle l={\sqrt {r^{2}+h^{2}}}}$ where ${\displaystyle h}$ is the height of the cone)

## References

1. ^ Geometry. Prentice Hall. p. 700.


text:
### Welcome to our site!

Thanks a lot...
dufferdev

#### audioguru

##### Well-Known Member
I don't know where in the world you are, but here in Canada I can buy a cheap Chinese FM radio at The Dollar Store for $1.00. It comes with earphones, a battery and has a light that is a flashlight. It uses a Chinese copy of a Philips radio IC. You don't tune it, it scans the dial and you can stop the scan at a station. Its reception is for strong local stations only, it picks up a lot of interference and its sound is very distorted. You might be able to make one for$25.00 in parts if you can buy all the parts somewhere.

#### audioguru

##### Well-Known Member
I saw some crystal FM radio circuits on the internet. They are so simple that I laughed.

##### Well-Known Member
- is this 1 US dollar, or dollar from another country?
- you mentioned that it must be a circuit. So I assume then that we must include a speaker for the sound to come out, but is this correct? Or is it ok to simply offer a connection for IPOD headphones or similar?
- must we include the pcb in this budget?
- must we include the connectors in this budget?
- must we include the power supply in this budget?

My starting assumption is that we will not find a complete FM receiver IC to fit into our budget, so the design may have to be a simple discrete design using low cost transistors like 2N2222 or the like. This is why it is best to find the simplest designs.

#### dufferdev

##### New Member

The budget is USD....

- The output will be in the HEAD-PHONES no speaker required.... Our cost $0 - We have to include PCB, however currently we shall use a general breadboard.... A 6cm by 6cm breadboard costs$ 0.2 in India

- Which connectors ?

#### dufferdev

##### New Member

Can we have a circuit for the Scan and Reset function as it is available with those cheap FM radio ? If you have any such circuit please help. I shall try to source the components from Asia which shall be cheaper....

Thanks
dufferdev

Status
Not open for further replies.

Replies
2
Views
2K
Replies
34
Views
19K
Replies
4
Views
1K
Replies
17
Views
2K
Replies
10
Views
3K


text:
The operating results for Duff Corporation, an S corporation for last year were as follows: a. C...

The operating results for Duff Corporation, an S corporation for last year were as follows: a. Calculate the corporation's taxable income b. Determine and state the corporation's separately stated items and prepare the first page of the form 1120 c. Calculate the corporation's financial accounting income if the loss on the sale of the land is only $12,000 and depreciation is$52,000 under financial accounting rules


text:
# How do you solve the equation 2/3(x+8)^2-66=0?

Jul 18, 2017

$x = - 8 \pm 3 \sqrt{11}$

#### Explanation:

$\frac{2}{3} {\left(x + 8\right)}^{2} - 66 = 0$

$\frac{2}{3} {\left(x + 8\right)}^{2} = 66$ $\to$add $66$ to both sides

$\textcolor{b l u e}{\frac{3}{2}} \cdot \frac{2}{3} {\left(x + 8\right)}^{2} = 66 \cdot \textcolor{b l u e}{\frac{3}{2}}$ $\to$multiply both sides by $\frac{3}{2}$

${\left(x + 8\right)}^{2} = 99$

$x + 8 = \pm \sqrt{99}$ $\to$take the square root of both sides

$x + 8 = \pm 3 \sqrt{11}$ $\to$simplify the radical

$x = - 8 \pm 3 \sqrt{11}$ $\to$subtract $8$ from both sides

Jul 18, 2017

color(magenta)(x=-8+3sqrt11 or x=-8-3sqrt11

#### Explanation:

$\frac{2}{3} {\left(x + 8\right)}^{2} - 66 = 0$

$\therefore \frac{2}{3} \left({x}^{2} + 16 x + 64\right) - 66 = 0$

$\therefore \frac{2}{3} {x}^{2} + \frac{32}{3} x + \frac{128}{3} - 66 = 0$

$\therefore \frac{2}{3} {x}^{2} + \frac{32}{3} x + \frac{128}{3} - \frac{198}{3} = 0$

$\therefore \frac{2}{3} {x}^{2} + \frac{32}{3} x - \frac{70}{3} = 0$

multiply both sides by 3

$\therefore 2 {x}^{2} + 32 x - 70 = 0$

$\therefore 2 \left({x}^{2} + 16 x - 35\right) = 0$

$\therefore {x}^{2} + 16 x - 35 = 0$

$\therefore x = \frac{- b \pm \sqrt{{b}^{2} - 4 a c}}{2 a}$

$a = 1 , b = 16 , c = - 35$

$\therefore x = \frac{- 16 \pm \sqrt{{16}^{2} - 4 \left(1\right) \left(- 35\right)}}{\left(2\right) \left(1\right)}$

$\therefore x = \frac{- 16 \pm \sqrt{256 + 140}}{2}$

$\therefore x = \frac{- 16 \pm \sqrt{396}}{2}$

$x = \frac{- 16 \pm \sqrt{11 \cdot 3 \cdot 3 \cdot 2 \cdot 2}}{2}$

$\therefore x = \frac{- 16 \pm 6 \sqrt{11}}{2}$

$\therefore {\cancel{- 16}}^{\textcolor{m a \ge n t a}{- 8}} / {\cancel{2}}^{\textcolor{m a \ge n t a}{1}} \pm \frac{{\cancel{- 6}}^{\textcolor{m a \ge n t a}{- 3}} \sqrt{11}}{\cancel{2}} ^ \textcolor{m a \ge n t a}{1}$

$\therefore - 8 \pm 3 \sqrt{11}$

:.color(magenta)(x=-8+3sqrt11,x=-8-3sqrt11


text:
Tab Content
• Yesterday, 20:06
MarkFL replied to a thread Related Rates in Calculus
Yes, the length of his shadow is increasing at a rate of 4/3 ft/s, and the tip of his shadow is moving at a rate of 10/3 ft/s. :D
8 replies | 96 view(s)
• Yesterday, 16:29
MarkFL replied to a thread Related Rates in Calculus
I got 4/3 ft/s where you got 3 ft/s. You have the line: \d{y}{t}=\frac{3}{2}\d{x}{t} Now, this implies: \d{x}{t}=\frac{2}{3}\d{y}{t}
8 replies | 96 view(s)
• Yesterday, 16:20
Nicely done, Chris! (Yes) This is the solution I found elsewhere: We are given to evaluate: \lim_{n\to\infty}\left(\sqrt{\int_0^1...
4 replies | 114 view(s)
• Yesterday, 15:15
MarkFL replied to a thread Related Rates in Calculus
Suppose we wanted to work the problem in general terms. Please refer to the following diagram: $P$ is the height of the pole, $M$ is the...
8 replies | 96 view(s)
• February 26th, 2017, 23:47
Compute the following limit: \lim_{n\to\infty}\left(\sqrt{\int_0^1 \left(1+x^n\right)^n\,dx}\right)
4 replies | 114 view(s)
• February 26th, 2017, 02:44
I am reading Anderson and Feil - A First Course in Abstract Algebra. I am currently focused on Ch. 8: Integral Domains and Fields ... I need...
1 replies | 42 view(s)
• February 26th, 2017, 02:02
I am reading Anderson and Feil - A First Course in Abstract Algebra. I am currently focused on Ch. 8: Integral Domains and Fields ... I need...
1 replies | 41 view(s)
• February 26th, 2017, 01:48
We could likely model this with a function of the form (which comes from Newton's Law Of Cooling): f(t)=c_1e^{-kt}+350 Since we know $f(0)=70$,...
1 replies | 77 view(s)
• February 26th, 2017, 01:01
My solution: Consider the objective function: F(A,B,C)=\sin(A)+\sin(B)+\sin(C)-\cos(A)-\cos(B)-\cos(C) Subject to the constraint: ...
3 replies | 114 view(s)
• February 25th, 2017, 23:35
Let's look at the behavior of the solutions near the equilibrium solution $y=3$. Suppose we have: y(0)=2.9 Then we see: y'<0 And for:
4 replies | 78 view(s)
• February 24th, 2017, 11:34
Hello and welcome to MHB, Jesse! (Wave) I removed the link to your email address, as that shouldn't be public, and people can just as easily...
3 replies | 94 view(s)
• February 24th, 2017, 11:28
There's no need to link to the commercial site you used. Thanks! :D
3 replies | 99 view(s)
• February 24th, 2017, 00:26
Let's find the $p$th term. I think what I would do is set: \frac{n(n+1)}{2}=p Now solve for $n$... n(n+1)=2p
1 replies | 58 view(s)
• February 23rd, 2017, 22:15
Perhaps you meant how can we derive this formula. Let's begin with the definition of average acceleration: \overline{a}=\frac{\Delta v}{\Delta...
6 replies | 103 view(s)
• February 23rd, 2017, 20:43
I made a mistake in my above post...the IVP should be: \displaystyle \d{F}{t}=2-\frac{5F}{4(50-t)} where $F(0)=0$. So, this leads to: ...
7 replies | 125 view(s)
• February 23rd, 2017, 16:33
I would let $F(t)$ be the amount of fertilizer in the tank at time $t$. From the problem statement, we get the following IVP: ...
7 replies | 125 view(s)
• February 23rd, 2017, 14:31
I don't follow...it's already in the form we want to give us the acceleration. :D
6 replies | 103 view(s)
• February 23rd, 2017, 14:28
The kinematic formula you want is: \overline{a}=\frac{v_f^2-v_i^2}{2\Delta x}
6 replies | 103 view(s)
• February 23rd, 2017, 02:40
Our goal is to make MHB stand out as much as we can from our competitors by offering features you just don't see anywhere else. If we were to package...
3 replies | 128 view(s)
• February 23rd, 2017, 02:35
MarkFL replied to a thread Tapatalk... in Chat Room
Personally, I don't think Tapatalk is ready for vBulletin, much less sites using $\LaTeX$. From what I've gathered, many vBulletin admins refuse to...
4 replies | 94 view(s)
• February 23rd, 2017, 01:29
MarkFL replied to a thread Tapatalk... in Chat Room
I'd be happy if Tapatalk didn't lazily bypass most vBulletin PHP plugin hooks. (Giggle)
4 replies | 94 view(s)
• February 22nd, 2017, 20:01
Well! ... that is a thoroughly convincing proof!!! Thanks HallsofIvy ... really helpful ... Peter
5 replies | 119 view(s)
• February 22nd, 2017, 04:08
Your solution to the given IVP isn't correct...we are given: \d{y}{x}=9x^8e^{-x^9} where y(0)=8 \int_{y(0)}^{y(x)}\,du=\int_0^x...
1 replies | 66 view(s)
• February 21st, 2017, 23:42
Hello, MHB Community! (Wave) We are pleased to announce a new method for embedding Desmos calculators in your posts, which we hope you will find...
3 replies | 128 view(s)
More Activity

### 13 Visitor Messages

1. Happy Birthday, Bacterius!! I hope you have a nice celebration with your close ones. . Thank you for being an active part of MHB!
2. Hey Bacterius, really liked the number theory challenge of yours. I wonder if that's an Olympiad problem. (PS : I have produced somewhat of a solution, but I am not sure if it's correct)
3. Hey, Happy Birthday! Here in the U.S., the 21st birthday in really the last milestone of youth, because that's the legal drinking age.
4. I hope that today is the beginning of a great year for you...Happy Birthday!
5. For what it's worth, you come across on the forums as someone with wisdom and maturity far beyond your years.
6. Happy Birthday!

You are no longer a teenager...it's all downhill from here.
7. Happy Birthday, mate!!
8. Happy birthday!
9. Hey Bacterius!

I thought your post in the sequence challenge was very enlightening and well-written! Excellent job!
10. Hey,

Just wanted to say that your post explaining the Taylor series was very well written! Marvelous job!
Showing Visitor Messages 1 to 10 of 13
Page 1 of 2 12 Last
Page 1 of 2 12 Last
About Bacterius

#### Basic Information

Date of Birth
August 1, 1993 (23)
Location:
New Zealand
Country Flag:
New Zealand

#### Statistics

Total Posts
643
Posts Per Day
0.35
##### Thanks Data
Thanks Given
668
Thanks Received
1,585
Thanks Received Per Post
2.465
##### Visitor Messages
Total Messages
13
Most Recent Message
August 1st, 2015 00:06
##### General Information
Last Activity
June 19th, 2016 20:52
Last Visit
April 17th, 2016 at 02:19
Last Post
April 10th, 2016 at 04:23
Join Date
January 26th, 2012

### 18 Friends

1. #### alane1994Offline

MHB Craftsman

2. #### caffeinemachineOffline

MHB Journeyman

MHB Master

4. #### Chris L T521Offline

MHB Journeyman

5. #### cooliceey002Offline

MHB Apprentice

6. #### ellyseperryOffline

MHB Apprentice

7. #### FantiniOnline

MHB Craftsman

8. #### JamesonOffline

Администратор

9. #### MarkFLOnline

Pessimist Singularitarian

10. #### mathbalarkaOffline

MHB Journeyman

Showing Friends 1 to 10 of 18
Page 1 of 2 12 Last
No results to show...
Ranks Showcase - 1 Ranks
Icon Image Description

Name: Chat Box Champion (Jan-June 2013)
 Issue time: July 1st, 2013 22:09 Issue reason:


text:
# Calculate volume from 3d_plot

I read this interesting kernel on kaggle. I learned that you can make a 3D plot starting from CT sections. Is there a way to calculate the volume of the lungs or other structures you plot using 3d_plot? I couldn't find any reference so I need help.

Interesting question! I'm assuming you have a bunch of points $(X, Y, Z)$ plotting the lung you want to measure. So your question is actually the measure of a convex hull.
Checkout scipy.spatial.ConvexHull. The ConvexHull object has a volumn attribute which indicates the column of the convex hull you want. It supports N dimensions, which your question is about three dimensions.


text:
# What is the complexity of this algorithm for sparse matrices?

I am reading a paper on sparse matrices and there is an algorithm for sparse lower triangular systems. In the below pseudo-code $l$ is a sparse matrix and $x,b$ are sparse vectors.

x=b
for j=0 to n-1 do
if x[j] != 0
for i=j+1 to n-1 do
if l[i][j] != 0
~5 operations.

The paper claims that this algorithm is $\mathcal{O}(n+|b|+f)$, where $|b|$ is the number of non-zero entries in $b$ and $f$ is the number of floating point operations carried out.

I do not understand why this algorithm is not $\mathcal{O}(n|b|f)$.

Reference: Chapter 3 - Section 2 of https://epubs.siam.org/doi/book/10.1137/1.9780898718881

Some context on sparse matrices: Sparse matrices usually appear when solving problems using numerical schemes such as the finite element method. These matrices usually have a very small number of non-zero entries per row. For instance, a $n$-by-$n$ matrix might only have $\mathcal{O}(1)$ entries per row. The situation might not be as good as this in the algorithm above (there is no mention of a definition of sparse matrices in the book, as far as I can see.)

• Since the content is not freely available (one has to pay for it), is there any relevant context from the book you can provide us? It might be hard to answer this question without reading that chapter of the book. It's possible it is just a typo or an error; that happens.
– D.W.
Jul 6, 2018 at 14:53
• What is $l$? What is $b$?
– orlp
Jul 6, 2018 at 15:13
• I don't think it is a typo: The paragraph in which the statement arises refers the result two further times and a following paragraph once again. What other context would be needed? I am sorry if that comes across a bit obtuse, but I am not a computer scientist.
– user88950
Jul 6, 2018 at 16:11
• @orlp I have update the post.
– user88950
Jul 6, 2018 at 16:12

They say the running time is $O(n+|b|+f)$, where $f$ is the number of floating point operations taken during the entire computation. This takes some careful reading to understand what they mean.

They mean, if you take the total number of floating-point operations that are performed inside the innermost loop, and call that $f$, then the running time is not too much more than $f$. In other words, they are saying we do $O(n+|b|)$ operations for updating the loop indices, plus the time taken by the body of the innermost loop (summed over all the iterations).

The body of the innermost loop might indeed take $5 n|b|$ operations in total (say), in which case the total running time would be $O(n + |b| + 5n|b|) = O(n|b|)$, as you had expected.

Note that $f$ is not the number of floating point operations in a single iteration of the innermost loop; it is the total across all iterations of the loop. So $f$ is not 5.

Why do they count it this way? Presumably, the floating point operations are ones that we'll have to do no matter what, because they are essential to the operation you're doing; they instead want to focus on the extra book-keeping needed to update loop indices etc. and make sure that doesn't add too much extra cost.

Finally, I'll note that you changed the pseudocode. The pseudocode they show is not the same as what you listed in the question. In the book the pseudocode of the innermost loop is not "for $i=j+1$ to...", but rather "for each $i>j$ for which $l_{ij}$". That's not the same. Since $l$ is a sparse matrix/vector, there are more efficient ways to iterate over the non-zero entries of $l$ than iterating over all possible values of $j$ and testing whether that element is non-zero. This also affects the running time.

• Ahh yes I did change the pseudo-code. It was some odd way to get my head around what was going on. I have never seen this business of lumping all floating point operations together before. Thanks for your detailed response.
– user88950
Jul 6, 2018 at 20:24