继上篇：
《肿瘤TMB的计算原理及数学模型》
，本篇主讲“shell脚本计算流程”。
假设我们有肿瘤测序数据SRR2496734，正常组织测序数据SRR2496734，这两组数据是MSKCC发布的，可以从NCBI SRA数据库下载。那么可以用下面的shell脚本计算得到SNP、INDEL，通过简单python脚本提取somatic位点并计算即可得到TMB值。这shell脚本结果跟www.genekras.com平台免费VarScanSomatic计算项目一致。
免费说明：我们参考broad 的做法：
At this time we are able to offer two services for running WDL workflows on Google Cloud using the Cromwell execution engine and the Google Pipelines API. Note that while access to both of these services is free, compute and storage are subject to Google Cloud pricing, charged directly by Google (Broad does not apply any markup).
计算脚本:
redup="true"
nt="6"
ref="hg19.fa"
bed="msk_impact.bed"
normal_R1_fq="SRR2496732_1.fastq.gz"
normal_R2_fq="SRR2496732_2.fastq.gz"
tumor_R1_fq="SRR2496734_1.fastq.gz"
tumor_R2_fq="SRR2496734_2.fastq.gz"
#[1]
bwa mem -t $nt -k 32 -M -R '@RG\tID:line1\tSM:normal\tLB:lib01\tPL:illumina' $ref $normal_R1_fq $normal_R2_fq |samtools view -b -S -o normal.sam
#
bwa mem -t $nt -k 32 -M -R '@RG\tID:line1\tSM:tumor\tLB:lib02\tPL:illumina' $ref $tumor_R1_fq $tumor_R2_fq |samtools view -b -S -o tumor.sam
#[2]
samtools sort -@ $nt normal.sam -o normal.bam -O bam
samtools sort -@ $nt tumor.sam -o tumor.bam -O bam
samtools index -@ $nt normal.bam
samtools index -@ $nt tumor.bam
##-----
if [ "$redup" == "true" ]; then
samtools rmdup normal.bam normalR.bam
samtools rmdup tumor.bam tumorR.bam
samtools index -@ $nt normalR.bam
samtools index -@ $nt tumorR.bam
next_normal="normalR.bam"
next_tumor="tumorR.bam"
else
next_normal="normal.bam"
next_tumor="tumor.bam"
fi
#[3]
samtools mpileup -B -Q 30 -C 50 -q 20 -d 50000 -l $bed -f $ref $next_normal $next_tumor > normal_tumor.mpileup
#[4]
java -jar VarScan.v2.4.2.jar somatic normal_tumor.mpileup --mpileup 1 --output-snp normal_tumor.snp --output-indel normal_tumor.indel --min-var-freq 0.008 --min-coverage 10 --min-coverage-normal 10 --min-coverage-tumor 8 --min-freq-for-hom 0.75 --normal-purity 1.0 --tumor-purity 1.0 --p-value 0.99 --somatic-p-value 0.001 --output-vcf 1 --strand-filter 1
要运行上述代码应先安装BWA、samtools、VarScan2软件，并对参考序列hg19.fa建索引。
流程说明：
[1] 把测序数据比对到参考序列，并通过samtools转换成bam格式文件。
[2] 因为测序文件reads是随机的，让bam文件reads按染色体位置排序。
[3] 生成碱基堆积文件作为VarScan2的输入文件。目前多数软件直接从bam文件读取并输出VCF格式的变异文件。
[4] VarScan2按设定的参数输出vcf格式的变异文件，本质上讲VarScan2仅是格式转换软件，把pileup格式转换成vcf格式，并采用一些过滤参数进行过滤。
也可把上述脚本封装成图形界面计算模式，把参数统一在一个页面上。下图是www.genekras.com上的计算界面，计算结果完全一致。
作者：刘远东：广州克拉斯基因科技有限公司创始人，算法研发工程师，有多年生信算法开发经验。广州克拉斯基因科技有限公司专注肿瘤云计算，国内首家提供UMI标签测序数据在线计算的系统平台!扫描下方二维码可随时与作者全方面交流。零成本生信计算助您成功走在同行前列！点击进行计算http://www.genekras.com。
如果你对本系列有兴趣可以添加助手拉你进群：
