精华内容
下载资源
问答
  • 基于磁性纳米粒子和化学发光的连接依赖性PCR的拷贝数变异分析
  • 09 拷贝数变异分析(GATK流程)

    千次阅读 2020-05-21 12:07:55
    09 拷贝数变异分析(GATK流程) 我们已经分析了 Somatic mutations,并进行了注释和可视化,接下来我们进行拷贝数变异的分析。 这里我们还是先从 GATK 的 somatic cnv 的最佳实践开始 拷 贝数变异(copy number ...

    09 拷贝数变异分析(GATK流程)

    我们已经分析了 Somatic mutations,并进行了注释和可视化,接下来我们进行拷贝数变异的分析。

    这里我们还是先从 GATK 的 somatic cnv 的最佳实践开始

    拷 贝数变异(copy number variations, CNVs)是属于基因组结构变异(structural variation, SV),是指 DNA 片 段长度在 1Kb-3Mb 的基因组结构变异。我们首先从 GATK 的 CNV 流程开始 CNV 的分析。

    首先是流程图,先从 bam 文件,结合坐标文件计算每个外显子的 reads counts 数,然后 call segment,最后是画图:
    在这里插入图片描述

    在下面这个链接中,给出了详细的教程:https://gatkforums.broadinstitute.org/gatk/discussion/11682#2

    或者最近更新的教程,一样的:https://gatk.broadinstitute.org/hc/en-us/articles/360035531092

    1. 外显子坐标的interval文件

    interval 文件其实也是一个坐标文件,类似 bed,只不过 bed 文件的坐标是从 0 开始记录,而 interval 文件的坐标是从1开始记录。使用基因组 interval 文件可以定义软件分析的分辨率。如果是全基因组测序,interval 文件就用全基因组坐标的等间隔区间就好。对于外显子组的数据,我们使用捕获试剂盒的目标区域,理论上应该返回原文查找对应试剂盒,去搜索其捕获外显子的区域,一般是外显子侧翼上下游 250bp 以内。我懒得去查,就用软件的默认参数 250bp。

    首先使用 BedToIntervalList 工具将 bed 转成 interval 格式(其实前面已经完成),然后用PreprocessIntervals工具获取 target 区间,即外显子侧翼上下游 250bp。

    GATK=~/wes_cancer/biosoft/gatk-4.1.4.1/gatk
    ref=~/wes_cancer/data/Homo_sapiens_assembly38.fasta
    bed=~/wes_cancer/data/hg38.exon.bed
    dict=~/wes_cancer/data/Homo_sapiens_assembly38.dict
    ref=~/wes_cancer/data/Homo_sapiens_assembly38.fasta
    
    ## bed to intervals_list
    $GATK BedToIntervalList -I ${bed} -O ~/wes_cancer/data/hg38.exon.interval_list -SD ${dict}
    
    ## Preprocess Intervals
    $GATK  PreprocessIntervals \
    -L ~/wes_cancer/data/hg38.exon.interval_list \
    --sequence-dictionary ${dict} \
    --reference ${ref}  \
    --padding 250 \
    --bin-length 0 \
    --interval-merging-rule OVERLAPPING_ONLY \
    --output ~/wes_cancer/data/targets.preprocessed.interval.list
    

    2. 获取样本的read counts

    这里分为两小步进行:

    • 首先是获取所有样本的 read counts,用到了工具是 CollectReadCounts, 其根据所提供的 interval 文件,对 bam 文件进行 reads 计数,可以简单理解为把 bam 文件转换成 interval 区间的 reads 数。最后会生成一个 HDF5 格式的文件,需要用第三方软件 HDFView 来查看,这里不做展示。这个文件记录了每个基因组interval 文件的 CONTIG,START,END 和原始 COUNT 值,并制成表格。
    • 然后构建正常样本的 CNV panel of normals,生成正常样本的 cnvponM.pon.hdf5 文件。对于外显子组捕获测序数据,捕获过程会引入一定的的噪音。因此后面需要降噪,该文件就是用于后面第 3 步 DenoiseReadCounts

    实际使用到的脚本如下:

    interval=~/wes_cancer/data/targets.preprocessed.interval.list
    GATK=~/wes_cancer/biosoft/gatk-4.1.4.1/gatk
    ref=~/wes_cancer/data/Homo_sapiens_assembly38.fasta
    
    cat config3 | while read id
    do
      i=./5.gatk/${id}_bqsr.bam
      echo ${i}
      ## step1 : CollectReadCounts
      time $GATK  --java-options "-Xmx20G -Djava.io.tmpdir=./"  CollectReadCounts \
      -I ${i} \
      -L ${interval} \
      -R ${ref} \
      --format HDF5  \
      --interval-merging-rule OVERLAPPING_ONLY \
      --output ./8.cnv/gatk/counts/${id}.clean_counts.hdf5
    
      ## step2 : Generate a CNV panel of normals:cnvponM.pon.hdf5
      $GATK  --java-options "-Xmx20G -Djava.io.tmpdir=./" CreateReadCountPanelOfNormals \
      --minimum-interval-median-percentile 5.0 \
      --output ./8.cnv/gatk/cnvponM.pon.hdf5 \
      --input ./8.cnv/gatk/counts/case1_germline.clean_counts.hdf5 \
      --input ./8.cnv/gatk/counts/case2_germline.clean_counts.hdf5 \
      --input ./8.cnv/gatk/counts/case3_germline.clean_counts.hdf5 \
      --input ./8.cnv/gatk/counts/case4_germline.clean_counts.hdf5 \
      --input ./8.cnv/gatk/counts/case5_germline.clean_counts.hdf5 \
      --input ./8.cnv/gatk/counts/case6_germline.clean_counts.hdf5 
    done
    

    这里的 input 重复了 6 次,为的是想让初学者容易理解。如果追求代码整洁,可以把这 6 句改成

    $(for i in {1..6} ;do echo "--input ./8.cnv/gatk/counts/case${i}_germline.clean_counts.hdf5" ;done)
    

    3. 降噪DenoiseReadCounts

    该步骤用到的工具是 DenoiseReadCounts,主要是做了一个标准化和降噪,会生成两个文件 ${id}.clean.standardizedCR.tsv${id}.clean.denoisedCR.tsv,该工具会根据 PoN 的 counts 中位数对输入文件 ${id}.clean_counts.hdf5 进行一个标准化,包括 log2 转换。然后使用 PoN 的主成分进行标准化后的 copy ratios 降噪。实际上这里只需要对 tumor 样本进行降噪的,为了比较,我就把 tumor 和 normal 样本都分析了一遍,后面可以做个对比。

    GATK=~/wes_cancer/biosoft/gatk-4.1.4.1/gatk
    cat config3 | while read id
    do
      i=./8.cnv/gatk/counts/${id}.clean_counts.hdf5
      $GATK  --java-options "-Xmx20g" DenoiseReadCounts \
      -I ${i} \
      --count-panel-of-normals ./8.cnv/gatk/cnvponM.pon.hdf5 \
      --standardized-copy-ratios ./8.cnv/gatk/standardizedCR/${id}.clean.standardizedCR.tsv \
      --denoised-copy-ratios ./8.cnv/gatk/denoisedCR/${id}.clean.denoisedCR.tsv
    done
    

    4. 可视化降噪后的copy ratios

    我们使用 PlotDenoisedCopyRatios 可视化标准化和去噪的 read counts。这些图可以直观地评估去噪的效果。

    GATK=~/wes_cancer/biosoft/gatk-4.1.4.1/gatk
    dict=~/wes_cancer/data/Homo_sapiens_assembly38.dict
    
    cat config3 | while read id
    do
      $GATK   --java-options "-Xmx20g" PlotDenoisedCopyRatios \
      --standardized-copy-ratios  ./8.cnv/gatk/standardizedCR/${id}.clean.standardizedCR.tsv \
      --denoised-copy-ratios ./8.cnv/gatk/denoisedCR/${id}.clean.denoisedCR.tsv \
      --sequence-dictionary ${dict} \
      --output ./8.cnv/gatk/cnv_plots \
      --output-prefix ${id}
    done
    

    对于每一个样本,如 case1_biorep_A_techrep,将生成 6 个文件:

    ./8.cnv/gatk/cnv_plots/case1_biorep_A_techrep.denoised.png
    ./8.cnv/gatk/cnv_plots/case1_biorep_A_techrep.denoisedLimit4.png  
    ./8.cnv/gatk/cnv_plots/case1_biorep_A_techrep.deltaMAD.txt        
    ./8.cnv/gatk/cnv_plots/case1_biorep_A_techrep.scaledDeltaMAD.txt
    ./8.cnv/gatk/cnv_plots/case1_biorep_A_techrep.standardizedMAD.txt
    ./8.cnv/gatk/cnv_plots/case1_biorep_A_techrep.denoisedMAD.txt
    

    其中4个文件是文本文件,里面就是一个数字,记录几个拷贝数变化比值 copy ratio 的中位数绝对偏差(median absolute deviation, MAD)

    ## 标准化后的 copy ratios 的 MAD
    $ cat ./8.cnv/gatk/cnv_plots/case1_biorep_A_techrep.standardizedMAD.txt
    0.229
    ## 降噪后的 copy ratios 的 MAD
    $ cat ./8.cnv/gatk/cnv_plots/case1_biorep_A_techrep.denoisedMAD.txt
    0.231
    ## 标准化后的 MAD 和降噪后的 MAD 的差
    $ cat ./8.cnv/gatk/cnv_plots/case1_biorep_A_techrep.deltaMAD.txt
    -0.002
    ## (降噪后的 MAD - 标准化后的 MAD ) / (标准化后的 MAD )
    $ cat ./8.cnv/gatk/cnv_plots/case1_biorep_A_techrep.scaledDeltaMAD.txt
    -0.01
    

    另 外两个文件是图片,表达同一个意思,标准化和降噪后的 copy ratios,底下还有中位数绝对偏差(median absolute deviation, MAD),只不过一张图片把 y 轴设置为 0 到 4。假如 tumor 样本的拷贝数没有发生变化,copy ratio 应该稳定在 1 附近。当然,要是发生了 CNV 事件,那应该就在 1 附近波动。

    在这里插入图片描述

    5. 计算常见的germline mutation位点

    这一步用到了 CollectAllelicCounts 工具,对输入的 bam 文件,根据指定的 interval 区间,进行 germline mutation 的检测(仅仅是 SNPs 位点,不包括 INDELs ),并计算该位点覆盖的 reads 数,即该位点的测序深度。值得注意的是,该工具一个默认参数是 MAPQ 值大于 20 的 reads 才会被纳入计数,最后生成一个 tsv 文件。

    GENOME=~/wes_cancer/data/Homo_sapiens_assembly38.fasta
    GATK=~/wes_cancer/biosoft/gatk-4.1.4.1/gatk
    interval=~/wes_cancer/data/targets.preprocessed.interval.list
    
    cat config3 | while read id
    do
      i=./5.gatk/${id}_bqsr.bam
      echo ${i}
      time $GATK  --java-options "-Xmx20G -Djava.io.tmpdir=./"  CollectAllelicCounts \
      -I ${i} \
      -L ${interval} \
      -R ${GENOME} \
      -O ./8.cnv/gatk/allelicCounts/${id}.allelicCounts.tsv
    done
    

    生成的 tsv 文件主要内容如下(这里我过滤掉了第六列为 N 的位点):

    $ less ./8.cnv/gatk/allelicCounts/case1_biorep_A_techrep.allelicCounts.tsv | grep -v ^@ | awk '{if($6 != "N") print $0}' |less
    CONTIG  POSITION        REF_COUNT       ALT_COUNT       REF_NUCLEOTIDE  ALT_NUCLEOTIDE
    chr1    925873  34      1       G       T
    chr1    925918  71      1       G       T
    chr1    925938  92      1       G       T
    chr1    925953  92      1       G       T
    chr1    925974  96      1       G       T
    chr1    925979  102     2       G       A
    

    6. ModelSegments

    在第 3 步我们拿到了标准化和降噪后的两个 tsv 文件,记录了某个区间的 LOG2_COPY_RATIO 值,内容大致如下:

    $ less ./8.cnv/gatk/denoisedCR/case1_biorep_A_techrep.clean.denoisedCR.tsv | grep -v ^@| less
    CONTIG  START   END     LOG2_COPY_RATIO
    chr1    925692  926262  -1.178123
    chr1    929905  930585  -0.569447
    chr1    930789  931338  -0.686712
    chr1    935522  936145  -0.404623
    chr1    938790  939201  -0.727205
    chr1    939202  939709  -0.882516
    

    而第 5 步拿到的记录等位基因测序深度的 tsv 文件已经在上面展示过了。接下来第 6 步将利用这两个结果进行 call segment,需要注意的是输入文件要求 tumor match normal。不过好像也可以不输入第 5 步 CollectAllelicCounts 的结果,等有时间再比较一下两者的区别吧。这一步用到的工具是ModelSegments,它根据去噪后的第三步的 reads counts 值对 copy ratios 进行分割,并根据第五步的CollectAllelicCounts等位基因计数对分割片段进行分类。代码如下:

    GATK=/~/wec_cancer/biosoft/gatk-4.1.4.1/gatk
    cat config3  | while read id
    do
      germline=${id:0:5}_germline
      ## ModelSegments
      $GATK   --java-options "-Xmx20g" ModelSegments \
      --denoised-copy-ratios ./8.cnv/gatk/denoisedCR/${id}.clean.denoisedCR.tsv \
      --allelic-counts ./8.cnv/gatk/allelicCounts/${id}.allelicCounts.tsv \
      --normal-allelic-counts ./8.cnv/gatk/allelicCounts/${germline}.allelicCounts.tsv \
      --output ./8.cnv/gatk/segments \
      --output-prefix ${id}
    
    done
    

    生成的文件有点多,每个样本生成 11 个文件。 param 文件包含用于 copy ratios(cr)和 allele fractions(af)的全局参数,而 seg 文件包含有关片段的数据。 具体说明可以查看这个链接:https://software.broadinstitute.org/gatk/documentation/tooldocs/4.1.4.0/org_broadinstitute_hellbender_tools_copynumber_ModelSegments.php

    ./8.cnv/gatk/segments/case1_biorep_A_techrep.af.igv.seg
    ./8.cnv/gatk/segments/case1_biorep_A_techrep.cr.igv.seg
    ./8.cnv/gatk/segments/case1_biorep_A_techrep.cr.seg
    ./8.cnv/gatk/segments/case1_biorep_A_techrep.modelBegin.af.param
    ./8.cnv/gatk/segments/case1_biorep_A_techrep.modelBegin.cr.param
    ./8.cnv/gatk/segments/case1_biorep_A_techrep.modelBegin.seg
    ./8.cnv/gatk/segments/case1_biorep_A_techrep.modelFinal.af.param
    ./8.cnv/gatk/segments/case1_biorep_A_techrep.modelFinal.cr.param
    ./8.cnv/gatk/segments/case1_biorep_A_techrep.modelFinal.seg
    ./8.cnv/gatk/segments/case1_biorep_A_techrep.hets.normal.tsv
    ./8.cnv/gatk/segments/case1_biorep_A_techrep.hets.tsv
    

    其不过上面拿到的文件中的 ${id}*.igv.seg 文件,可以直接载入到 IGV 中进行可视化,如:

    (不知道为什么,case4_biorep_B_techrepcase4_techrep_2的CNV事件是碎片化的,而 case5 病人的 X 染色体直接缺失)
    在这里插入图片描述

    7. CallCopyRatioSegments

    这一步用来判断 copy ratio segments 是扩增、缺失、还是正常的可能性。对上一步拿到的${id}.cr.seg进行推断,得到${id}.clean.called.seg文件会增加一列 CALL,用 +-0分别表示扩增、缺失和正常。基本上MEAN_LOG2_COPY_RATIO大于 0.14 就是扩增,小于 -0.15 就是缺失,其他的为正常。

    GATK=~/wes_cancer/biosoft/gatk-4.1.4.1/gatk
    cat config3 | while read id
    do
      $GATK   --java-options "-Xmx20g" CallCopyRatioSegments \
      -I ./8.cnv/gatk/segments/${id}.cr.seg \
      -O ./8.cnv/gatk/segments/${id}.clean.called.seg
    done
    

    8. 可视化CNV结果

    通过上面的分析,我们拿到了最后建模的 copy ratios 和 allele fractions segment,接下来用一个工具进行可视化:PlotModeledSegments

    dict=~/wes_cancer/data/Homo_sapiens_assembly38.dict
    GATK=~/wes_cancer/biosoft/gatk-4.1.4.1/gatk
    
    cat config3 | while read id
    do
      $GATK   --java-options "-Xmx20g" PlotModeledSegments \
      --denoised-copy-ratios ./8.cnv/gatk/denoisedCR/${id}.clean.denoisedCR.tsv \
      --allelic-counts ./8.cnv/gatk/segments/${id}.hets.tsv \
      --segments ./8.cnv/gatk/segments/${id}.modelFinal.seg \
      --sequence-dictionary ${dict} \
      --output ./8.cnv/gatk/cnv_plots \
      --output-prefix ${id}.clean
    done
    

    在这里插入图片描述

    整个gatk cnv流程

    上面整个流程的代码,其实可以合并为一个脚本 gatk_cnv.sh:

    GENOME=~/wes_cancer/data/Homo_sapiens_assembly38.fasta
    dict=~/wes_cancer/data/Homo_sapiens_assembly38.dict
    INDEX=~/wes_cancer/data/bwa_index/gatk_hg38
    GATK=~/wes_cancer/biosoft/gatk-4.1.4.1/gatk
    DBSNP=~/wes_cancer/data/dbsnp_146.hg38.vcf.gz
    kgSNP=~/wes_cancer/data/1000G_phase1.snps.high_confidence.hg38.vcf.gz
    kgINDEL=~/wes_cancer/data/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz
    interval=~/wes_cancer/data/targets.preprocessed.interval.list
    
    cd ~/wes_cancer/project/8.cnv/gatk
    ###################################
    #### 把bam文件转为外显子reads数 ######
    ###################################
    
    
    cat ~/wes_cancer/project/config3 | while read id
    do
      i=~/wes_cancer/project/5.gatk/${id}_bqsr.bam
      echo ${i}
      ## step1 : CollectReadCounts
      time $GATK  --java-options "-Xmx20G -Djava.io.tmpdir=./"  CollectReadCounts \
      -I ${i} \
      -L ${interval} \
      -R ${GENOME} \
      --format HDF5  \
      --interval-merging-rule OVERLAPPING_ONLY \
      --output ${id}.clean_counts.hdf5
      ## step2 : CollectAllelicCounts
      time $GATK  --java-options "-Xmx20G -Djava.io.tmpdir=./"  CollectAllelicCounts \
      -I ${i} \
      -L ${interval} \
      -R ${GENOME} \
      -O ${id}.allelicCounts.tsv
    done
    ### 注意这个CollectAllelicCounts步骤非常耗时,而且占空间
    
    mkdir allelicCounts
    mv *.allelicCounts.tsv ./allelicCounts
    mkdir counts
    mv *.clean_counts.hdf5  ./counts
    ##################################################
    # 接着合并所有的normal样本的数据创建 cnvponM.pon.hdf5 #
    ##################################################
    
    $GATK  --java-options "-Xmx20g" CreateReadCountPanelOfNormals \
    --minimum-interval-median-percentile 5.0 \
    --output cnvponM.pon.hdf5 \
    --input counts/case1_germline.clean_counts.hdf5 \
    --input counts/case2_germline.clean_counts.hdf5 \
    --input counts/case3_germline.clean_counts.hdf5 \
    --input counts/case4_germline.clean_counts.hdf5 \
    --input counts/case5_germline.clean_counts.hdf5 \
    --input counts/case6_germline.clean_counts.hdf5
    
    
    ############################################
    ############# 最后走真正的CNV流程 #############
    ############################################
    
    cat config | while read id
    do
      i=./counts/${id}.clean_counts.hdf5
      $GATK  --java-options "-Xmx20g" DenoiseReadCounts \
      -I $i \
      --count-panel-of-normals cnvponM.pon.hdf5 \
      --standardized-copy-ratios ${id}.clean.standardizedCR.tsv \
      --denoised-copy-ratios ${id}.clean.denoisedCR.tsv
    done
    
    mkdir denoisedCR standardizedCR segments cnv_plots
    mv *denoisedCR.tsv ./denoisedCR
    mv *standardizedCR.tsv ./standardizedCR
    
    cat config | while read id
    do
      i=./denoisedCR/${id}.clean.denoisedCR.tsv
      ## ModelSegments的时候有两个策略,是否利用CollectAllelicCounts的结果
      $GATK   --java-options "-Xmx20g" ModelSegments \
      --denoised-copy-ratios $i \
      --output segments \
      --output-prefix ${id}
      ## 如果要利用CollectAllelicCounts的结果就需要增加两个参数,这里就不讲解了。
    
      $GATK   --java-options "-Xmx20g" CallCopyRatioSegments \
      -I segments/${id}.cr.seg \
      -O segments/${id}.clean.called.seg
    
    
      ## 这里面有两个绘图函数,PlotDenoisedCopyRatios 和 PlotModeledSegments ,可以选择性运行。
    
      $GATK   --java-options "-Xmx20g" PlotDenoisedCopyRatios \
      --standardized-copy-ratios  ./standardizedCR/${id}.clean.standardizedCR.tsv \
      --denoised-copy-ratios $i \
      --sequence-dictionary ${dict} \
      --output cnv_plots \
      --output-prefix ${id}
    
      $GATK   --java-options "-Xmx20g" PlotModeledSegments \
      --denoised-copy-ratios $i \
      --segments segments/${id}.modelFinal.seg \
      --sequence-dictionary ${dict} \
      --output cnv_plots \
      --output-prefix ${id}.clean
    done
    

    对于每一个样本,就会拿到拷贝数变异的结果,如:

    在这里插入图片描述
    在这里插入图片描述在这里插入图片描述

    展开全文
  • CNV拷贝数变异分析是什么?贴一段TCGA官网的介绍 “The copy number variation (CNV) pipeline uses Affymetrix SNP 6.0 array data to identify genomic regions that are repeated and infer the copy number of ...

    CNV拷贝数变异分析是什么?贴一段TCGA官网的介绍

    “The copy number variation (CNV) pipeline uses Affymetrix SNP 6.0 array data to identify genomic regions that are repeated and infer the copy number of these repeats. This pipeline is built onto the existing TCGA level 2 data generated by Birdsuite and uses the DNAcopy R-package to perform a circular binary segmentation (CBS) analysis. CBS translates noisy intensity measurements into chromosomal regions of equal copy number. The final output files are segmented into genomic regions with the estimated copy number for each region. The GDC further transforms these copy number values into segment mean values, which are equal to log2(copy-number/ 2). Diploid regions will have a segment mean of zero, amplified regions will have positive values, and deletions will have negative values.”

    1. segment file数据下载和处理

    1.1 从TCGA下载数据

    下载文件类型:
    Copy Number Segment:A table that associates contiguous chromosomal segments with genomic coordinates, mean array intensity, and the number of probes that bind to each segment.

    Masked Copy Number Segment:A table with the same information as the Copy Number Segment except that segments with probes known to contain germline mutations are removed

    这里我用Masked Copy Number Segment做示范

    rm(list = ls())
    options(stringsAsFactors = F)
    options(scipen = 200)
    
    library(SummarizedExperiment)
    library(TCGAbiolinks)
    
    query <- GDCquery(project = "TCGA-BLCA",
                      data.category = "Copy Number Variation",
                      data.type = "Masked Copy Number Segment")
    GDCdownload(query,method = "api")
    BLCA_CNV_download <- GDCprepare(query = query, save = TRUE, save.filename = "BLCA_CNV_download.rda")
    

    1.2 数据处理

    #读取rda文件
    A=load("C:/Users/Meredith/Desktop/BLCA_CNV_download.rda")
    tumorCNV <- eval(parse(text = A))
    
    #改名
    tumorCNV=tumorCNV[,2:7]
    tumorCNV=tumorCNV[,c('Sample','Chromosome','Start','End','Num_Probes','Segment_Mean')]
    write.table(tumorCNV,file = 'BLCA_CNV.txt',sep = '\t',quote = F,row.names = F)
    

    BLCA_CNV.txt

    #提取01A结尾的样本(这里我用了python,小伙伴们可以用R来做)
    filename = 'BLCA_CNV.txt'
    finalResultName = 'segment_file.txt'
    
    read_file = open(filename)
    out_file = open(finalResultName,"r+")
    for line in read_file.readlines():
        data = line.split()
        x = data[0][13:16]
        if x == '01A':
            out_file.write(data[0])
            out_file.write('\t')
            out_file.write(data[1])
            out_file.write('\t')
            out_file.write(data[2])
            out_file.write('\t')
            out_file.write(data[3])
            out_file.write('\t')
            out_file.write(data[4])
            out_file.write('\t')
            out_file.write(data[5])
            out_file.write('\t')
            out_file.write('\n')
    

    segment_file

    2. marker file数据下载和处理

    2.1 从TCGA下载数据

    TCGA现在的参考基因组版本是hg38,需要从官网下载marker file。下载地址:GDC Reference File Website,选择最新版本的“SNP6 GRCh38 Remapped Probeset File for Copy Number Variation Analysis”文件,并注意“If you are using Masked Copy Number Segment for GISTIC analysis, please only keep probesets with freqcnv =FALSE

    官网下载marker file

    2.2 提取freqcnv=FALSE数据,并且整理成标准格式

    #这里我也是用的python(因为我电脑太菜用R要跑很久)
    filename = "snp6.na35.remap.hg38.subset.txt"
    finalResultName = "marker_file.txt"
    
    read_file = open(filename)
    out_file = open(finalResultName,"r+")
    for line in read_file.readlines():
        data = line.split()
        if data[5]=='FALSE':
            out_file.write(data[0])
            out_file.write('\t')
            out_file.write(data[1])
            out_file.write('\t')
            out_file.write(data[2])
            out_file.write('\t')
            out_file.write('\n')
    

    marker_file

    3. GenePattern GISTIC_2.0在线分析

    GenePattern
    GISTIC_2.0
    refgene file小伙伴们根据需要选择,这里我用的是TCGA下载的数据,所以选择hg38。将segment_file跟marker_file分别拖到seg file跟markers file区域。置信区间系统默认0.9,可以根据需要调整。点击RUN。
    要等半个小时左右。
    GISTIC_2.0输出结果

    我们来看下其中两个文件
    amp_qplot

    del_qplot

    4. maftools可视化GISTIC结果

    rm(list = ls())
    options(stringsAsFactors = F)
    
    BiocManager::install("PoisonAlien/maftools") #建议下载github最新版本的maftools包
    library(maftools)
    
    #读入GISTIC文件
    laml.gistic = readGistic(gisticAllLesionsFile = 'C:/Users/Meredith/Desktop/maftools/all_lesions.conf_99.txt',
                             gisticAmpGenesFile = 'C:/Users/Meredith/Desktop/maftools/amp_genes.conf_99.txt', 
                             gisticDelGenesFile = 'C:/Users/Meredith/Desktop/maftools/del_genes.conf_99.txt', 
                             gisticScoresFile = 'C:/Users/Meredith/Desktop/maftools/scores.gistic', 
                             isTCGA = T)
    

    readGistic里isTCGA=T作者给出的解释:please note that setting isTCGA=TRUE truncates sample names to first twelve characters - meaning primary samples and their metastatic counterparts (or from other sources) will all be collapsed into one single sample type. 所以在准备GISTIC输入文件时要将正常样品跟肿瘤样品分开。

    ##绘图
    #genome plot
    gisticChromPlot(gistic = laml.gistic, ref.build = 'hg38')
    

    A
    A. G-scores assigned by GISTIC for every cytoband plotted along the chromosome.

    #Bubble plot
    gisticBubblePlot(gistic = laml.gistic)
    

    B
    B. GISTIC results plotted as function of altered cytobands, mutated samples, and genes involved within the cytoband. Size of each bubble is according to -log10 transformed q values.

    #oncoplot
    gisticOncoPlot(gistic = laml.gistic, 
                   sortByAnnotation = TRUE, top = 10)
    

    C
    C. Oncoplot displays most frequently altered (amplifications or deletions) copy number events ordered according to the frequency. Each columns represents a sample and each row represent a CNV segment.

    #Visualizing CBS segments
    plotCBSsegments(cbsFile = 'C:/Users/Meredith/Desktop/maftools/segment.file.txt', #这里的segment file的sample name只保留前12位
                    tsb = 'ALL', #如果想指定某个样本,就输样本的名字,比如:tsb = 'TCGA-C4-A0F7'
                    ref.build = 'hg38')
    

    D
    D. Plots segmented copy number data.

    展开全文
  • 12 拷贝数变异分析(非GATK) CNVkit CNVkit 的用法比较简单,可以参考官网的教程:https://cnvkit.readthedocs.io/en/stable/index.html,写得很 详细,这里我不再说明。下面是我用到的脚本(在这里我比较关心的是...

    12 拷贝数变异分析(非GATK)

    CNVkit

    CNVkit 的用法比较简单,可以参考官网的教程:https://cnvkit.readthedocs.io/en/stable/index.html,写得很 详细,这里我不再说明。下面是我用到的脚本(在这里我比较关心的是最后拿到的 segment 文件,即 final.seg,该文件可以用直接载入到 IGV 进行可视化,当然也可以用 cnvkit.py export 生成这样的 segments 文件,两者基本一样,然后和上一节 GATK 的结果进行比较)

    ## cnvkit.sh
    GENOME=~/wes_cancer/data/Homo_sapiens_assembly38.fasta
    bed=~/wes_cancer/data/hg38.exon.bed
    
    cnvkit.py batch ./5.gatk/*[ABC]*bqsr.bam   ./5.gatk/*techrep_2_bqsr.bam  \
    --normal  ./5.gatk/*germline_bqsr.bam \
    --targets  ${bed} \
    --fasta $GENOME  \
    --drop-low-coverage --scatter --diagram --method amplicon \
    --output-reference my_reference.cnn --output-dir ./8.cnv/cnvkit
    
    cd ./8.cnv/cnvkit
    
    cnvkit.py export seg *bqsr.cns -o gistic.segments
    sed 's/_bqsr//' gistic.segments
    
    awk '{print FILENAME"\t"$0}' *bqsr.cns  | grep -v chromosome |sed 's/_bqsr.cns//g' |awk '{print $1"\t"$2"\t"$3"\t"$4"\t"$8"\t"$6}' >final.seg
    

    下面两张图是 CNVkit 结果和上一节 GATK CNV 流程的结果,对比两个图,发现这两个软件拿到的结果基本一样,同样的 case4 技术重复的两个样本拷贝数缺失事件碎片化。

    在这里插入图片描述
    在这里插入图片描述

    GISTIC2

    GISTIC2 的安装比较复杂,需要安装 MATLAB,具体可以参考曾老师的博客:用 GISTIC 多个 segment 文件来找SCNA 变异 http://www.bio-info-trainee.com/1648.html

    cd ~/wes_cancer/project/8.cnv/gistic
    ## 下载解压
    wget -c ftp://ftp.broadinstitute.org/pub/GISTIC2.0/GISTIC_2_0_23.tar.gz
    tar zvxf GISTIC_2_0_23.tar.gz
    
    ## 安装MATLAB
    cd MCR_Installer/
    unzip MCRInstaller.zip 
    ./install -mode silent -agreeToLicense yes -destinationFolder ~/wes_cancer/project/8.cnv/gistic/MATLAB_Compiler_Runtime/
    ## 安装后的文件结构
    
    ## 添加环境变量
    echo "export XAPPLRESDIR=~/wes_cancer/project/8.cnv/gistic/MATLAB_Compiler_Runtime/v83/X11/app-defaults:\$XAPPLRESDIR" >> ~/.bashrc
    source ~/.bashrc
    

    安装完成之后,新建一个目录,把上一步 cnvkit 拿到的 gistic.segments 文件拷贝到当前目录

    mkdir SRP070662_gistic
    cp  ~/wes_cancer/project/8.cnv/cnvkit/gistic.segments  ./SRP070662_gistic
    

    然后把将脚本 run_gistic_example 做了一定的修改,实际运行 GISTIC 的脚本 my.gistic.sh 如下:

    #!/bin/sh
    ## run example GISTIC analysis
    
    ## output directory
    echo --- creating output directory ---
    basedir=`pwd`/SRP070662_gistic
    # mkdir -p $basedir 
    
    echo --- running GISTIC ---
    ## input file definitions
    segfile=`pwd`/SRP070662_gistic/gistic.segments
    #markersfile=`pwd`/examplefiles/markersfile.txt
    refgenefile=`pwd`/refgenefiles/hg38.UCSC.add_miR.160920.refgene.mat
    #alf=`pwd`/examplefiles/arraylistfile.txt
    #cnvfile=`pwd`/examplefiles/cnvfile.txt
    ## call script that sets MCR environment and calls GISTIC executable 
    ./gistic2 -b $basedir -seg $segfile  -refgene $refgenefile  -genegistic 1 -smallmem 1 -broad 1 -brlen 0.5 -conf 0.90 -armpeel 1 -savegene 1 -gcm extreme
    

    生成的文件有

    all_data_by_genes.txt
    all_lesions.conf_90.txt
    all_thresholded.by_genes.txt
    amp_genes.conf_90.txt
    amp_qplot.pdf
    amp_qplot.png
    broad_data_by_genes.txt
    broad_significance_results.txt
    broad_values_by_arm.txt
    D.cap1.5.mat
    del_genes.conf_90.txt
    del_qplot.pdf
    del_qplot.png
    focal_dat.0.5.mat
    focal_data_by_genes.txt
    freqarms_vs_ngenes.pdf
    gistic_inputs.mat
    peak_regs.mat
    perm_ads.mat
    raw_copy_number.pdf
    raw_copy_number.png
    regions_track.conf_90.bed
    sample_cutoffs.txt
    sample_seg_counts.txt
    scores.0.5.mat
    scores.gistic
    wide_peak_regs.mat
    

    可 视化结果生成的 3 张图,3张图均未展示性染色体的拷贝数变异信息,图 1 垂直方向是基因组,水平方向是样本,如果倒过来,就和前面的 cnvkit 结果基本一致。图 2,3 表示的是拷贝数增加或者缺失的情况,水平方向下面是q值,绿线代表显着性阈值(q 值 = 0.25)。图上面是 G-score,某个区域的 G-score 分数越高,该区域的发生拷贝数变异事件的可能性就越大,偶然性越低。可以在输出文件 amp_genes.conf_90.txtdel_genes.conf_90.txt查看每一个区域所包含的基因。

    在这里插入图片描述
    在这里,拷贝数缺失事件明显比拷贝数增加事件要多,可能是 case4 技术重复异常带来的影响,我试着删除掉 case4 技术重复的两个样本,对比一下结果(如下)发现,拷贝数增加和缺失事件分布:
    在这里插入图片描述

    facet

    facet是一个 R 包,但不在 CRAN 中收录,安装的时候有点麻烦,在 Mac 和 Linux 上安装没问题,在 windows 上安装没成功。 facet 的输入文件是处理后的 vcf 文件,在这里我用的是 GATK 的HaplotypeCaller工具生成的结果并走了gvcf流程拿到的 merge.vcf,最后载入 R 中进行分析, R 语言分析部分的代码如下:

    if(F){
      install.packages('devtools')
      devtools::install_github("mskcc/facets", build_vignettes = TRUE)
      devtools::install_github("mskcc/pctGCdata")
    }
    options( stringsAsFactors = F )
    library("facets")
    # check if .Random.seed exists
    seedexists <- exists(".Random.seed")
    # save seed
    if(seedexists) 
      oldSeed <- .Random.seed
    # Alway use the same random seed
    set.seed(0xfade)
    # read the data
    if(F){
      library(vcfR)
      vcf_file='./5.gatk/gvcf/merge.vcf'
      ### 直接读取群体gvcf文件即可
      vcf <- read.vcfR( vcf_file, verbose = FALSE )
      save(vcf,file = 'input_vcf.Rdata')
    }
    load(file = 'input_vcf.Rdata')
    vcf@fix[1:4,1:4]
    vcf@gt[1:14,1:4]
    colnames(vcf@gt)
    library(stringr)
    tmp_f <- function(x){
      #x=vcf@gt[,2]
      gt=str_split(x,':',simplify = T)[,1]
      ad=str_split(str_split(x,':',simplify = T)[,2],',',simplify = T)[,2]
      dp=str_split(x,':',simplify = T)[,3]
      return(data.frame(gt=gt,dp=as.numeric(dp),ad=as.numeric(ad)))
    }
    colnames(vcf@gt)
    tms=colnames(vcf@gt)[!grepl('_germline',colnames(vcf@gt))][-1]
    tmp <- lapply(tms, function(tm){
      #tm=tms[1]
      print(tm)
      nm=paste0(strsplit(tm,'_')[[1]][1],'_germline')
      print(nm)
      if(nm %in% colnames(vcf@gt)){
        snpm=cbind(vcf@fix[,1:2],
                   tmp_f(vcf@gt[,nm]),
                   tmp_f(vcf@gt[,tm]))
        ## only keep tumor or normal have mutations
        kp=snpm[,3] %in% c('1/1','0/1') | snpm[,6] %in% c('1/1','0/1')  
        table(kp)
        snpm=snpm[kp,]
        ## remove those show ./. positions
        kp=snpm[,3] == './.' | snpm[,6]== './.'  
        print(table(!kp))
        snpm=snpm[!kp,c(1,2,4,5,7,8)]
        rcmat=snpm
        rcmat$POS=as.numeric(rcmat$POS)
        dim(rcmat)
        rcmat=na.omit(rcmat)
        colnames(rcmat)=c("Chromosome", "Position","NOR.DP","NOR.RD","TUM.DP","TUM.RD")
        rcmat[,1]=gsub('chr','',rcmat$Chrom)
        ## fit segmentation tree
        xx = preProcSample(rcmat)
        ## estimate allele specific copy numbers
        oo=procSample(xx,cval=150)
        oo$dipLogR
        ## EM fit version 1
        fit=emcncf(oo)
        tmp=fit$cncf
        write.table(tmp,file = paste0('facets_cnv_',tm,'.seg.txt'),
                    row.names = F,col.names = F,quote = F)
        head(fit$cncf)
        fit$purity
        fit$ploidy
        png(paste0('facets_cnv_',tm,'.png'),res=150,width = 880,height = 880)
        plotSample(x=oo,emfit=fit, sname=tm)
        dev.off()
        if(F){
          fit2=emcncf2(oo)
          head(fit2$cncf)
          fit2$purity
          fit2$ploidy
          png(paste0('facets_cnv2_',tm,'.png'),res=150,width = 880,height = 880)
          plotSample(x=oo,emfit=fit2, sname=tm)
          dev.off()
    
        }
        return(c( fit$dipLogR,fit$ploidy,fit$purity))
      }
    
    })
    tmp <- do.call(rbind,tmp)
    rownames(tmp)=tms
    colnames(tmp)=c('dipLogR', 'ploidy', 'purity')
    write.csv(tmp,'ploidy_and_purity.csv')
    

    我们就拿case1_biorep_A_techrep样 本的结果进行展示,主要关心图的第三个 panel ,绘制每个 segment 的主要(黑色 =2 )和次要(红色 = 1)拷贝数。底部栏显示相关的细胞分数(cf)。深蓝色表示高 cf。浅蓝色表示低 cf。米色表示正常。这个结果就和上面略有不同了,比如这个样本的 19 号染色体在第三个 panel(第一、二个panel 好像还是正常的)中显示是拷贝数增加,但在上面几个软件的结果中 19 号染色体的拷贝数均没有太大的变化。毕竟 facet 的输入数据是 germline mutations,和前面即个软件都不太一样。
    在这里插入图片描述可能是软件基于不同类型的输入数据,算法不一样,facet 对 case4 技术重复的结果就比较正常了

    在这里插入图片描述

    展开全文
  • 今天我们学习一个拷贝数变异的整合软件——GISTIC2。注意,这和软件本身并不做CNV calling,而是主要用于检测一组样品中显着扩增或缺失的基因组区域(明白一点说就是你需要提供一批样本中的每个样本的CNV检测结果,...

    今天我们学习一个拷贝数变异的整合软件——GISTIC2。注意,这和软件本身并不做CNV calling,而是主要用于检测一组样品中显着扩增或缺失的基因组区域(明白一点说就是你需要提供一批样本中的每个样本的CNV检测结果,这个软件经过呼啦呼啦显著性计算会告诉你这一批样本中显著扩增和缺失的是哪些区域)。这个是癌症基因组CNV分析中十分常见也十分必要的内容。

     

    1.软件安装

    注意:

    a.软件包没有打包在一个文件夹下,所以第2步新建了一个GISTIC2文件夹,请在该文件夹下解压源文件;

    b.第5、6步安装Matlab环境的时候,-destinationFolder后面接的是绝对路径,请将“PATH”替换为GISTIC2文件夹所在的绝对路径。

    c.下面步骤编译好后/PATH/GISTIC2/gistic2即为可执行文件。

    wget -c ftp://ftp.broadinstitute.org/pub/GISTIC2.0/GISTIC_2_0_23.tar.gz
    mkdir GISTIC2 && mv GISTIC_2_0_23.tar.gz GISTIC2/ && cd GISTIC2/
    tar -xzvf GISTIC_2_0_23.tar.gz
    cd MCR_Installer/ && unzip MCRInstaller.zip
    ./install -mode silent -agreeToLicense yes -destinationFolder /PATH/GISTIC2/MATLAB_Compiler_Runtime/
    export XAPPLRESDIR=/PATH/GISTIC2/MATLAB_Compiler_Runtime/v83/X11/app-defaults:$XAPPLRESDIR

     

    2.输入文件准备

    a.Segmentation File

    文件应含有所有样本的CNV信息,共6列,以tab键分割:

    (1)  Sample:样本名称

    (2)  Chromosome:染色体名称

    (3)  Start Position:起始位置

    (4)  End Position :终止位置

    (5)  Num markers:在区域内的标记数目

    (6)  Seg.CN:拷贝数取log2()-1

    文件示例:

    b.Reference Genome File

    参考基因组文件,软件包里的refgenefiles文件夹下提供了不同版本的参考基因组mat文件,这里我是用的hg19参考基因组,因此选择hg19.mat文件作为输入的参考基因组文件。

     

    3.软件使用

    软件的参数有很多,详细可以见软件包下的document,这里介绍几个常用的:

    -b  输出目录,后面接一个输出文件前缀prefix

    -seg  上面讲到的Segmentation文件

    -refgene  上面讲到的参考基因组文件

    -conf  置信水平,默认是0.75

    gistic2  -b /outdir/prefix -seg Segmentation_File -refgene hg19.mat -conf 0.9

     

    4.结果文件解读

    运行结束后,会在输出目录下生成一些文件和图片:

    a.all_lesions.conf_XX.txt(XX是置信度)

    该文件是gistic全部结果的整合,包含了扩增和缺失区域,以及每个区域中扩增或缺失来源于哪些样本。从第10列开始往后是单样本信息,下面解释前9列:

    (1)  Unique Name:  鉴定出的扩增或缺失区域名称;

    (2)  Descriptor:  位于基因组染色体臂的位置;

    (3)  Wide Peak Limits:  最可能包含目标基因的参考基因组边界范围;

    (4)  Peak Limits:  最大扩增或缺失区域的参考基因组边界范围;

    (5)  Region Limits:  显著扩增或缺失区域的参考基因组边界范围;

    (6)  q-values:  峰值区域的q值;

    (7)  Residual q-values:  去除与相同染色体中其他更显著的峰区域重叠的扩增或缺失后,峰区域的q值;

    (8) Broad or Focal:  根据置信度将事件分为“broad”和“focal”;

    (9) Amplitude Threshold:  文件10列及以后列的每个样本给出值的含义解释;

    b.amp_genes.conf_XX.txt(XX是置信度)

    该文件是扩增区域的结果,每一列是一个扩增峰的结果,文件共4行:

    (1)  cytoband

    (2)  q-value

    (3)  residual q-value

    (4)  wide peak boundaries

    c.del_genes.conf_XX.txt(XX是置信度)

    该文件是缺失区域的结果,每一列是一个缺失峰的结果,文件共4行,格式和Amplification Genes File相同。

    d.Gistic Scores File

    GISTIC打分结果文件,文件列出了整个基因组中扩增和缺失的q值[表示为-log10(q)],G得分,异常样本之间的平均幅度以及畸变频率。文件可以导入到IGV中。

    e.raw_copy_number.pdf

    一个pdf热图文件, 基因组沿垂直轴表示,样品水平排列,红色为拷贝数扩增,蓝色为拷贝数缺失。

    f.amp_qplot.pdf、del_qplot.pdf

    扩增或缺失区域单独的展示。

    更多生信小知识关注:

    展开全文

空空如也

空空如也

1 2 3 4 5
收藏数 85
精华内容 34
关键字:

拷贝数变异分析