AstraZeneca-NGS/Seq2C

Seq2C

Seq2C is Copy-Number Variations (CNV) variant caller. Seq2C normalizes the coverage from targeted sequencing to CNV log2ratio, and detects CNVs in regions that have abnormally high or low read depths compared to the rest.

Getting started

Requirements:

1. Program 'samtools' need to be in your PATH.
2. Perl version 5.8 and higher.
3. You need to add PATH variables to Seq2C scripts, description inthis section.

Getting source code

The Seq2C source code is located at https://github.com/AstraZeneca-NGS/Seq2C.

To load the project, execute the following command:

git clone https://github.com/AstraZeneca-NGS/Seq2C.git

Add PATH variables for scripts

Add short paths to the scripts seq2cov.pl, cov2lr.pl, lr2gene.pl, bam2reads.pl and others to your OS PATH. In Linux it can be done by adding this line to your .bashrc file:

export PATH=path_to_Seq2C_scripts_folder:$PATH

For using plotting script, also add the plotting folder:
export PATH=path_to_Seq2C_scripts_folder/plotting:$PATH

Run Seq2C

Usage:

seq2c.sh sample2bam.txt bed [control_sample1:[control_sample2]]

Parameters:

• sample2bam.txt - required, file must contains name of sample and paths to bam files with TAB delimiter (\t). Must contain one empty line in end of file. The example of the file:

NA12878.unmapped    /samples/NA12878.unmapped.bam
NA12889.mapped  /samples/NA12889.mapped.bam

• bed - required, path to BED file with TAB delimiter. Can be simple or amplicon. For simple mode it must contain 4 columns:

2   45904950    45904999    gene

For amplicon mode it must contain 8 columns. 7 column must contain amplicon start (must be more than start of region), 8 column must contain amplicon end (must be less than end of region).

2   45904950    45904999    gene    .   .   45904958    45904994

• control_samples - optional, for multiple controls, separate them using colon (:).

seq2c.sh script starts scripts from this repository in this order:

1. seq2cov.pl - Calculate candidate variance for a given region(s) in an indexed BAM file.
2. bam2reads.pl - Reads each read from sample and BAM file with samtools, creates statistics in read_stats.txt file.
3. cov2lr.pl - Normalizes the coverage from targeted sequencing to CNV log2 ratio.
4. lr2gene.pl - Converts a coverage file with CNV log2 ratio to copy number profile.

Output result (seq2c_results.txt file) contains columns:

1. Sample - sample name from original sample2bam file
2. Gene - gene name from BED file/region info
3. Chr - chromosome name from BED file/region info
4. Start - start of the segment/gene from BED file/region info
5. End - end of the segment/gene from BED file/region info
6. Length - length of the segment/gene
7. Log2ratio - CNV log2 ratio of normalized median depth by sample
8. Sig - significance (0 if AMP or DEL)
9. BP_Whole - type of CNV: Whole if AMP or DEL, else empty
10. Amp_Del - AMP (amplified) if log2 ratio more then -A option, DEL (deleted) if log2 ratio less then -D option
11. Ab_Seg - affected segments on gene
12. Total_Seg - total count of segments on gene
13. Ab_log2ratio - log2 ratio of normalized median depth by sample
14. Log2r_Diff - difference between CNV log2 ratio
15. Ab_Seg_Loc - segment location: ALL if AMP or DEL, else empty
16. Ab_Samples - abberation samples count for segments. Output only for BP calls that are “good”, i.e. wasn’t filtered by thresholds.
17. Ab_Samples_Pcnt - abberation samples fraction from total number of samples.

Command line options for scripts

You can run scripts separately if you want transfer additional command (for example, gender files or tumor purities). Parameters for scripts are listed in this chapter.

seq2cov.pl

Usage:
seq2cov.pl [-hz] [-n name_reg] [-b bam] [-c chr] [-S start] [-E end] [-s seg_starts] [-e seg_ends] [-x #_nu] [-g gene] [-o ori] [-d depth] region_info

The program will calculate candidate variance for a given region(s) in an indexed BAM file. The default input is IGV's one or more entries in refGene.txt, but can be regions (BED file).

Arguments are:

• region_info
  Required. IGV's one or more entries in refGene.txt, but can be regions (BED file).

Options are:

• -h
  Print this help
• -a int:float
  Indicate that it's PCR amplicon based calling. Each line in region_info represents a PCR amplicon (including primers). Two numbers in option are parameter to decide whether a particular read or pairs belongs to the amplicon.
  First is the number of base pairs.
  The second is the fraction of overlapped portion to the length of read or pairs.
  If the edges of reads (paired for Illumina) are within defined base pairs to the edges of amplicons and overlapped portion greater then the fraction, it's considered belonging to the amplicon.
  Suggested values are: 10:0.95. When given a 8 column amplicon format BED files, it'll be set to 10:0.95 automatically, but can still be overwritten by -a option.
• -n name_reg
  The regular expression to extract sample name from bam filename
• -N name
  Mutual exclusive to -n. Set the sample name to name
• -b bam
  The indexed BAM file
• -c chr
  The column for chr
• -S start
  The column for region start, e.g. gene start
• -E end
  The column for region end, e.g. gene end
• -s seg_starts
  The column for segment starts in the region, e.g. exon starts
• -e seg_ends
  The column for segment ends in the region, e.g. exon ends
• -g gene
  The column for gene name
• -x num
  The number of nucleotide to extend for each segment, default: 0
• -z
  Indicate whether it's zero based numbering, default is 1-based

bam2reads.pl

The program produces a file containing number of mapped or sequenced reads for samples.
Usage:
bam2reads.pl sample2bam.txt

Arguments are:

• sample2bam.txt file must contains name of sample and paths to bam files with TAB delimiter (\t).

cov2lr.pl

Usage: cov2lr.pl [-aH] [-c control] mapping_reads.txt coverage.txt

The program will convert a coverage file to copy number profile.

Arguments are:

• mapping_reads.txt
  Required. A file containing # of mapped or sequenced reads for samples. At least two columns. Can be created by bam2reads.pl script. First is the sample name, 2nd is the number of mapped or sequenced reads.
• coverage.txt
  The coverage output file from seq2cov.pl script. Can also take from standard in or more than one file.

Options are:

• -a
  Indicate this is amplicon or exon based calling. By default, it'll aggregate at gene level.
• -M
  Indicate to adjust the MAD when transforming the distribution. Default: no, or just simple linear function. If not sure, don't use this option. It might have better performance when cohort size is over 30.
• -c sample(s)
  Specify the control sample(s), if aplicable. Multiple controls are allowed, which are separated by ":"
• -F double
  The failed factor for individual amplicons. If (the 80th percentile of an amplicon depth)/(the global median depth) is less than the argument, the amplicon is considered failed and won't be used in calculation. Default: 0.2.
• -G Gender
  Take a file of gender information. Two columns, first is sample name, second is either M or F. If not provided, the program will try to guess.
• -Y chrYratio
  For gender testing, if chrY is designed. Default: 0.15. If chrY is carefully designed, such as Foundation's assay, default is good. For exome, the number should be higher, such as 0.3.
• -Z
  Indicate to output the frozen_file and all parameters into file Seq2C.frozen.txt

lr2gene.pl

The program will convert a coverage file to copy number profile. The default parameters are designed for detecting such aberrations for high tumor purity samples, such as cancer cell lines.
For clinical samples, many parameters need to be adjusted.

Usage:
lr2gene.pl [-aPH] [-cy] [-F float] [-s min_amplicon_#] [-A float] [-D float] cov2lr.txt

Arguments are:

• cov2lr.txt The coverage output file from cov2lr.pl script. Can also take from standard in or more than one file.

Options are:

• -c
  Indicate that control sample is used for normalization
• -y
  Debugging mode
• -s int
  The minimum consecutive amplicons to look for deletions and amplifications. Default: 1. Use with caution when it's less than 3. There might be more false positives. Though it has been successfully applied with option "-s 1" and identified one-exon deletion of PTEN and TP53 that were confirmed by RNA-seq.

For whole gene:

• -A float Minimum log2 ratio for a whole gene to be considered amplified. Default: 1.50
• -D float Minimum log2 ratio for a whole gene to be considered deleted. Default: -2.00

For < 3 exons:

• -E float Minimum mean log2 ratio difference for <3 exon deletion/amplification to be called. Default: 1.25
• -M float When considering partial deletions less than 3 exons/amplicons, the minimum MAD value, in addition to -E, before considering it to be amplified or deleted. Default: 10
• -d float When considering >=3 exons deletion/amplification within a gene, the minimum differences between the log2 of two segments. Default: 0.5
• -p float (0-1) The p-value for t-test when consecutive exons/amplicons are >= 3. Default: 0.00001

For break point in the middle of the gene:

• -e float When considering breakpoint in the middle of a gene, the minimum number of exons. Default: 5
• -t float When considering breakpoint in the middle of a gene, the minimum differences between the log2 of two segments. Default: 0.4
• -P float (0-1) The p-value for t-test when the breakpoint is in the middle with min exons/amplicons >= [-e]. Default: 0.000001

For cohort level aberrations:

• -R float (0-1) If a breakpoint has been detected more than "float" fraction of samples, it's considered false positive and removed. Default: 0.03, or 3%. Use in combination with -N
• -N int If a breakpoint has been detected more than "int" samples, it's considered false positives and removed. Default: 5. Use in combination with -R.

Additional scripts

Here is the information about scripts that can be used to generate gender files or convert the results of Seq2C.

seq2c2fm.pl

The program will parse seq2c output and make calls for each gene and output in the format compatible with OncoPrint. It has the option to provide the purity so that log2ratio thresholds will be adjusted accordingly.
By default, it calls genes that are homozygously deleted or amplified with >= 6 copies.

Usage:
seq2c2fm.pl [-g] [-e exons] [-n reg] [-N num] [-A amp] [-a amp] [-D del] [-d del] [-p purity_file] [-P purity] lr2gene_output

Arguments are:

• lr2gene_output
  The output file from seq2c.sh script (or lr2gene.pl). Can also take from standard in or more than one file.

Options:

• -k
  Print header
• -g
  Whether to output copy gains [4-5] copies. Default: no.
• -p file
  A file contains the tumor purity for all samples. Two columns, first is sample name, second is the purity in % or fraction [0-1].
• -P double
  The purity. Default: 1 or 100%, as is for cell lines. If set, all samples will assume to have the same purity.
• -n regex
  The regular expression to extract sample names. Default: none.
• -N num
  If an breakpoint has been called in >= num of samples, it's deemed false positive. Default: 5
• -e exons
  Minimum number of exons/amplicon. Default: 1

For whole gene:

• -D log2ratio
  The log2ratio to determine that a gene is homozygously deleted. Default: -2.0
• -A log2ratio
  The log2ratio to determine that a gene is amplified. Default: 1.45.

For < 3 exons:

• -d log2ratio
  The minimum log2ratio to determine that 1-2 exons are deleted. Should be lower than [-d] to reduce false positives. Default: -2.5
• -a log2ratio
  The minimum log2ratio to determine that 1-2 exons are amplified. Should be larger than [-a] to reduce false positives. Default: 1.75

For gains:

• -G Genes
  List of genes, seperated by ":", for which gain will be captured. Default: MYC

Output:
Some columns are filled with blanks, "NA" or "-" for compatibility with OncoPrint format.

1. Sample
2. Always ""
3. Variant_type: always "copy_number_alteration"
4. Gene
5. Always "NA"
6. Always "-"
7. Always "-"
8. Segment "chr:start"
9. Always "-"
10. Always "-"
11. Transf_LogRatio: transformed log_ratio as (2^LogRatio)*2
12. Segments: for BP - affected segments of total segments, for Whole - both numbers are total segments.
13. LogRatio
14. Alteration in format "amplification/gain/loss"
15. Always "-"
16. Always "-"
17. Always "-"
18. Always "-"
19. Always "-"
20. Always "-"
21. Always "-"
22. Alteration in format "Amplification/Gain/Deletion"

testGender.pl

The program will calculate the chrY coverage and make prediction of genders. The BAM file need to be targeted with chrY genes, exome or WGS. Otherwise, it might make wrong prediction.
Output can be used in cov2lr.pl script with -G option.

Usage:
testGender.pl [-hHy] [-B gender_BED] [-d dir] [-x cov] [-L len] [-b bam] [-n regex] [-N name] sample2bam.txt

Output: The output contains 7 columns:

1. The name of the sample
2. Predicted gender. Possible values are: 2.1 Male It's likely a male 2.2 Female It's likely a female 2.3 X,-Y It's likely only one X chr, but without Y chr. 2.4 Unknown The information is not enough to determine the gender
3. ChrY median depth
4. ChrX median depth
5. ChrA (autosome) median depth
6. p_value of t-test between chrX and chrA
7. The ratio between chrA/chrX (0.75-1.25 for Female, 1.65-2.35 for Male, roughly)

Limitations:

1. The bed file used for gender test should be targeted for targeted panels. No issue for WXS or WGS.
2. The chrY in bed file should be unique to chrY
3. For pure samples (e.g. cancer cell lines), males lost chrY or females lost one chrX can't not be differentiated for type "X,-Y". It should be less an issue for clinical sequencing, as there're almost always have normal cells in the sample.
4. For pure samples, males who lose chrY and but with amplified chrX can be mistaken as female.

Options:

• -H Print this help usage.
• -h Print the header.
• -y
  Print intermediate results (debuging purpose only).
• -B chrY_BED
  The BED file with CDS of chrY specific genes. Default to hg19_gender.bed in seq2c base directory.
• -d dir
  The installed seq2c base directory, if it's not in the path
• -x coverage
  The expected depth of coverage. Specify it unless the depth is > 10x. Default: 10.
• -L read_length
  The read length. Default: 100
• -b bam_file (optional)
  The bam file. For single sample prediction only, and no need to have sample2bam.txt file.
• -n regex
  Used with -b option. The regular expression to extract sample name
• -N string
  The sample name

EXAMPLES

1. Given a BAM file aligned to hg19, sample.bam, determine the gender: testGender.pl -b sample.bam -B hg19_gender.bed -N sample_name
2. If you have a file (samples.txt) with list of BAM files (tab delimited), then run the following command, assuming aligned to hg19:
   cat samples.txt | testGender.pl -B hg19_gender.bed

testrand.pl

The program will create matrix consists needed number of lines containg 5 random elements in each line. Sum of elements in each line will be the requested.

Usage:
testrand.pl sum_of_elements lines

Arguments are:

• sum_of_elements
  Sum of random elements in each line.
• lines
  Number of lines in constructed matrix.

Written by Zhongwu Lai, AstraZeneca, Boston, USA