암유전체에서차세대시퀀싱기반의 DNA 카피수변화발굴을위한 SOP (Standard Operating Protocols for Identification of NGS-Based DNA Copy Number Alterations in Cancer Genomes) 1
목차 1. 준비사항 (1) 배경 (3p) (2) 시퀀싱관련준비사항 (5p) 2. 시퀀싱데이터의전처리 (1) 유전체 alignment를위한주의사항 (8p) (2) Alignment를위한레퍼런스유전체 indexing (8p) (3) 레퍼런스유전체에대한 fastq파일의시퀀싱리드의 alignment (9p) (4) PicardTool기반전처리 (10p) (5) GATK기반전처리 (preprocessing) - #1 Local realignment (10p) (6) GATK기반전처리 (preprocessing) - #2 Score recalibration (11p) 3. VarScan을이용한 DNA카피수변화측정 (1) bam파일에서 Samtools를이용한 mpileup준비 (13p) (2) VarScan을이용한 tumor/normal리드수차이계산 (14p) (3) GC correction (15p) 4. Segmentation (CBS/circular binary segmentation) (16p) 5. Big-seq2를이용한 DNA카피수변화측정 (1) Modified Samtools를이용하여 uniquely mapping 된 read 탐색 (18p) (2) BicSeq-norm을이용한시퀀싱데이터의정규화 (normalization) (19p) (3) BicSeq-seg를이용한카피수 segmentation (20p) 6. Visualization (IGV browser) (22p) 7. GISTIC 분석 (24p) 8. 참고문헌 (28p) 2
1. 준비사항 (1) 배경 본 SOP에서다루는파이프라인은두종류의알고리즘 [VarScan2 (https://sourceforge.net/projects/varscan/files/) [1] 및 CBS (https://bioconductor.org/packages/release/bioc/html/dnacopy.html) [2] 를조합하여암유전체시퀀싱데이터에서 DNA카피수변화 (DNA copy number alterations) 를발굴하는방법에대한것으로, 전장유전체시퀀싱 (wholegenome sequencing혹은 WGS) 및유전체부분시퀀싱 ( 전장엑솜시퀀싱 - whole-exome sequencing/wes 및캡쳐리시퀀싱-captured resequencing) 데이터에폭넓게이용될수있음. 또다른방법으로전장유전체시퀀싱데이터로부터 DNA 카피수변화를얻 기위해구현된 Bic-seq2 [3] 를이용하여분석하기위한파이프라인도함 께다루고있음. 그림 1. DNA 카피수변화를측정하기위한전통적인마이크로어레이기법의도식화. Test/tumor DNA 를 Cy5 및 Reference/normal DNA 를 Cy3 로염색하여동시에 hybridization 할경우 (2-dye) tumor genome 의상대적인유전체증폭은높은 log2(tumor/normal) 을갖는 segment 로나타나며유전체결손은낮은 log2ratio 를갖는영역으로나타남. 3
유전체내에이미많이존재하는성선돌연변이 (germline alterations) 을효과적으로필터하고암유전체에특이적인체성돌연변이 (somatic mutation) 을발굴하기위해서종양유전체및해당환자의 matched normal유전체의시퀀싱쌍의조합이필수적임 ( 그림 1에서 test DNA 및 reference가각각 tumor DNA 및 matched normal DNA에해당함 ). 시퀀싱데이터를이용한 DNA카피수변화의발굴은크게두가지기법으로나눌수있음. 첫번째는소위 read-depth의변화를발굴하는것으로구획된유전체영역 (genomic bins) 에매핑되는종양유전체기원및정상유전체기원시퀀싱리드수의불균형 (biased presentation) 을측정하는것으로실제 VarScan2에서이용하는기법임 ( 그림 2). 이외, 유전체의구조적이상을시퀀싱리드에서직접적으로알아내는기법으로, paired read의 gap size 의불균형을통한유전체증폭및결손감지기법및 split리드의분석을통한염색체재조합의직접적인발굴기법임. 후자의경우, 본 SOP에서다루지않음. 그림 2. Read depth 의차이에의한 DNA 카피수의변화에따른 tumor 및 normal genome 의예. Read depth/coverage 의차이가염색체증폭 (copy gain) 및결손 (copy loss) 에따른변화양상을나타냄. 돌연변이 (point mutation/indel등 ) 와의차이점으로 DNA가아닌전사체 (RNA-seq등) 기반의시퀀싱데이터에서 DNA카피수를발굴하는것은일반적이지않으며본 SOP에서는 DNA(WGS및 WES) 시퀀싱데이터를대상으로함. 4
(2) 시퀀싱관련준비사항 본 SOP의경우종양유전체및해당환자의 matched normal유전체에서수행된전장유전체시퀀싱을대상으로하나, 필요에따라본 SOP는전장엑솜시퀀싱및리시퀀싱데이터를대상으로수행될수있음. 단, 종양유전체및정상유전체는동일한플랫픔으로시퀀싱이수행되어야하며, 또한시퀀싱규모 (depth혹은 coverage) 의경우종양-정상간차이가클경우통계적bias 를유발할수있음. 본 SOP에서는 30X의종양-정상전장유전체시퀀싱및 60-100X의종양-정상엑솜시퀀싱데이터를표준으로함 ( 그림 3). 그림 3. 시퀀싱 depth 가정해진 WGS/WES 기반의데이터에서특정영역에 align 된시퀀싱리드가주변혹은유전체전체영역에비해낮을경우 single copy (n=1) 혹은 homozygous (n=0) loss 를고려할수있으며이와반대로시퀀싱리드수가주변혹은유전체전체에비해높은경우유전체증폭을추론할수있다. 두종류의 n=1 혹은 n=0 결손만이가능하나, 유전체증폭의경우, single copy gain 에서부터특정영역에수백배수준으로증폭된 focal, high-level amplification 이가능할수있음. 5
본 SOP의경우현재표준적으로이용되고있는 Illumina사의 HiSeq기반의시퀀싱자료를표준대상으로함. 본 SOP를다른시퀀싱데이터 (IonTorrent/Proton외기타 NGS자료 ) 에적용시에도그에맞는변수조절이필요할수있음. 유전체작업에서사용되는유전체버전 (e.g., hg19) 을통일하는것은중요함. 현재표준적으로이용되는버전은 hg19 혹은그이후버전인 hg38임. 해당버전에맞는 reference genome을 alignment에서추후다양한 processing에이용하며, 특히 GATK bundle [4] 등의경우해당버전에맞는 variant reference등을제공하고있기때문에반드시통일된유전체버전을사용해야함. 본표준프로토콜은데이터의종류, 목적, 수행시기에따라많은변동사항 이있을수있으므로, 작성된프로토콜은모든연구를대변하는방법론이 될수는없음 본프로토콜은 2016년 10월현재를기준으로작성된것으로추후알고리즘버전업데이트등에의해프로토콜이수정되어야할수있으며, 보다좋은성능을보이는새로운알고리즘이개발되는등의경우에는프로토콜이변경될수있음. 그림 4 는 VarScan 을이용하여시퀀싱데이터로부터카피수변화를얻기 위한전체과정에대한모식도를보임. 6
그림 4. VarScan 과 CBS 를이용하여시퀀싱데이터로부터카피수변화를얻기위한과정모식도 7
2. 시퀀싱데이터의전처리 (1) 유전체 alignment를위한주의사항 유전체의 alignment를위해 BWA (Burrows-Wheeler aligner) [5] 등의표준 NGS aligner를사용함. 필요시, Bowtie, NovoAlign등의 BW-transformation 혹은 hash기반의 aligner로대체될수있음. 필요프로그램 : 본 SOP에서는 BWA에의한통상의 alignment를표준으로함. Alignment의경우 (1) reference sequence의 index file생성및 (2) 실제 fastq데이터의 alignment 후 SAM/BAM파일생성의두단계로나눌수있음. 본 SOP 에서는통상적으로이용되는 paired-end sequencing 의 alignment 를표준으로하나 single-end sequencing 의경우 alignment 과정의차이를 제외하고동일한 SOP 를이용할수있음. alignment 를포함한전처리과정은 WGS 과 WEX 에서 mutation 을발굴하 기위해수행하는초기과정과동일함. (2) Alignment 를위한레퍼런스유전체 indexing (BW transformation) Reference sequence 의 BW transformation 및 index 파일생성명령어 bwa index reference.fasta 해당파일은표준fasta포맷으로단일파일안에해당버전의전체유전체 (chr1, chr2,...) 가모두들어가도록함. hg19의유전체의경우 GATK bundle 의 ucsc.hg19.fasta를사용할수있음. GATK bundle의경우다양한유전체버전에맞는시퀀싱관련 reference file을제공하고있음. 해당 ftp의주소및억세스정보는다음과같음. 8
location: ftp.broadinstitute.org username: gsapubftp-anonymous password: <blank> GATK bundle의 hg19 유전체를사용할경우, "ucsc.hg19.fasta", " ucsc.hg19.fasta.fai ( 인덱스 )", " ucsc.hg19.dict" 의 3파일을다운로드하고해당 reference 파일을추후 GATK등의 processing에도이용하게함 ( 주의. alignment및 preprocessing의 reference가같아야함. 같지않은경우 preprocessing과정에서에러가발생할있음. 예, alignment에서 ucsc.hg19.fasta를이용하였으나염색체구성이다른 fasta파일을 preprocessing에서이용하는경우 ). (3) 레퍼런스유전체에대한 fastq파일의시퀀싱리드의 alignment Paired-end시퀀싱의경우 ( 본 SOP의표준포맷 ) 두개의쌍으로생성된 fastq를각각 alignment한후, 하나의 SAM파일로합침. SAM생성명령어는다음과같음. 예를들어, paired-end시퀀싱으로생성된 raw-fastq1.fq및 raw-fastq2.fq의경우다음의명령어를수행함. Single-end sequencing의경우단일 fastq를 mapping한후, bwa samse의명령어를수행함. bwa aln -f raw_1.sai $ref_dir/ucsc.hg19.fasta raw-fastq1.fq bwa aln -f raw_2.sai $ref_dir/ucsc.hg19.fasta raw-fastq2.fq bwa sampe -f output.sam $ref_dir/ucsc.hg19.fasta raw_1.sai raw_2.sai raw-fastq1.fq raw-fastq2.fq 상기명령어의조합은 bwa-mem을사용하여단일커맨드로수행할수있음. Initial alignment가수행된 SAM은이후, GATK에의한 preprocessing 을거침. 단, Read depth기반의 DNA카피수변화발굴의경우 SNV및 indel 의발굴처럼정밀한 read-position및 remapping을요구하지않는다는점에서 GATK의 preprocessing이필수적인과정은아님. 하지만암유전체분석의많은부분이 (SNV/indel calling 등 ) processing이끝난 BAM을요구한 9
다는점에서 preprocessing 이진행된 bam 파일에서 CNA 발굴및 mutation calling 을포괄적으로수행하는것이일반적인과정임. (4) PicardTool 기반전처리 (preprocessing) 준비물 : PICARD 최신버전 BWA 로 (local) align 된 SAM 파일의 (1) sorting, (2) bam 전환및 index 생성의 일련의과정을수행함. 이과정은 PICARD 툴을이용하며명령어는다음과 같음. java -jar $picard_dir/sortsam.jar SO=coordinate I=input.sam O=output.bam CREATE_INDEX=true PCR duplicate(technical 중복 ) reads 의삭제. 역시 PICARD 툴의 MarkDuplicate 옵션을사용하며명령어는다음과같음 java -jar $picard_dir/markduplicates.jar I=input.bam O=output.bam REMOVE_DUPLICATES=true (5) GATK 기반전처리 (preprocessing) - #1 Local realignment 크게 local realignment및 score recalibration과정으로이루어지며최신버전의 GATK및관련파일 ( 역시 GATK bundle로서얻을수있음 ) 이필요함. "Local realignment" 는 BWA기반의 genome-wide혹은 global realignment 에대응되는용어로서이미 BWA등의표준 aligner에의해 genome-align 된 read를주변에존재하는 indel정보에의거하여 local realignment를수행. GATK 의 local realignment 의경우 germline indel 을 input 파일로요구하며, 이미알려진 indel 의위치정보에기준하며 read 를 local 하게 realignment 함. 10
다음의명령어로수행할수있음 ( 첫번째명령어는 target_interval 을생성하 고, 이결과를이용하여두번째명령어로실제 realignment 를수행함 ). java -jar $GATK_dir/GenomeAnalysisTK.jar \ -R $ref_dir/ucsc.hg19.fasta \ -T RealignerTargetCreator \ -known $GATK_bundle/1000G_phase1.indels.hg19.vcf \ -known $GATK_bundle/Mills_and_1000G_gold_standard.indels.hg19.vcf \ -I input.bam \ -o output_intervals java -jar $GATK_dir/GenomeAnalysisTK.jar \ -R $ref_dir/ucsc.hg19.fasta \ -T IndelRealigner \ -targetintervals input_intervals \ -I input.bam \ -o output.bam 예제의경우, GATK 의표준 practice 에의거한 2 개의 indel file (-known option 과사용함. vcf 형태의 indel calls) 의예로다음의두파일역시 GATK bundle 로서얻을수있음. - 1000G_phase1.indels.hg19.vcf - Mills_and_1000G_gold_standard.indels.hg19.vcf (6) GATK 기반전처리 (preprocessing) - #2 Score recalibration 알려진SNP위치기준으로시퀀싱리드의매핑스코어를재조정함. 역시 GATK에기반한 2개의일련의명령어를통해수행함. 알려진 SNP정보로 GATK bundle에서제공하는 "dbsnp_137.hg19.excluding_sites_after_129.vcf" 파일을사용함. 중간단계로 recalibration을위한 covariate정보 (recal.grp) 를생성함. 11
java -jar $GATK_dir/GenomeAnalysisTK.jar \ -R $ref_dir/ucsc.hg19.fasta \ -T BaseRecalibrator \ -knownsites $GATK_bundle/dbsnp_137.hg19.excluding_sites_after_129.v cf \ -I input.bam \ -o recal.grp java -jar $GATK_dir/GenomeAnalysisTK.jar \ -R $ref_dir/ucsc.hg19.fasta \ -T PrintReads \ -BQSR recal.grp \ -I input.bam \ -o output.bam 12
3. VarScan 을이용한 DNA 카피수변화측정 (1) bam파일에서 Samtools를이용한 mpileup준비 (Samtools [6] 필요 ) VarScan은 bam형태의시퀀싱데이터에서 somatic/germline SNP/indel/CNA (DNA카피수변화) 를발굴하는통합프로그램패키지이며본 SOP에서는특히 somatic CNA발굴기능을다룸. 체성 (somatic) 돌연변이발굴을위해종양유전체및해당환자의정상유전체기반의시퀀싱데이터가쌍 (pair) 으로필요함. VarScan의 CNA발굴원리는 read depth의 local차이에기반한것으로 read depth를계산하기위해시퀀싱데이터내의각각의시퀀싱위치정보가필요함. Bam파일에서개별시퀀싱리드의위치정보추출을위해 Samtools의 mpileup을이용하여종양유전체및정상유전체에서각각시퀀싱리드의위치정보를추출하며해당명령어는다음과같음. samtools mpileup -q 1 -f $ref_dir/ucsc.hg19.fasta normal.bam >normal_mpileup samtools mpileup -q 1 -f $ref_dir/ucsc.hg19.fasta tumor.bam >tumor_mpileup normal.bam및 tumor.bam은비교분석하고자하는종양유전체및동일환자의정상유전체에서추출되고상기과정들에의해서 preprocessing이완료된 bam파일임. 특히, position/coordinate-sorted되어있어야함 (preprocessing과정참조 ). -q옵션을통해일정이하의 mapping quality를갖는시퀀싱리드를제거할수있음. 필요에따라 1(mapping quality 0만제거 - VarScan default) - 30 범위내에서설정할수있음. Score cutoff가높은경우 spurious reads의제거가가능하고노이즈를줄일수있으나, 충분한 coverage를얻지못할수있으므로 high-depth시퀀싱데이터에만적용하는것이권장됨. 13
(2) VarScan 을이용한 tumor/normal 리드수차이계산 두개의 mpileup 결과에대해서다음의명령을수행함. java -jar VarScan.jar copynumber normal_mpileup tumor_mpileup varscan.copynumber 상기명령에대한조정가능한옵션은다음과같음. 특히, genomic bin 사 이즈의조정이필요한경우 max-segment-size를 1000(1kb) 까지상향조정이가능함. genomic bin사이즈의경우, 클수록 local noise에 robust해지나적은크기의 focal amplification/deletion을무시할수있음. VarScan의해당과정에서조절할수있는옵션은다음과같음 ( 괄호안은 default수치임 ). --min-base-qual - Minimum base quality to count for coverage [20] --min-map-qual - Minimum read mapping quality to count for coverage [20] --min-coverage - Minimum coverage threshold for copynumber segments [20] --min-segment-size - Minimum number of consecutive bases to report a segment [10] --max-segment-size - Max size before a new segment is made [100] --p-value - P-value threshold for significant copynumber change-point [0.01] --data-ratio - The normal/tumor input data ratio for copynumber adjustment [1.0] 생성된 varscan.copynumer 파일의예는다음과같음. chrom chr_start chr_stop num_ positions normal_ depth tumor_ depth log2_ ratio gc_ content chr1 131367 131556 190 17.5 17.6 0.007 60.5 chr1 133395 133495 101 12 17.4 0.541 69.3 상기예는실제 log2ratio(tumor/normal) 이계산된 genomic bin 의일부를 나타내는것으로 genomic bin별로 (chrom/ 염색체 - chr_start/ 세그멘터시작위치 - chr_stop/ 세그먼트끝위치 ) normal depth( 매칭정상유전체에서해당 bin에관찰되는평균리드수 ), tumor depth( 종양유전체의해당bin에서관찰되는평균리드수 ), log2 ratio및해당 bin의시퀀스정보에기반한 GC content를나타냄. 14
(3) GC correction 생성된파일은 genomic bin별로 tumor및 normal genome기원의시퀀싱 depth/count가기록되어있으며실제 DNA카피수변화발굴을위해서다양한 genomic factor를조정할수있는데 VarScan에서는 GC correction을기본적으로권장하고있음. GC correction을위한명령어는다음과같음. java -jar VarScan.jar copycaller [varscan.copynumber] 상기명령어의결과로생성된파일의예는다음과같음. genomic bin 에대 한 GC fraction 정보및실제이러한 GC fraction 에의해보정된 genomic bin 의 log2(tumor/normal) read ratio 임. chrom chr_start chr_stop num_ positions normal_ depth tumor_ depth adjusted_ log_ratio gc_ content region_ call raw_ ratio chr1 657872 658079 208 21.8 18.8 0.049 54.3 neutral -0.212 chr1 761977 762364 388 84.6 92 0.369 53.9 amp 0.121 Genomic bin별로 GC content에따라 adjusted된 log ratio를제공하며, bin 별로기본적인 copy number call (amp - neutral - del) 을제시함. 실제 bin 별 copy number call은참고용으로사용되게되며, GC-corrected, adjusted log_ratio에대한 smoothing및 segmentation을별도의알고리즘으로수행하게됨. 15
4. Segmentation (CBS/circular binary segmentation) GC-correction ratio of genomic bin의경우다양한 segmentation알고리즘을통해 smoothing 및 CN (copy number) calling을할수있음. 현재까지알려진다양한 smooting/segmentation알고리즘중가장 performance가좋은것으로알려진 CBS(circular binary segmentation; https://bioconductor.org/packages/release/bioc/html/dnacopy.html) 알고리즘이일반적으로권장됨 ( 그림 5). 그림 5. 다양한 smoothing/segmention 알고리즘의 simulation 결과에서의 performance. Simulation data 를통해가장높은 performance 를보이는 smoothing/segmentation 알고리즘은 CBS 로알려져있고 TCGA 프로젝트등에서폭넓게사용됨. 단, Affymetrix SNP6.0 등의 noisy profile 의경우, GLAD(https://www.bioconductor.org/packages/release/bioc/html/GLAD. html) 를사용할경우좀더낳은결과를보일수있음. CBS/GLAD 모두 R package 가제공됨. CBS의경우 R package가제공되며 (DNAcopy R package; https://bioconductor.org/packages/release/bioc/html/dnacopy.html) Bioconductor에서얻을수있음. CBS의경우, 마이크로어레이기반및시퀀싱기반 (VarScan2 output) 을 smoothing/segmentation하며, noise를제거한 segmentation형태의결과를제공한다. CBS를이용해 VarScan GCcorrected file에직접적으로이용할수있는 R code는다음과같음. library(dnacopy) cn <- read.table("gc-corrected-varscan-output",header=f) CNA.object <-CNA( genomdat = cn[,6], chrom = cn[,1], 16
maploc = cn[,2], data.type = 'logratio') CNA.smoothed <- smooth.cna(cna.object) segs <- segment(cna.smoothed, verbose=0, min.width=2) segs2 = segs$output write.table(segs2[,2:6], file="out.file", row.names=f, col.names=f, quote=f, sep="\t") 상기명령어의결과는카피수프로파일의기본포맷인.seg 형태로 IGV(Integrated genome viewer; www.broadinstitute.org/igv/) 등의표준 browser 로시각화가가능함. *.seg 포맷형태의설명은그림 6 과같음. 그림 6. *seg 포맷형태. 단일카피수를보이는 segment (smoothing/segmentation 알고리즘적용후 ) 에속하는수십 - 수천개의 probe 를하나로묶어 single line 으로제공하며, 개별 probe 수준의 ( 여기서는 genomic bin) log2 profile 을기능적으로압축하는형태로, DNA 카피수프로파일의표준포맷으로사용됨. 해당예에서는 p1-p5, p6-10 가각각단일 seg line 으로묶임. seg.mean 의경우해당 probe 들의 log ratio 의평균값임. 데이터의 quality에따라서 noise를충분히제거할필요가있는경우상대적으로 smoothing이강한 segmentation알고리즘 ( 예, GLAD등 ) 을사용할수있음. GLAD역시 R package로제공되나 (GLAD), 생성되는 default output이 *.seg형태와다르기때문에변환용 script가필요함. 17
5. Bic-Seq2 를이용한 WGS 로부터의 DNA 카피수변화측정 그림 7. Bic-Seq2 를이용하여전장유전체시퀀싱에서카피수변화를얻기위한과정모식도 (1) Modified Samtools 를이용하여 uniquely mapping 된 read 탐색 (modified Samtools 필요 ) Bam 파일에서 uniquely mapping 된 read 들을추출하기위해 bic-seq2 에서 함께제공하는 modified samtools 이용 (samtools-0.1.7a_getunique-0.1.3). 종양유전체및정상유전체에서각각에해당하는명령어는다음과같음. samtools view -U BWA,dir_normal/,N,N normal.bam samtools view -U BWA,dir_tumor/,N,N tumor.bam 본예시는 BWA 를이용하여 mapping 한경우를설명하며, Bowtie 를이용 하여 mappiing 한경우에는 BWA 등 Bowtie 로적어주어야함. 18
(2) BicSeq-norm 을이용한시퀀싱데이터의정규화 (normalization) BicSeq2 에서는 GC-contents, nucleotide composition of short reads, mappability 정보를이용하여정규화를수행한후 copy number 를계산함. 정규화를위한명령어는다음과같음. BICseq2-norm.pl [options] <configfile><output> Options: -l=<int>: read length -s=<int>: fragment size -p=<float>: a subsample percentage (Default 0.0002). -b=<int>: bin the expected and observed as< int> bp bins (Default 100). --gc_bin: if specified, report the GC-content in the bins --NoMapBin: if specified, do NOT bin the reads according to the mappability --bin_only: only bin the reads without normalization --fig=<string>: plot the read count VS GC figure in the specified file (in PDF format) --title=<string>: title of the figure --tmp=<string>: the temp directory configfile 은다음과같이구성됨. chromname fafile MapFile readposfile binfilenorm chr1 chr1.fa hg19.50mer.crc.chr1.txt chr1.seq chr1.norm.bin chr2 chr2.fa hg19.50mer.crc.chr2.txt chr2.seq chr2.norm.bin 기본 (default) 옵션을이용하여 config 파일명을 'normalize_config' 라하고, output 명을 'normalize_output' 이라고할때, Big-seq 을이용한정규화의 실제수행예는다음과같음 perl NBICseq-norm.pl normalize_config normalize_output 19
(3) BicSeq-seg 를이용한카피수 segmentation BicSeq2 에서 segmentation 방법은기존 Bic-seq 과유사하게이웃한 bin 을 merging 하는것에의해수행됨. segmentation 을위한명령어는다음과 같음. BICseq2-seg.pl [options] <configfile> <output> Options: --lambda=<float>: the (positive) penalty used for BICseq2 --tmp=<string>: the temp directory --help: print this message --fig=<string>: plot the CNV profile in a PNG file --title=<string>: the title of the figure --nrm: do not remove likely germline CNVs (with a matched normal) or segments with bad mappability (without a matched normal) --bootstrap: perform bootstrap test to assign confidence (only for one sample case) --noscale: do not automatically adjust the lambda parameter according to the noise level in the data --strict: if specified, use a more stringent method to adjust the lambda parameter --control: the data has a control genome --detail: if specified, print the detailed segmentation result (for multisample only) paired 데이터에대한 configfile 은다음과같은형식으로구성됨. chromname binfilenorm.case binfilenorm.control chr1 CaseChr1.norm.bin ControlChr1.norm.bin chr2 CaseChr1.norm.bin ControlChr1.norm.bin normal과 tumor paired 데이터에대해서는반드시 --control 옵션을명시적으로적어주어야만 paired 데이터로인지하고수행됨. 기본 (default) 옵션을이용한 paired 데이터에대하여 config 파일의이름을 'seg_config', ouput 이름을 'seg_output' 이라고가정했을때실제수행방법의예는다음과같음 20
perl NBICseq-seg.pl --control seg_config seg_output 생성된 output 파일은표준 seg 파일의형식은아니므로, 필요시향후수 행될 visualization 이나 GISTIC 분석등을수행하기위하여표준 seg 파일 로변환과정이추가로필요할수있음. 21
6. Visualization (IGV browser) 생성된표준 *.seg는현재 genomic data의표준 browser인 IGV [7] 로변환없이 visualization이가능하며, 코호트데이터의경우 recurrent CNA를 GISTIC으로발굴할수있음. IGV browser의 snapshot은다음과같음 ( 그림 8). 그림 8. IGV browser snapshot. DNA 카피수변화프로파일의기본포맷인 *.seg 를직접적으로읽을수있음. Copy number 영역의 individual row 가특정환자 / 샘플을지칭하며, genome-wide 혹은 focal 염색체변화를 color change(red/gain - blue/loss) 혹은 bar graph 형태로시각화할수있음. 현재유전체데이터를읽는표준 browser로이용되고있는 IGV browser는 *.seg를직접적으로읽을수있음. IGV browser에서유전체버전 (hg19, hg38등 ) 을맞추고 drag-drop형식으로파일을읽음. IGV browser는해당홈페이지에서무료로얻을수있음 (http://software.broadinstitute.org/software/igv/). Log2형태의 raw 22
profile 의시각화가필요한경우, cn (copy number) 형태로전환후, IGV 로 시각화가가능함. IGV browser 는 genome-wide 및다양한수준의확대버전을제공함으로 써 genome-/chromosome-wide DNA 카피수변화및유전자수준의카피 수변화를단일혹은다수의샘플에서관찰할수있음. segmentation의수준을체크하기위해원래의 probe수준의데이터를그대로보존한 *.cn형태의파일을만들수있음. IGV browser에서는파일확장자에따라파일의형태를파악하기때문에 N (probe) x M (sample) 의 2D matrix형태를 txt로만들어 *cn형태로만들어 IGV browser에서읽을수있음. 보통암유전체분석을 DNA카피수분석외에 mutation(snv, indel등 ) 을동시에분석하게됨. BAM파일의시각화를통한 mutation의존재여부를 IGV로 manual검사를수행하며, 동시에 DNA카피수변화의검사가가능하며, TP53등의유전자에서흔히관찰되는 'double-hit/biallelic inactivation' 등의현상을관찰할수있음. 23
7. GISTIC 분석 다수의샘플에서 DNA카피수프로파일이확보될경우, GISTIC [8] 분석을통해유의하게반복적인소위 cancer driver의후보군을발견할수있음. GISTIC알고리즘은다수의샘플에서유전영역 /probe별로 log2ratio를합한 score를구하고, 이를 permutation을통해유의도를구하는알고리즘으로기본 schematics는그림 9와같음. 그림 9. GISTIC 알고리즘. 다수의샘플 (A, B, C) 에서나타나는 DNA 카피수변화를 summation 하고, 이를 permutation 을기반으로유의도를구하여, 유의하게반복적인 DNA 카피수변화를발굴하고, background random 변화와구분함. 현재 GISTIC 2.0은 MatLab에서수행할수있으며 (https://www.broadinstitute.org/cancer/cga/gistic) 또한 MatLab이없는경우에도 GenePattern을통해예제파일과함께 web기반으로손쉽게수행할수있음 (software.broadinstitute.org/cancer/software/genepattern/ 에서 GISTIC모듈선택 ). GISTIC을수행하기위해 *seg파일과 probe정보가필요함. 마이크로어레이데이터의경우해당플랫폼의 probe정보를그대로이용 (probe별염색체 / 위치정보 ) 할수있으나, 시퀀싱기반의 *seg파일의경우이러한정보가없으므로직접생성해야함. 보통 *.seg에존재하는모든 coordinate를추출하여이를염색체 / 위치정보별로 sorting하여 probe정보로이용할수있음. 24
GISTIC결과는주어진 probe별로염색체증가 / 감소를구분하여, 각 probe 별로염색체증가 / 감소에대한유의도 (FDR) 를구함. 보통 FDR < 0.25를표준값으로염색체증가 / 감소를구분하며특히 GISTIC 2.0의경우염색체수준 (arm-level) 및국소적 (focal) 변화에대한결과를제공함. GenePattern 하에서 GISTIC2.0 을돌리기위한인터페이스는다음과같음. *seg에대응하는유전체버전 (hg19등) 을선택한후, seg및 marker파일을업로드하여 GISTIC분석을수행할수있음. array list file은옵션으로 *.seg 파일에들어있는샘플중 subset을선택할수있고, cnv file의경우지정될경우, germline cnv를 filter할수있음. gene gistic을선택할경우, gene level의 DNA카피수 call을지정하여추후효과적인 gene-level의매칭분석을수행할수있음. seg 파일과 marker 파일을업로드한후 Run 버튼을클릭하면 GenePattern 하에서 GISTIC 을수행할수있음. GISTIC 실행결과다음과 같은파일이생성됨. 25
각파일에대한설명은다음과같음. o all_lesions 파일 : GISTIC 결과에대한요약파일로 amplification과 deletion 정보를포함한다. o amp_genes 파일 : GISTIC에의해 amplification된지역과이지역에관련된유전자정보를보인다. o del_genes 파일 : GISTIC에의해 deletion된지역과이지역에관련된유전자정보를보인다. o GISTIC score 파일 : amplification과 deletion 지역에대해 FDR 계산에의해얻어진 q value 값을 -log10(q) 형태로보인다. 또한 G- score, average amplitude, aberration frequency 등도함께보인다. o raw_copy_number 파일 : 입력으로받은 segmented copy number 데이터에대한 heatmap 그림을보인다. x 축은샘플이며, y 축은 genome을의미한다. o amp_qplot 파일 : amplification 지역에대해 GISTIC score와 q value 값을그림으로보인다. o del_qplot 파일 : deletion 지역에대해 GISTIC score와 q value 값을보인다. 26
GISTIC 결과에대한요약파일인 all_lesions 의결과는다음과같은형식으 로구성되어있음 Unique Name Amplification Peak 1 Amplification Peak 2 Amplification Peak 3 Descrip tor Wide Peak Limits 1q32.1 chr1:2015115 10-201914639(pr obes 7360:7368) 4q12 chr4:5459109 5-55217249(pro bes 29488:29500) 6p21.1 chr6:4261866 3-43350861(pro bes 46219:46225) Peak Limits chr1:2015121 99-201894643(pr obes 7361:7367) chr4:5460303 9-55203560(pro bes 29489:29499) chr6:4266481 7-43203014(pro bes 46220:46224) Region Limits chr1:2008732 70-202230155(pr obes 7337:7377) chr4:4883393 8-57740681(pro bes 29408:29589) chr6:4266481 7-43700259(pro bes 46220:46227) o Unique Name: 각영역에대항하는이름 o Descriptor: 각각에대한유전체영역설명 q values 4.06E -07 1.41E -10 0.219 51 Resid ual q values 4.06E -07 1.41E -10 0.219 51 Broa d or Foc al o Wide Peak Limits: 주로타겟유전자를포함하는피크영역 Amplitude Threshold 0: t<0.1; 1:0.1<t< 0.9; 2 t>0.9 0: t<0.1 1:0.1<t< 0.9; 2: t>0.9 0: t<0.1; 1:0.1<t< 0.9; 2: t>0.9 o Peak Limits: maximal amplification 혹은 deletion 이나타나는영역 o Region Limits: amplification 혹은 deletion 에대한유의한영역에 대한경계 o q values: 피크영역에대한 q value o Residual q-values: 오버랩등이일어난영역을제거한후에계산된 q-value o Broad or Focal: 넓은영역에서일어나는지여부 o Amplitude Threshold: 각샘플에대해 0, 1, 2 로표기되기위한 threshold 값 27
8. 참고문헌 [1] Koboldt DC, Zhang Q, Larson DE, Shen D, McLellan MD, Lin L, Miller CA, Mardis ER, Ding L, Wilson RK., VarScan 2: somatic mutation and copy number alteration discovery in cancer by exome sequencing, Genome Res. 22(3):568-576, 2012. [2] Olshen AB, Venkatraman ES, Lucito R, Wigler M. Circular binary segmentation for the analysis of array-based DNA copy number data. Biostatistics. 5(4):557-572, 2004 [3] Xi R, Lee S, Xia Y, Kim TM, Park PJ., Copy number analysis of whole-genome data using BIC-seq2 and its application to detection of cancer susceptibility variants. Nucleic Acids Res. 44(13):6274-6286, 2016. [4] McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, Garimella K, Altshuler D, Gabriel S, Daly M, DePristo MA. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 20(9):1297-1303, 2010. [5] Li H, Durbin R., Fast and accurate short read alignment with Burrows-Wheeler transform., Bioinformatics 25(14):1754-1760, 2009 [6] Li H., Handsaker B., Wysoker A., Fennell T., Ruan J., Homer N., Marth G., Abecasis G., Durbin R. and 1000 Genome Project Data Processing Subgroup, The Sequence alignment/map (SAM) format and SAMtools. Bioinformatics 25:2078-2079, 2009. [7] Robinson JT, Thorvaldsdóttir H, Winckler W, Guttman M, Lander ES, Getz G, Mesirov JP. Integrative genomics viewer. Nat Biotechnol. 29(1):24-26, 2011. [8] Mermel CH, Schumacher SE, Hill B, Meyerson ML, Beroukhim R, Getz G. 28
GISTIC2.0 facilitates sensitive and confident localization of the targets of focal somatic copy-number alteration in human cancers. Genome Biol., 12(4):R41, 2011. 29