715 lines
30 KiB
Bash
715 lines
30 KiB
Bash
#!/usr/bin/bash
|
||
|
||
################################################################
|
||
## step1: recording all version informations ##
|
||
################################################################
|
||
|
||
################# pipeline_version ########################
|
||
# pipeline: ${pipeline_version}
|
||
|
||
#################### soft_version ########################
|
||
# fastp: ${fastp_version}
|
||
# bwa: ${bwa_version}
|
||
# samtools: ${samtools_version}
|
||
# sentieon: ${sentieon_version}
|
||
# bcftools: ${bcftools_version}
|
||
# annovar: ${annovar_version}
|
||
# vep: ${vep_version}
|
||
# msisensor2: ${msisensor2_version}
|
||
# genefuse: ${genefuse_version}
|
||
# cnvkit: ${cnvkit_version}
|
||
# bamdst: ${bamdst_version}
|
||
# seqtk: ${seqtk_version}
|
||
# htslib: ${htslib_version}
|
||
# python2: ${python2_version}
|
||
# python3: ${python3_version}
|
||
# delly: ${delly_version}
|
||
|
||
#################### bed_version ########################
|
||
# bed: ${bed_version}
|
||
|
||
################# database_version ########################
|
||
# reference: ${reference_version}
|
||
# refgene: ${refgene_version}
|
||
# cytoBand: ${cytoBand_version}
|
||
# 1000g: ${g1000_version}
|
||
# esp: ${esp_version}
|
||
# exac: ${exac_version}
|
||
# cosmic: ${cosmic_version}
|
||
# clinvar: ${clinvar_version}
|
||
# intervar: ${intervar_version}
|
||
# dbnsfp: ${dbnsfp_version}
|
||
# avsnp: ${avsnp_version}
|
||
# private_database: ${private_database}
|
||
# pon_database: ${pon_database}
|
||
|
||
|
||
|
||
################################################################
|
||
## step2: pipeline for analysis ##
|
||
################################################################
|
||
|
||
# Update with the location of the Sentieon software package and license file
|
||
export SENTIEON_LICENSE=${SENTIEON_LICENSE}
|
||
|
||
|
||
threads=${THREAD}
|
||
|
||
cd $outDir
|
||
echo "$sampleID -- beginning at: $(date)"
|
||
|
||
$python3 ${checkUMI} -f $fq1 -s ${sampleType}
|
||
################# fastp ########################
|
||
# fastp is used for remove lower quality bases
|
||
echo "fastp"
|
||
echo "$(date)"
|
||
$fastp -w ${threads} -i $fq1 -o ${sampleID}.R1_001.trim.fq.gz -I $fq2 -O ${sampleID}.R2_001.trim.fq.gz -l 35 -c --detect_adapter_for_pe -j ${sampleID}.jion -h ${sampleID}.html
|
||
|
||
|
||
################# UMI + alignment + dedup ########################
|
||
# extract umi info and alignment to reference
|
||
if [[ ${bedType} =~ ^(${umi_panel})$ ]];
|
||
then
|
||
# remove UMI from following panels
|
||
echo $'\nsentieon umi extract'
|
||
echo "$(date)"
|
||
$sentieon umi extract 21S+T,12M11S+T ${sampleID}.R1_001.trim.fq.gz ${sampleID}.R2_001.trim.fq.gz | \
|
||
( $sentieon bwa mem -t ${threads} -p -C -M -R "@RG\\tID:${sampleID}\\tSM:${sampleID}\\tPL:ILLUMINA\\tLB:lib" -K 10000000 ${hg19fa} - || echo -n 'error') | \
|
||
$sentieon util sort -t ${threads} --sam2bam -i - -o ${sampleID}.sorted.bam
|
||
|
||
echo $'\nsentieon umi consensus'
|
||
echo "$(date)"
|
||
$sentieon umi consensus -i ${sampleID}.sorted.bam --input_format BAM -o ${sampleID}.umi_consensus.fq.gz
|
||
|
||
echo $'\nbwa'
|
||
echo "$(date)"
|
||
( $sentieon bwa mem -t ${threads} -p -C -M -R "@RG\\tID:${sampleID}\\tSM:${sampleID}\\tPL:ILLUMINA\\tLB:lib" -K 10000000 ${hg19fa} ${sampleID}.umi_consensus.fq.gz || echo -n 'error' ) | \
|
||
$sentieon util sort --umi_post_process --sam2bam -i - -o ${sampleID}.sorted.dedup.bam
|
||
|
||
elif [[ ${bedType} =~ (${probe_panel}) && ${sampleType} =~ (ctdna) ]];
|
||
then
|
||
# remove UMI from following panels
|
||
echo $'\nsentieon umi'
|
||
echo "$(date)"
|
||
$sentieon umi extract -d 6M1S+T,6M1S+T ${sampleID}.R1_001.trim.fq.gz ${sampleID}.R2_001.trim.fq.gz | \
|
||
( $sentieon bwa mem -t ${threads} -p -C -M -R "@RG\\tID:${sampleID}\\tSM:${sampleID}\\tPL:ILLUMINA\\tLB:lib" -K 10000000 ${hg19fa} - || echo -n 'error') | \
|
||
$sentieon util sort -t ${threads} --sam2bam -i - -o ${sampleID}.sorted.bam
|
||
|
||
echo $'\nsentieon umi consensus'
|
||
echo "$(date)"
|
||
$sentieon umi consensus -i ${sampleID}.sorted.bam --input_format BAM -o ${sampleID}.umi_consensus.fq.gz
|
||
|
||
echo $'\nbwa'
|
||
echo "$(date)"
|
||
( $sentieon bwa mem -t ${threads} -p -C -M -R "@RG\\tID:${sampleID}\\tSM:${sampleID}\\tPL:ILLUMINA\\tLB:lib" -K 10000000 ${hg19fa} ${sampleID}.umi_consensus.fq.gz || echo -n 'error' ) | \
|
||
$sentieon util sort --umi_post_process --sam2bam -i - -o ${sampleID}.sorted.dedup.bam
|
||
|
||
else
|
||
# following panels does not contain UMI
|
||
echo $'\nbwa'
|
||
echo "$(date)"
|
||
( $sentieon bwa mem -t ${threads} -M -R "@RG\\tID:${sampleID}\\tSM:${sampleID}\\tPL:ILLUMINA\\tLB:lib" -K 10000000 ${hg19fa} ${sampleID}.R1_001.trim.fq.gz ${sampleID}.R2_001.trim.fq.gz || echo -n 'error' ) | \
|
||
$sentieon util sort -o ${sampleID}.sorted.bam -t ${threads} --sam2bam -i -
|
||
|
||
if [[ ${bedType} =~ ^(${probe_panel})$ ]]; then
|
||
echo $'\nRemove Duplicate'
|
||
echo "$(date)"
|
||
$sentieon driver -t ${threads} -i ${sampleID}.sorted.bam --algo LocusCollector --fun score_info ${sampleID}.score.txt
|
||
$sentieon driver -t ${threads} -i ${sampleID}.sorted.bam --algo Dedup --score_info ${sampleID}.score.txt --metrics ${sampleID}.dedup_metrics.txt ${sampleID}.sorted.dedup.bam
|
||
rm ${sampleID}.score.txt* ${sampleID}.dedup_metrics.txt
|
||
fi
|
||
fi
|
||
|
||
|
||
################## Variant Calling ######################
|
||
# TNscope is used for somatic SNVs & small INDELs calling
|
||
if [[ ${bedType} =~ ^(${pcr_panel})$ ]]; then
|
||
|
||
echo $'\nHaplotyper'
|
||
echo "$(date)"
|
||
$sentieon driver -r ${hg19fa} -t ${threads} -i ${sampleID}.sorted.bam \
|
||
--interval ${bedPath}/${bedType}.bed \
|
||
--algo Haplotyper --emit_conf=10 --call_conf=10 \
|
||
${sampleID}.germline.raw.vcf
|
||
|
||
echo $'\nTNscope'
|
||
echo "$(date)"
|
||
|
||
if [[ ${bedType} =~ (AT|CA) ]];
|
||
then
|
||
$sentieon driver -r ${hg19fa} -t ${threads} -i ${sampleID}.sorted.bam \
|
||
--interval ${bedPath}/${bedType}.bed --interval_padding 10 \
|
||
--algo TNscope \
|
||
--tumor_sample ${sampleID} \
|
||
--dbsnp ${dbsnp138vcf} \
|
||
--disable_detector sv \
|
||
--max_fisher_pv_active 0.05 \
|
||
--min_tumor_allele_frac 0.009 \
|
||
--sv_mask_ext 10 \
|
||
--prune_factor 0 \
|
||
--assemble_mode 2 \
|
||
--resample_depth 100000 \
|
||
--pon ${pon} \
|
||
${sampleID}.output_tnscope.pre_filter.vcf
|
||
else
|
||
$sentieon driver -r ${hg19fa} -t ${threads} -i ${sampleID}.sorted.bam \
|
||
--interval ${bedPath}/${bedType}.bed --interval_padding 10 \
|
||
--algo TNscope \
|
||
--tumor_sample ${sampleID} \
|
||
--dbsnp ${dbsnp138vcf} \
|
||
--disable_detector sv \
|
||
--max_fisher_pv_active 0.05 \
|
||
--min_tumor_allele_frac 0.009 \
|
||
--sv_mask_ext 10 \
|
||
--prune_factor 0 \
|
||
--assemble_mode 2 \
|
||
--resample_depth 100000 \
|
||
--trim_primer ${bedPath}/${bedType}.amplicon.bed \
|
||
--pon ${pon} \
|
||
${sampleID}.output_tnscope.pre_filter.vcf
|
||
fi
|
||
|
||
echo $'\nbcftools filter'
|
||
echo "$(date)"
|
||
$bcftools annotate -x "FILTER/triallelic_site" ${sampleID}.output_tnscope.pre_filter.vcf | \
|
||
$bcftools filter -m + -s "insignificant" -e "(PV>0.25 && PV2>0.25)" | \
|
||
$bcftools filter -m + -s "insignificant" -e "(INFO/STR == 1 && PV>0.05)" | \
|
||
$bcftools filter -m + -s "low_qual" -e "QUAL < 20" | \
|
||
$bcftools filter -m + -s "short_tandem_repeat" -e "RPA[0]>=10" | \
|
||
$bcftools filter -m + -s "noisy_region" -e "ECNT>5" | \
|
||
$sentieon util vcfconvert - ${sampleID}.somatic.raw.vcf
|
||
rm ${sampleID}.output_tnscope.pre_filter.vcf*
|
||
|
||
elif [[ ${bedType} =~ ^(${probe_panel})$ && ${sampleType} != "ctdna" ]]; then
|
||
|
||
echo $'\nHaplotyper'
|
||
echo "$(date)"
|
||
$sentieon driver -r ${hg19fa} -t ${threads} -i ${sampleID}.sorted.dedup.bam \
|
||
--interval ${bedPath}/${bedType}.bed \
|
||
--algo Haplotyper --emit_conf=10 --call_conf=10 \
|
||
${sampleID}.germline.raw.vcf
|
||
|
||
echo $'\nTNscope'
|
||
echo "$(date)"
|
||
$sentieon driver -r ${hg19fa} -t ${threads} -i ${sampleID}.sorted.dedup.bam \
|
||
--interval ${bedPath}/${bedType}.bed --interval_padding 10 \
|
||
--algo TNscope \
|
||
--tumor_sample ${sampleID} \
|
||
--dbsnp ${dbsnp138vcf} \
|
||
--disable_detector sv \
|
||
--max_fisher_pv_active 0.05 \
|
||
--min_tumor_allele_frac 0.005 \
|
||
--filter_t_alt_frac 0.005 \
|
||
--resample_depth 100000 \
|
||
--assemble_mode 4 \
|
||
--pon ${pon} \
|
||
${sampleID}.output_tnscope.pre_filter.vcf
|
||
|
||
echo $'\nbcftools filter'
|
||
echo "$(date)"
|
||
echo $'\nbcftools filter'
|
||
echo "$(date)"
|
||
$bcftools annotate -x "FILTER/triallelic_site" ${sampleID}.output_tnscope.pre_filter.vcf | \
|
||
$bcftools filter -m + -s "insignificant" -e "(PV>0.25 && PV2>0.25)" | \
|
||
$bcftools filter -m + -s "insignificant" -e "(INFO/STR == 1 && PV>0.05)" | \
|
||
$bcftools filter -m + -s "orientation_bias" -e "FMT/FOXOG[0] == 1" | \
|
||
$bcftools filter -m + -s "strand_bias" -e "SOR > 3" | \
|
||
$bcftools filter -m + -s "low_qual" -e "QUAL < 20" | \
|
||
$bcftools filter -m + -s "short_tandem_repeat" -e "RPA[0]>=10" | \
|
||
$bcftools filter -m + -s "noisy_region" -e "ECNT>5" | \
|
||
$bcftools filter -m + -s "read_pos_bias" -e "FMT/ReadPosRankSumPS[0] < -8" | \
|
||
$sentieon util vcfconvert - ${sampleID}.somatic.raw.vcf
|
||
rm ${sampleID}.output_tnscope.pre_filter.vcf*
|
||
|
||
else
|
||
echo $'\nHaplotyper'
|
||
echo "$(date)"
|
||
$sentieon driver -r ${hg19fa} -t ${threads} -i ${sampleID}.sorted.dedup.bam \
|
||
--interval ${bedPath}/${bedType}.bed \
|
||
--algo Haplotyper --emit_conf=10 --call_conf=10 \
|
||
${sampleID}.germline.raw.vcf
|
||
|
||
echo $'\nTNscope'
|
||
echo "$(date)"
|
||
$sentieon driver -r ${hg19fa} -t ${threads} -i ${sampleID}.sorted.dedup.bam --interval ${bedPath}/${bedType}.bed \
|
||
--algo TNscope \
|
||
--tumor_sample ${sampleID} \
|
||
--dbsnp ${dbsnp138vcf} \
|
||
--pcr_indel_model NONE \
|
||
--disable_detector sv \
|
||
--min_tumor_allele_frac 0.0003 \
|
||
--filter_t_alt_frac 0.0003 \
|
||
--max_error_per_read 3 \
|
||
--min_init_tumor_lod 3.0 \
|
||
--min_tumor_lod 3.0 \
|
||
--min_base_qual 40 \
|
||
--min_base_qual_asm 40 \
|
||
--resample_depth 100000 \
|
||
--assemble_mode 4 \
|
||
--pon ${pon} \
|
||
${sampleID}.output_tnscope.pre_filter.vcf
|
||
|
||
echo $'\nbcftools filter'
|
||
echo "$(date)"
|
||
$bcftools annotate -x "FILTER/triallelic_site" ${sampleID}.output_tnscope.pre_filter.vcf | \
|
||
$bcftools filter -m + -s "low_qual" -e "QUAL < 10" | \
|
||
$bcftools filter -m + -s "short_tandem_repeat" -e "RPA[0]>=10" | \
|
||
$bcftools filter -m + -s "read_pos_bias" -e "FMT/ReadPosRankSumPS[0] < -5" | \
|
||
$bcftools filter -m + -s "base_qual_bias" -e "FMT/BaseQRankSumPS[0] < -5" | \
|
||
$sentieon util vcfconvert - ${sampleID}.somatic.raw.vcf
|
||
rm ${sampleID}.output_tnscope.pre_filter.vcf*
|
||
fi
|
||
|
||
|
||
echo $'\nfilter repeat regions'
|
||
echo "$(date)"
|
||
grep "#" ${sampleID}.somatic.raw.vcf > ${sampleID}.somatic.filter.vcf
|
||
grep -v "#" ${sampleID}.somatic.raw.vcf | awk '{len=split($5,a,",");if(len<3 || (len>=3 && $7 !~ "panel_of_normals")) print $0}' >> ${sampleID}.somatic.filter.vcf
|
||
|
||
|
||
################# 1p19q ########################
|
||
|
||
if [[ ${bedType} =~ (AT) ]];
|
||
then
|
||
$samtools mpileup -l ${bedPath}/AT.force.bed -f ${hg19fa} -d 30000 -Q 0 -O -s ${sampleID}.sorted.bam > ${sampleID}.mpileup.snp.txt
|
||
$python3 ${dig_pileup} ${sampleID}.mpileup.snp.txt -o ${sampleID}.mpileup.snp.vcf
|
||
$python3 ${p1q19_paint} ${sampleID}.mpileup.snp.vcf
|
||
rm ${sampleID}.mpileup.snp.txt
|
||
fi
|
||
|
||
################# bcftools ########################
|
||
# bcftools is used for retrieving mutations of interest, or mutations with higher confidence
|
||
|
||
echo $'\nbcftools'
|
||
echo "$(date)"
|
||
#germline
|
||
$bcftools norm -m -both ${sampleID}.germline.raw.vcf -o ${sampleID}.germline.raw.split.vcf
|
||
$bcftools norm -f ${hg19fa} ${sampleID}.germline.raw.split.vcf -o ${sampleID}.germline.split.norm.vcf
|
||
rm ${sampleID}.germline.raw.vcf* ${sampleID}.germline.raw.split.vcf
|
||
|
||
# #somatic
|
||
$bcftools norm -m -both ${sampleID}.somatic.filter.vcf -o ${sampleID}.somatic.raw.split.vcf
|
||
$bcftools norm -f ${hg19fa} ${sampleID}.somatic.raw.split.vcf -o ${sampleID}.somatic.split.norm.vcf
|
||
rm ${sampleID}.somatic.filter.vcf ${sampleID}.somatic.raw.split.vcf
|
||
|
||
|
||
###################### ANNOVAR ########################
|
||
# used for annotation
|
||
|
||
###################### germline ########################
|
||
echo $'\ngermline'
|
||
echo "$(date)"
|
||
mkdir anno
|
||
${annovar} ${sampleID}.germline.split.norm.vcf ${humandb} \
|
||
--outfile anno/${sampleID}.germline.anno --buildver hg19 \
|
||
--protocol refGene,cytoBand,1000g2015aug_all,esp6500siv2_all,exac03nonpsych,cosmic70,clinvar_20220624,intervar_20180118,dbnsfp35a,avsnp150 \
|
||
--operation g,r,f,f,f,f,f,f,f,f --vcfinput --remove
|
||
rm anno/${sampleID}.germline.anno.avinput
|
||
|
||
if ([[ "$bedType" =~ (${germline_panel1}) ]]) || ([[ "$bedType" =~ (${germline_panel2}) ]]); then
|
||
$python3 ${germline_filter} anno/${sampleID}.germline.anno.hg19_multianno.vcf -o anno/${sampleID}.germline.inital.filter.vcf > anno/${sampleID}.germline.filter_reason.txt
|
||
|
||
# vep注释（主要为HGVS）
|
||
${vep} --dir ${vep_dir} --cache --offline --cache_version 98 \
|
||
--use_given_ref --refseq --assembly GRCh37 --format vcf \
|
||
--fa ${hg19fa} --force_overwrite --vcf --variant_class --gene_phenotype --vcf_info_field ANN \
|
||
--hgvs --hgvsg --transcript_version \
|
||
-i anno/${sampleID}.germline.inital.filter.vcf -o anno/${sampleID}.germline.vep_anno.vcf
|
||
rm anno/${sampleID}.germline.inital.filter.vcf
|
||
|
||
# HGVS命名注释
|
||
$python3 ${anno_hgvs} anno/${sampleID}.germline.vep_anno.vcf ${clinic_transcripts} ${refFlat} -o anno/${sampleID}.germline.hgvs_anno.vcf
|
||
rm anno/${sampleID}.germline.vep_anno.vcf*
|
||
|
||
# VEP修复后假阳性删除
|
||
$python3 ${snv_second_filter} anno/${sampleID}.germline.hgvs_anno.vcf -k germline -o anno/${sampleID}.germline.final.filter.vcf >>anno/${sampleID}.germline.filter_reason.txt
|
||
fi
|
||
|
||
|
||
###################### somatic ########################
|
||
echo $'\nsomatic'
|
||
echo "$(date)"
|
||
|
||
# 过滤 nopass 突变
|
||
grep "#" ${sampleID}.somatic.split.norm.vcf > anno/${sampleID}.somatic.split.filter.vcf
|
||
grep -v "#" ${sampleID}.somatic.split.norm.vcf | awk '{len=split($7,a,";");split($NF,b,":");if(len<2 || (len>=2 && b[3] >=0.01)) print $0}' >> anno/${sampleID}.somatic.split.filter.vcf
|
||
|
||
# annovar注释
|
||
${annovar} anno/${sampleID}.somatic.split.filter.vcf ${humandb} \
|
||
--outfile anno/${sampleID}.somatic.annovar --buildver hg19 \
|
||
--protocol refGene,1000g2015aug_all,esp6500siv2_all,exac03nonpsych,cosmic70,clinvar_20220624,intervar_20180118,dbnsfp35a,avsnp150,dbscsnv11 \
|
||
--operation g,f,f,f,f,f,f,f,f,f --vcfinput --remove
|
||
|
||
# vep注释(主要为HGVS)
|
||
${vep} --dir ${vep_dir} --cache --offline --cache_version 98 \
|
||
--use_given_ref --refseq --assembly GRCh37 --format vcf \
|
||
--fa ${hg19fa} --force_overwrite --vcf --variant_class --gene_phenotype --vcf_info_field ANN \
|
||
--hgvs --hgvsg --transcript_version \
|
||
-i anno/${sampleID}.somatic.annovar.hg19_multianno.vcf -o anno/${sampleID}.somatic.vep_anno.vcf
|
||
|
||
# HGVS命名注释
|
||
$python3 ${anno_hgvs} anno/${sampleID}.somatic.vep_anno.vcf ${clinic_transcripts} ${refFlat} -o anno/${sampleID}.somatic.hgvs_anno.vcf
|
||
|
||
# 初步过滤(假阳性等)
|
||
$python3 ${snv_init_filter} anno/${sampleID}.somatic.hgvs_anno.vcf \
|
||
-o anno/${sampleID}.somatic.inital.filter.vcf -t ${sampleType} -m ${mode} > anno/${sampleID}.somatic.filter_reason.txt
|
||
|
||
# 合并同一单倍型上的复杂突变
|
||
$python3 ${merge_mnv} anno/${sampleID}.somatic.inital.filter.vcf $hg19fa -o anno/${sampleID}.somatic.merge_mnv.vcf
|
||
|
||
# 将合并后的mnp位置左移校准
|
||
$bcftools norm -m -both -f ${hg19fa} anno/${sampleID}.somatic.merge_mnv.vcf -o anno/${sampleID}.somatic.merged_norm.vcf
|
||
|
||
# 拆分MNV突变
|
||
$python3 ${find_mnp} anno/${sampleID}.somatic.merged_norm.vcf
|
||
|
||
if [ -f "anno/${sampleID}.somatic.mnp.vcf" ];
|
||
then
|
||
# annovar注释MNV
|
||
${annovar} anno/${sampleID}.somatic.mnp.vcf ${humandb} \
|
||
--outfile anno/${sampleID}.somatic.mnp.annovar --buildver hg19 \
|
||
--protocol refGene,1000g2015aug_all,esp6500siv2_all,exac03nonpsych,cosmic70,clinvar_20220624,intervar_20180118,dbnsfp35a,avsnp150,dbscsnv11 \
|
||
--operation g,f,f,f,f,f,f,f,f,f --vcfinput --remove
|
||
|
||
# vep注释MNV(主要为HGVS)
|
||
${vep} --dir ${vep_dir} --cache --offline --cache_version 98 \
|
||
--use_given_ref --refseq --assembly GRCh37 --format vcf \
|
||
--fa ${hg19fa} --force_overwrite --vcf --variant_class --gene_phenotype --vcf_info_field ANN \
|
||
--hgvs --hgvsg --transcript_version \
|
||
-i anno/${sampleID}.somatic.mnp.annovar.hg19_multianno.vcf -o anno/${sampleID}.somatic.mnp.vep_anno.vcf
|
||
|
||
# HGVS命名注释MNV
|
||
$python3 ${anno_hgvs} anno/${sampleID}.somatic.mnp.vep_anno.vcf ${clinic_transcripts} ${refFlat} -o anno/${sampleID}.somatic.mnp.hgvs_anno.vcf
|
||
|
||
# 合并突变
|
||
$bgzip -f anno/${sampleID}.somatic.common.vcf
|
||
$tabix -p vcf anno/${sampleID}.somatic.common.vcf.gz
|
||
|
||
$bgzip -f anno/${sampleID}.somatic.mnp.hgvs_anno.vcf
|
||
$tabix -p vcf anno/${sampleID}.somatic.mnp.hgvs_anno.vcf.gz
|
||
|
||
$bcftools concat anno/${sampleID}.somatic.common.vcf.gz anno/${sampleID}.somatic.mnp.hgvs_anno.vcf.gz -a -d all -o anno/${sampleID}.somatic.merge.vcf
|
||
|
||
else
|
||
mv anno/${sampleID}.somatic.merged_norm.vcf anno/${sampleID}.somatic.merge.vcf
|
||
fi
|
||
|
||
# VEP修复后假阳性及良性变异删除
|
||
$python3 ${snv_second_filter} anno/${sampleID}.somatic.merge.vcf -k somatic -o anno/${sampleID}.somatic.final.filter.vcf >>anno/${sampleID}.somatic.filter_reason.txt
|
||
|
||
# 统计变异删除的原因
|
||
$python3 ${get_filter_reason} . -k somatic
|
||
|
||
|
||
###################### tmb ########################
|
||
if [[ ${bedType} =~ ^(${tmb_panel})$ ]]; then
|
||
echo $'\nTMB'
|
||
echo "$(date)"
|
||
|
||
# 过滤低频假阳性
|
||
grep "#" ${sampleID}.somatic.split.norm.vcf > anno/${sampleID}.tmb.filter_1.vcf
|
||
if [[ ${sampleType} =~ (ctdna) ]]; then
|
||
grep -v "#" ${sampleID}.somatic.split.norm.vcf | awk '{split($NF,a,":");if(a[3]>0.05 && a[3] <= 0.5) print $0}' >> anno/${sampleID}.tmb.filter_1.vcf
|
||
else
|
||
grep -v "#" ${sampleID}.somatic.split.norm.vcf | awk '{split($NF,a,":");if(a[3]>0.05) print $0}' >> anno/${sampleID}.tmb.filter_1.vcf
|
||
fi
|
||
|
||
# 过滤filter信息大于2的突变
|
||
grep "#" anno/${sampleID}.tmb.filter_1.vcf > anno/${sampleID}.tmb.filter_2.vcf
|
||
grep -v "#" anno/${sampleID}.tmb.filter_1.vcf | awk '{len=split($7,a,";");if(len<2) print $0}' >> anno/${sampleID}.tmb.filter_2.vcf
|
||
|
||
# 合并同一单倍型上的复杂突变
|
||
$python3 ${merge_mnv} anno/${sampleID}.tmb.filter_2.vcf $hg19fa -o anno/${sampleID}.tmb.merged.vcf
|
||
|
||
# 将合并后的mnp位置左移校准
|
||
$bcftools norm -m -both -f ${hg19fa} anno/${sampleID}.tmb.merged.vcf -o anno/${sampleID}.tmb.merged_norm.vcf
|
||
rm anno/${sampleID}.tmb.merged.vcf
|
||
|
||
# annovar 注释相关数据库
|
||
${annovar} anno/${sampleID}.tmb.merged_norm.vcf ${humandb} \
|
||
--outfile anno/${sampleID}.tmb.annovar --buildver hg19 \
|
||
--protocol refGene,cytoBand,1000g2015aug_all,esp6500siv2_all,exac03nonpsych,cosmic70,clinvar_20220624,intervar_20180118,dbnsfp35a,avsnp150 \
|
||
--operation g,r,f,f,f,f,f,f,f,f --vcfinput --remove
|
||
rm anno/${sampleID}.tmb.merged_norm.vcf
|
||
|
||
# vep注释（主要为HGVS）
|
||
${vep} --dir ${vep_dir} --cache --offline --cache_version 98 \
|
||
--use_given_ref --refseq --assembly GRCh37 --format vcf \
|
||
--fa ${hg19fa} --force_overwrite --vcf --variant_class --gene_phenotype --vcf_info_field ANN \
|
||
--hgvs --hgvsg --transcript_version \
|
||
-i anno/${sampleID}.tmb.annovar.hg19_multianno.vcf -o anno/${sampleID}.tmb.vep_anno.vcf
|
||
rm anno/${sampleID}.tmb.annovar*
|
||
|
||
# HGVS命名注释
|
||
$python3 ${anno_hgvs} anno/${sampleID}.tmb.vep_anno.vcf ${clinic_transcripts} ${refFlat} -o anno/${sampleID}.tmb.hgvs_anno.vcf
|
||
rm anno/${sampleID}.tmb.vep_anno.vcf*
|
||
|
||
# 过滤其他假阳性，germline等
|
||
$python3 ${tmb_filter} anno/${sampleID}.tmb.hgvs_anno.vcf -t ${sampleType} -o anno/${sampleID}.tmb.final.filter.vcf -m ${mode} > anno/${sampleID}.tmb.filter_reason.txt
|
||
rm anno/${sampleID}.tmb.hgvs_anno.vcf
|
||
|
||
# 统计变异删除的原因
|
||
$python3 ${get_filter_reason} . -k tmb
|
||
rm anno/${sampleID}.tmb.filter_1.vcf anno/${sampleID}.tmb.filter_2.vcf
|
||
fi
|
||
rm ${sampleID}.somatic.split.norm.vcf
|
||
|
||
|
||
###################### concat vcf ########################
|
||
if ([[ "$bedType" =~ (${germline_panel1}) ]]) || ([[ "$bedType" =~ (QR) && ${sampleType} != "blood" && ${sampleType} != "dna" ]]);then
|
||
echo $'\nmerge_vcf'
|
||
echo "$(date)"
|
||
|
||
cp anno/${sampleID}.germline.final.filter.vcf anno/${sampleID}.germline.final.filter.concat.vcf
|
||
$bgzip -f anno/${sampleID}.germline.final.filter.concat.vcf
|
||
$tabix -p vcf anno/${sampleID}.germline.final.filter.concat.vcf.gz
|
||
|
||
cp anno/${sampleID}.somatic.final.filter.vcf anno/${sampleID}.somatic.final.filter.concat.vcf
|
||
$bgzip -f anno/${sampleID}.somatic.final.filter.concat.vcf
|
||
$tabix -p vcf anno/${sampleID}.somatic.final.filter.concat.vcf.gz
|
||
|
||
if [[ "$bedType" =~ (${tmb_panel}) ]];then
|
||
cp anno/${sampleID}.tmb.final.filter.vcf anno/${sampleID}.tmb.final.filter.concat.vcf
|
||
$bgzip -f anno/${sampleID}.tmb.final.filter.concat.vcf
|
||
$tabix -p vcf anno/${sampleID}.tmb.final.filter.concat.vcf.gz
|
||
$bcftools concat anno/${sampleID}.somatic.final.filter.concat.vcf.gz anno/${sampleID}.germline.final.filter.concat.vcf.gz anno/${sampleID}.tmb.final.filter.concat.vcf.gz -a -d all -o anno/${sampleID}.somatic.concat.vcf
|
||
rm anno/${sampleID}.tmb.final.filter.concat.vcf.gz*
|
||
else
|
||
$bcftools concat anno/${sampleID}.somatic.final.filter.concat.vcf.gz anno/${sampleID}.germline.final.filter.concat.vcf.gz -a -d all -o anno/${sampleID}.somatic.concat.vcf
|
||
fi
|
||
rm anno/${sampleID}.somatic.final.filter.concat.vcf.gz* anno/${sampleID}.germline.final.filter.concat.vcf.gz*
|
||
fi
|
||
|
||
###################### anno drugs ########################
|
||
echo $'\nanno_drugs'
|
||
echo "$(date)"
|
||
|
||
# 注释药物相关数据库
|
||
if ([[ "$bedType" =~ (${germline_panel1}) ]]) || ([[ "$bedType" =~ (QR) && ${sampleType} != "blood" && ${sampleType} != "dna" ]]);then
|
||
$python3 ${snv_out} anno/${sampleID}.somatic.concat.vcf -o anno/${sampleID}.somatic.xlsx -t ${sampleType} >> anno/${sampleID}.somatic.filter_reason.txt
|
||
|
||
elif [[ "$bedType" =~ (${germline_panel2}) && "$sampleType" =~ (blood|dna) ]];then
|
||
$python3 ${snv_out} anno/${sampleID}.germline.final.filter.vcf -o anno/${sampleID}.somatic.xlsx -t ${sampleType} >> anno/${sampleID}.somatic.filter_reason.txt
|
||
|
||
else
|
||
$python3 ${snv_out} anno/${sampleID}.somatic.final.filter.vcf -o anno/${sampleID}.somatic.xlsx -t ${sampleType} >> anno/${sampleID}.somatic.filter_reason.txt
|
||
fi
|
||
|
||
###################### add_pipe_tag ######################
|
||
# 添加流程pipe tag
|
||
$python3 ${add_pipe_tag} anno/${sampleID}.somatic.xlsx
|
||
|
||
#################### sanger_determine ####################
|
||
# 确定germline突变是否需要验证
|
||
$python3 ${sanger_determine} anno/${sampleID}.somatic.xlsx
|
||
|
||
#################### cis_trans ####################
|
||
# 确定T790M和G797S共突变的顺式或反式
|
||
$python3 ${cis_trans} anno/${sampleID}.somatic.xlsx
|
||
|
||
###################### core_filter ########################
|
||
if [[ "$bedType" =~ (${core_panel}) ]];then
|
||
$python3 ${core_filter} anno/${sampleID}.somatic.xlsx
|
||
fi
|
||
|
||
###################### qc_bamdst ########################
|
||
echo $'\nqc_bamdst'
|
||
echo "$(date)"
|
||
|
||
if [[ ${bedType} =~ ^(${pcr_panel})$ || ${bedType} =~ ^(${umi_panel})$ ]]; then
|
||
mkdir -p qc/bamdst
|
||
$bamdst --uncover 30 -p ${bedPath}/${bedType}.bed --cutoffdepth 500 -o qc/bamdst ${sampleID}.sorted.bam
|
||
ls qc/bamdst/ | xargs -I{} mv qc/bamdst/{} qc/bamdst/${sampleID}.{}
|
||
|
||
elif [[ ${bedType} =~ ^(${probe_panel})$ && ${sampleType} != "ctdna" ]]; then
|
||
mkdir -p qc/bamdst
|
||
$bamdst --uncover 30 -p ${bedPath}/${bedType}.bed --cutoffdepth 1000 -o qc/bamdst ${sampleID}.sorted.dedup.bam
|
||
ls qc/bamdst/ | xargs -I{} mv qc/bamdst/{} qc/bamdst/${sampleID}.{}
|
||
|
||
else
|
||
mkdir -p qc/bamdst/raw
|
||
if [[ "$sampleType" =~ (ctdna) ]]; then
|
||
$bamdst --uncover 30 -p ${bedPath}/${bedType}.bed --cutoffdepth 1000 -o qc/bamdst/raw ${sampleID}.sorted.bam
|
||
else
|
||
$bamdst --uncover 30 -p ${bedPath}/${bedType}.bed --cutoffdepth 500 -o qc/bamdst/raw ${sampleID}.sorted.bam
|
||
fi
|
||
ls qc/bamdst/raw/ | xargs -I{} mv qc/bamdst/raw/{} qc/bamdst/raw/${sampleID}.raw.{}
|
||
|
||
mkdir -p qc/bamdst/dedup
|
||
if [[ "$sampleType" =~ (ctdna) ]]; then
|
||
$bamdst --uncover 30 -p ${bedPath}/${bedType}.bed --cutoffdepth 1000 -o qc/bamdst/dedup ${sampleID}.sorted.dedup.bam
|
||
else
|
||
$bamdst --uncover 30 -p ${bedPath}/${bedType}.bed --cutoffdepth 500 -o qc/bamdst/dedup ${sampleID}.sorted.dedup.bam
|
||
fi
|
||
ls qc/bamdst/dedup | xargs -I{} mv qc/bamdst/dedup/{} qc/bamdst/dedup/${sampleID}.dedup.{}
|
||
fi
|
||
|
||
|
||
###################### chemical ########################
|
||
|
||
echo $'\nchemical'
|
||
echo "$(date)"
|
||
|
||
if [[ ${bedType} =~ ^(DA)$ ]]; then
|
||
$samtools bedcov ${bedPath}/chemical_DA.bed ${sampleID}.sorted.bam >${sampleID}.chemical.depth.txt
|
||
else
|
||
$samtools bedcov ${bedPath}/chemical_DF.bed ${sampleID}.sorted.bam >${sampleID}.chemical.depth.txt
|
||
fi
|
||
|
||
if [[ "$bedType" =~ (${chemical_panel}) ]]; then
|
||
$python3 ${germline_chemical} anno/${sampleID}.germline.anno.hg19_multianno.txt >anno/${sampleID}.chemical.txt
|
||
$python3 ${check_sample_data} anno/${sampleID}.chemical.txt
|
||
rm ${sampleID}.chemical.depth.txt anno/${sampleID}.chemical.txt
|
||
fi
|
||
|
||
|
||
###################### Copy Number Aberration ######################
|
||
### hybrid capture: DG|DS|DH|DT|DL|DW|DR|DJ|DF|DQ|DE
|
||
### targeted amplicon seq: DP|HH|HT|HW|HR|HC|HF|ZJ|ZG|YW|QJ|YR|YZ|YW|YF
|
||
### 需要做cnv的panel：DF, DG, DJ, DQ, DR, DS, DW, DE
|
||
|
||
##################### CNVkit ########################
|
||
if [[ ${bedType} =~ ^(${cnv_panel1})$ ]];
|
||
then
|
||
echo "CNA starts at -- `date` "
|
||
|
||
if [[ ${bedType} =~ ^(${cnv_panel2})$ ]];
|
||
then
|
||
$python3 $cnvkit coverage ${sampleID}.sorted.dedup.bam ${cnv_target}/${bedType}.target.bed -p ${threads} -o ${sampleID}.targetcoverage.cnn
|
||
$python3 $cnvkit fix ${sampleID}.targetcoverage.cnn ${cnv_ref}/empty.cnn ${cnv_ref}/${bedType}.flat.cnn -o ${sampleID}.cnr
|
||
$python3 $cnvkit segment -p ${threads} ${sampleID}.cnr --drop-low-coverage -o ${sampleID}.cns
|
||
|
||
elif [[ ${bedType} =~ ^(${cnv_panel3})$ ]];
|
||
then
|
||
$python3 $cnvkit coverage ${sampleID}.sorted.bam ${cnv_target}/${bedType}.target.bed -p ${threads} -o ${sampleID}.targetcoverage.cnn
|
||
$python3 $cnvkit fix ${sampleID}.targetcoverage.cnn ${cnv_ref}/empty.cnn ${cnv_ref}/${bedType}.flat.cnn --no-edge -o ${sampleID}.cnr
|
||
$python3 $cnvkit segment -p ${threads} ${sampleID}.cnr --drop-low-coverage -o ${sampleID}.cns
|
||
|
||
else
|
||
$python3 $cnvkit coverage ${sampleID}.sorted.dedup.bam ${cnv_target}/${bedType}.target.bed -p ${threads} -o ${sampleID}.targetcoverage.cnn
|
||
$python3 $cnvkit coverage ${sampleID}.sorted.dedup.bam ${cnv_anti}/${bedType}.antitarget.bed -p ${threads} -o ${sampleID}.antitargetcoverage.cnn
|
||
$python3 $cnvkit fix ${sampleID}.targetcoverage.cnn ${sampleID}.antitargetcoverage.cnn ${cnv_ref}/${bedType}.flat.cnn -o ${sampleID}.cnr
|
||
$python3 $cnvkit segment -p ${threads} ${sampleID}.cnr --drop-low-coverage -o ${sampleID}.cns
|
||
fi
|
||
|
||
$python3 $cnvkit call ${sampleID}.cns --center --drop-low-coverage -o ${sampleID}.call.cns
|
||
$python3 $cnvkit export vcf ${sampleID}.call.cns -i "${sampleID}" -o ${sampleID}.cnv.vcf
|
||
$python3 $annotate_vcf -i ${sampleID}.cnv.vcf -r $refFlat -o ${sampleID}.annotated.cnv.csv
|
||
|
||
$python3 ${cna_out} ${sampleID}.annotated.cnv.csv
|
||
|
||
##################### THeTA ########################
|
||
# THeTA -- calculate sample purity (tumour/normal ratio)
|
||
$python3 $cnvkit export theta ${sampleID}.cns -r ${cnv_ref}/${bedType}.flat.cnn
|
||
${theta} ${sampleID}.interval_count --NUM_PROCESSES ${threads}
|
||
|
||
|
||
echo "CNA analysis completed at: `date` "
|
||
rm ${sampleID}.*coverage.cnn ${sampleID}.cnr ${sampleID}.cnv.vcf ${sampleID}.cns ${sampleID}.call.cns ${sampleID}.interval_count ${sampleID}.n2.graph.pdf ${sampleID}.n2.withBounds
|
||
fi
|
||
|
||
|
||
################### Microsatelite instability ######################
|
||
# ###################### MSIsensor2 ######################
|
||
# # Detect msi on below panel using msisensor
|
||
|
||
if [[ ${bedType} =~ ^(${msi_panel})$ ]];
|
||
then
|
||
echo $'\nmsisensor2'
|
||
echo "$(date)"
|
||
$msisensor2 msi -M ${msi_models} -t ${sampleID}.sorted.dedup.bam -o ${sampleID}.prefix
|
||
fi
|
||
|
||
|
||
########################## Fusion of genes ########################
|
||
##################### genefuse ########################
|
||
# Detect potential gene fusion on below panels.
|
||
if [[ ${bedType} =~ ^(${fusion_panel})$ ]]; then
|
||
|
||
echo $'\ngenefuse'
|
||
echo "$(date)"
|
||
$python3 ${make_fusion_gene} ${sampleID}
|
||
if [[ "$sampleType" =~ (ctdna) ]]; then
|
||
$genefuse -t ${threads} -r ${hg19fa} -f ${sampleID}.fusion_genes.csv -1 ${sampleID}.umi_consensus.fq.gz -h ${sampleID}.fusion.html -j ${sampleID}.fusion.json 1> ${sampleID}.fusion.txt
|
||
rm ${sampleID}.umi_consensus.fq.gz
|
||
else
|
||
$genefuse -t ${threads} -r ${hg19fa} -f ${sampleID}.fusion_genes.csv -1 ${sampleID}.R1_001.trim.fq.gz -2 ${sampleID}.R2_001.trim.fq.gz -h ${sampleID}.fusion.html -j ${sampleID}.fusion.json 1> ${sampleID}.fusion.txt
|
||
fi
|
||
|
||
echo $'\ndelly'
|
||
echo "$(date)"
|
||
$delly call -x ${delly_low_regions} -g ${hg19fa} -o ${sampleID}.bcf ${sampleID}.sorted.dedup.bam
|
||
$bcftools view ${sampleID}.bcf > ${sampleID}.sv.raw.vcf
|
||
rm ${sampleID}.bcf*
|
||
$python3 $delly_sv_filter ${sampleID}.sv.raw.vcf
|
||
rm ${sampleID}.delly*
|
||
|
||
$python3 $fusion_filter ${sampleID}.fusion.json -o ${sampleID}.fusion_filtered.json
|
||
$python3 ${fusion_out} ${sampleID}.fusion_filtered.json
|
||
$python3 ${check_fusion_result} ${sampleID}.fusion.xlsx
|
||
rm ${sampleID}.fusion.txt
|
||
fi
|
||
|
||
|
||
########################## Quality control #########################
|
||
##################### Fastp ########################
|
||
# This step is used for quality control by fastqc.
|
||
|
||
mkdir -p qc/fastp
|
||
echo $'\nQC_fastp'
|
||
echo "$(date)"
|
||
mv ${sampleID}.jion ${sampleID}.html qc/fastp/
|
||
|
||
##################### samtools ########################
|
||
# This step is used for quality control by samtools.
|
||
|
||
mkdir -p qc/samtools
|
||
echo $'\nQC_samtools'
|
||
echo "$(date)"
|
||
$samtools stats -@ ${threads} ${sampleID}.sorted.bam >qc/samtools/${sampleID}.sorted.bam.stats.txt
|
||
|
||
|
||
##################### hla_type_neo_predict ########################
|
||
|
||
# This step is used for hla_type_neo_predict by hlahd and MuPeXI.
|
||
if [[ "$bedType" =~ ^(${hla_panel})$ ]]; then
|
||
echo $'\nHLA'
|
||
echo "$(date)"
|
||
cp anno/${sampleID}.tmb.final.filter.vcf anno/${sampleID}.tmb_copy.vcf
|
||
grep '#' anno/${sampleID}.tmb_copy.vcf > anno/${sampleID}.header
|
||
grep -v '#' anno/${sampleID}.tmb_copy.vcf > anno/${sampleID}.body
|
||
|
||
gunzip -k ${sampleID}.R1_001.trim.fq.gz
|
||
gunzip -k ${sampleID}.R2_001.trim.fq.gz
|
||
# get 300w reads, about 1GB to to hla typing
|
||
$seqtk sample -s100 ${sampleID}.R1_001.trim.fq 3000000 > ${sampleID}.R1_001.trim.sub.fq
|
||
$seqtk sample -s100 ${sampleID}.R2_001.trim.fq 3000000 > ${sampleID}.R2_001.trim.sub.fq
|
||
|
||
$python2 ${hla_optitype} -sampleID $sampleID -outDir $outDir -results_TMB_txt anno/${sampleID}.body -fq1 ${sampleID}.R1_001.trim.sub.fq -fq2 ${sampleID}.R2_001.trim.sub.fq
|
||
|
||
# 拼接vcf
|
||
cat anno/${sampleID}.header anno/${sampleID}.select_variant_from_TMB.txt > anno/${sampleID}.mupexi.vcf
|
||
|
||
$python2 ${hla_optitype} -sampleID $sampleID -outDir $outDir -vcf anno/${sampleID}.mupexi.vcf -determine_tnb anno/${sampleID}.body
|
||
rm ${sampleID}.R*_001.trim.fq ${sampleID}.R*_001.trim.sub.fq \
|
||
anno/${sampleID}.body anno/${sampleID}.tmb_copy.vcf anno/${sampleID}.header anno/${sampleID}.select_variant_from_TMB.txt \
|
||
neo/${sampleID}_coverage_plot.pdf neo/${sampleID}.fasta neo/${sampleID}.g \
|
||
neo/${sampleID}_netMHCpan* neo/${sampleID}_pepmatch_* neo/${sampleID}_peptide_netMHCinput.txt neo/${sampleID}_reference_peptide_* \
|
||
neo/${sampleID}_vep_compatible_vcf.vcf neo/${sampleID}_vep.vep
|
||
fi
|
||
|
||
|
||
########################## merge_data #########################
|
||
$python3 ${single_sample_qc} . ${sampleID} ${sampleType}
|
||
$python3 ${parser_single_sample_data} . ${sampleID}
|
||
$python3 ${collect_single_sample_data} . ${sampleID}
|
||
|
||
|
||
rm ${sampleID}.R*_001.trim.fq.gz
|
||
if [[ ${bedType} =~ ^(${probe_panel})$ ]]; then
|
||
rm ${sampleID}.sorted.bam*
|
||
fi
|
||
|
||
echo "$sampleID -- ending at: $(date)"
|
||
echo "pipeline finished!" >finish.txt
|