COSINE使用指南_0.1

一、概述

肿瘤DNA序列的亚克隆重建已成为肿瘤进化研究中的重要组成部分，为研究变异与突变过程的相对顺序提供了新的思路。在以往的研究中，多数研究是通过复现前人研究结果，在此基础上运行新的数据集，以验证分析软件的泛化能力，精确度和灵敏度，等等。然而，近年来逐渐增加的肿瘤进化分析软件，也为医生和研究人员的选择带来困难。如何选择最有效的一个或几个分析软件，如何评价这些软件的优劣，成为亟待解决的热门问题。为评估不断增长的肿瘤进化分析软件、给业内人士提供可靠的、清晰的数据分析信息，我们扩展了前人的工作，将12个肿瘤进化分析软件集成在同一平台内，并为12种软件提供了相应的输入数据的生成方式。使用者可以通过原始数据，即fastq文件，生成12种分析软件所需要的输入文件。此项工作旨在保证软件评估中数据来源相同，评价标准相同。在软件使用中，使用者可以通过参考12种软件的输出结果，对同一份样本测序文件进行更为全面的分析。

二、输入数据

使用者可以向服务器上传待处理的数据文件，包括fastq、bam、vcf等文件。也可以使用服务器提供的数据文件用作测试。另外，服务器内包含有hg38，hg19，mm10三种全基因组参考序列，以及相应的已知变异集。无需使用者再次上传。

Fastq文件：

使用自备的待处理的原始测序数据或在公开数据库中下载可以获得，可上传至服务器，在运算平台中可通过gatk全基因组数据分析对其进行处理，获得处理后的bam文件，vcf文件等。

Bam文件：

通过已经上传的fastq文件，处理后生成bam文件。或通过互联网中公开的数据集也包含处理后的bam文件。在运算平台中可以通过变异检测、进化分析软件等工具对其进行处理，进一步获得CNV、SV、SNV和重构进化树等信息。

Vcf文件：

对bam文件使用varscan2、mutect2或者strelka进行变异检测。Vcf文件一般记录了基因组中单位点变异的详细信息。该文件一般作为计算平台中肿瘤进化分析的软件的输入数据。

其他文件：

在计算平台中，可以对其他类型的文件进行处理，主要包含12种分析软件对应的输入文件，或每种软件的中间生成文件，其格式包括但不限于txt、csv和sam等等。

三、预处理

第一步：测序与参考序列的比对并排序

在将测序序列与参考序列的对比过程中，经常会使用到的软件是BWA软件。BWA是一种能够将差异度较小的序列比对到一个较大的参考基因组上的软件包。我们这里使用的是BWA-MEM算法。

对于一般的双端测序数据需要使用两个fastq测序文件和一个fasta参考序列文件进行比对。比对结果以字符的形式保存在sam格式文件中，由于sam文件过大体积，给计算机带来更多压力，一般的，需要使用samtools将其转换为bam二进制文件

bwa mem -t 24 -M -Y -R \ “@RG\tID:${sampleID}\tPL:illumina\tLB:WGS\tSM:${sampleID}” \ $reference $fasgq1 $fastq2 > ${sampleID}.WGS.sam

运行结果文件，输出sam文件的片段。

文本 描述已自动生成

通过BWA文件获得的比对结果需要通过samtools进行进一步处理，samtools是一个用于操作sam和bam格式文件的程序。它可以对sam文件进行格式的导入导出，排序，合并和索引等。

首先使用samtools的view功能，对sam文件格式进行二进制转换，压缩其大小，然后使用sort功能对测序比对结果进行排序，-@ 设定程序并行参数，-o 设定输出文件名称

samtools view -Sb ${sampleID}.WGS.sam > ${sampleID}.WGS.bam samtools sort -@ 24 -o ${sampleID}.WGS.sorted.bam ${sampleID}.WGS.bam

第二步：MarkDuplicates，标记并去除重复序列

在测序过程中，由于PCR扩增过程产生的误差，会导致一些读段被重复读取，进一步导致测序数据的错误。为保证读数的精确性，需要对数据的重复序列进行标记，所使用的软件是GATK。GATK(The Genome Analysis Toolkit)是用于第二代测序数据分析的一款软件，是基因分析的工具集。

使用gatk软件中MarkDuplication功能，可以标记测序数据中由PCR扩增等原因导致的重复序列。具体命令中，-I参数控制输入文件，-O参数控制输出文件，-REMOVE_DUPLICATES参数控制是否将重复序列删除。

gatk MarkDuplicates -I ${sampleID}.WGS.sorted.bam \ -O ${sampleID}.WGS.sorted.markdup.bam \ -M ${sampleID}.WGS.sorted.markdup.txt \ -REMOVE_DUPLICATES true

第三步：BQSR

在使用软件检测的变异位点，即与参考基因不同的部分。通常是由于原始数据中存在一些测序仪器产生的系统性误差，进而导致假阳性变异位点。这一步骤的主要目的是分辨真实的变异位点和不可信的变异位点。

在gatk的BQSR过程分为两步，首先，使用已知变异信息过滤测序结果数据，生成table格式的中间文件，再使用table文件与测序数据，参考序列合并，输出最终去除不可信变异位点的测序数据。

gatk BaseRecalibrator -R $reference \ -I ${sampleID}.WGS.sorted.markdup.bam \ –known-sites $known_indel1 –known-sites $known_snp \ –known-sites $known_indel2 -O ${sampleID}.table   gatk ApplyBQSR –bqsr-recal-file ${sampleID}.table \ -R $reference -I ${sampleID}.WGS.sorted.markdup.bam \ -O ${sampleID}.WGS.sorted.markdup.bqsr.bam

第四步：变异检测-单位点变异（doing）

使用软件：mutect2

mutect2: gatk Mutect2 -R $reference -I ${sampleID}.WGS.sorted.markdup.bqsr.bam \ -L $interval_list -O ${sampleID}.WGS.mutect2.vcf gatk FilterMutectCalls -V ${sampleID}.WGS.mutect2.vcf \ -O ${sampleID}.WGS.somatic.vcf -R $reference  

运行结果文件：内容截图（截图用vcf文件片段），

该步骤的意义：（简单介绍）

第五步：变异检测-拷贝数变异及其他（doing）

使用软件：varscan2 或 ？

Varscan2: varscan somatic < (samtools mpileup \ –no-BAQ -f $reference \ $normal_contral.sorted.markdup.bam \ ${sampleID}.WGS.sorted.markdup.bam) \ $work_dir –mpileup 1 –output-vcf –min-coverage-tumor 15 –min-coverage-norml 12 –somatic-p-value 0.01 varscan processSomatic ${sampleID}.WGS.snp.vcf varscan processSomatic ${sampleID}.WGS.indel.vcf

第六步：分析文件格式转换。

使用mutect2或varscan对处理后的原始数据变异进行检测，结果文件以vcf的形式输出，但集成的12种进化分析软件所使用的的输入数据并不相同。因此，需要使用脚本文件对结果文件进行转换。目前，我们提供sclust、pyclone、FastClone三种进化分析软件输入的格式转换脚本，随后会更新其他脚本转换。

通过对输入文件的转换，完成了变异检测与进化分析的对接，完善了肿瘤细胞信息的分析流程。对使用相同数据对不同软件性能进行评估，提供了可靠的、稳定的保障。

varscan2输出文件转换至sclust输入文件：-i参数输入varscan软件的indel检测结果；-s输入varscan软件的snp检测结果；-o参数输入转换结果文件名称

python $map_dir/bash/varscanTosclust.py -i $varscan2.out.indel.vcf \                                             -s $varscan2.out.snp.vcf \ -o ${sclust}.input.vcf

Mutect2输出文件转换至sclust输入文件：-i 参数输入mutect2软件的indel检测结果；-o参数输入转换结果文件名称。

python $map_dir/bash/varscanTosclust.py -i $varscan2.out.indel.vcf \ -o ${sclust}.input.vcf

sclust输出文件转换至pyclone输入文件：-i 参数输入等位基因信息文件，可由sclust软件获得；-n 输入采样名字符串，控制输出文件名称；-v参数输入用于sclust的vcf信息，与前两个脚本的输出文件相同。

python $map_dir/bash/sclust_to_pyclone.py -i ${sclust}.allelic.txt \                                               -n $sample_name \                                               -v ${sclust}.input.vcf \ -o $output path of Directory

sclust输出文件转换至FastClone输入文件：同sclust输出文件转换至pyclone输入文件。

python $map_dir/bash/ sclust_to_fastclone.py -i ${sclust}.allelic.txt \                                               -n $sample_name \                                               -v ${sclust}.input.vcf \ -o $output path of Directory

运行结果文件：

1、输入文件：varscan的结果文件。

2、输出文件：转换后的文件。

四、肿瘤进化分析

Sclust进化分析

  1. bamprocess:

该步骤负责分析测序数据的读深信息，等位基因信息，GC含量以及测序数据每个未知的详细计数。分别存储在${sampleID}_rcount.txt和${sampleID}_snps.txt两个文件中。特别的，该步骤需要分别对每条染色体进行处理，最后使用-i和-o参数合并所有中间文件，得出最终结果。

Sclust bamprocess -t ${sampleID}.WGS.sorted.markdup.bqsr.bam \ -n ${sampleID}.WGS.sorted.markdup.bqsr.bam \ -o ${sampleID} -part 2 -build hg38 -r chr1 Sclust bamprocess -t ${sampleID}.WGS.sorted.markdup.bqsr.bam \ -n ${sampleID}.WGS.sorted.markdup.bqsr.bam \ -o ${sampleID} -part 2 -build hg38 -r chr2 … … Sclust bamprocess -t ${sampleID}.WGS.sorted.markdup.bqsr.bam \ -n ${sampleID}.WGS.sorted.markdup.bqsr.bam \ -o ${sampleID} -part 2 -build hg38 -r chrY Sclust bamprocess -i ${sampleID}. -o ${sampleID}.
  • cn

该步骤负责计算整个基因组中拷贝数变异信息，亚克隆拷贝数的分数，倍性，纯度和聚类所使用的输入信息。使用6个文件存储。

Sclust cn -rc ${sampleID}_rcount.txt -snp ${sampleID}_snps.txt \ -vcf ${sclust}.input.vcf -o ${sampleID} –ns 1000
  • clust

命令行命令：（就用咱们自己用的参数，注意代码颜色最好保留，直接拷贝）

Pyclone进化分析

  1. analysis pipeline

运行Pyclone全部功能，得出肿瘤细胞进化树。

  • setup and run analysis

Pyclone通过setup_analysis功能指定了输入样本和结果文件输出的位置，同时生成yaml格式的配置文件。

PyClone setup_analysis –in_files ${sampleID}/*.tsv –working_dir ${sampleID}_out_dir

run_analysis功能指定了yaml格式的配置文件输出的位置，同时该命令会输出yaml格式文件。

PyClone run_analysis –config_file ${sampleID}_out_dir/config.yaml
  • plot clusters

plot_clusters功能生成了聚类后的图片，以pdf格式输出，并指定了pdf文件的位置。

PyClone plot_clusters –config_file ${sampleID}_out_dir/config.yaml \ –plot_file ${sampleID}_out_dir/cluster –plot_type density
  • plot loci

plot_clusters功能生成了生成loci的图片文件，并以pdf格式输出

PyClone plot_loci –config_file ${sampleID}_out_dir/config.yaml \ –plot_file ${sampleID}_out_dir/cluster –plot_type density

Fastclone进化分析

  1. run fastclone

fastclone包含三个关键参数，load-pyclon、load-pyclone_truth和load-mutect_battenberg。使用load-python参数，代表载入的文件是跟pyclone的输入文件相同的。load-pyclone_truth参数代表载入的文件是从pyclone输入文件中挑选过的。load-mutect_battenberg代表载入的文件是MuTect和Battenberg方法得到的文件。其具体的运行代码如下

for sample_file in ${sampleID}_dir do (fastclone load-pyclone prop ${sample_file} None solve \ ${sample_file}_output) >> output_record_update 2>&1 done

Author: Y. Cun

Computational biologist