Home Explore Help
Sign In
chaowen.xi
/
602panel_test
1
0
Fork 0
Files Issues 0 Pull Requests 0 Wiki
Branch: master
Branches Tags
602panel_test
/
s2_noUMI_sv_calling_varscan_pair_20220705.pbs

s2_noUMI_sv_calling_varscan_pair_20220705.pbs 8.7 KB
Permalink History Raw

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143 
  1. #!/bin/sh
    2. 

  2. #PBS -N noUMI_varscan_pair
    3. 

  3. #PBS -j oe
    4. 

  4. #PBS -l ncpus=12
    5. 

  5. #PBS -l nodes=1
    6. 

  6. #PBS -l mem=20G
    7. 

  7. 

  8. #inputpath=$1
    9. 

  9. #sample=$2
    10. 

  10. source  /home/liuxiangqiong/miniconda3/bin/activate base
    11. 

  11. hg19=/cgdata/Database/GATK/b37/human_g1k_v37_decoy.fasta
    12. 

  12. target=/cgdata/liuxiangqiong/work62pancancer/Client/v0/script/refdata/NanOnco_Plus_Panel_v2.0_Covered_b37_cg.parY2X.sort.bed
    13. 

  13. #min_af=0.005 for blood,0.01 for tissue
    14. 

  14. outputdir=${inputpath}/1SV_varscan_pair
    15. 

  15. #######################1.pileup file########################################################################
    16. 

  16. samtools mpileup -B -q 10 -d 1000000 -f ${hg19} -l $target ${bam_dir}/${tumor}_clean.bam > ${outputdir}/${tumor}.pileup
    17. 

  17. samtools mpileup -B -q 10 -d 1000000 -f ${hg19} -l $target ${bam_dir}/${normal}_clean.bam > ${outputdir}/${normal}.pileup
    18. 

  18. 

  19. ########################2.pair########################################################################
    20. 

  20. varscan somatic  ${outputdir}/${normal}.pileup  ${outputdir}/${tumor}.pileup --output-snp  ${outputdir}/${sample}.snp --output-indel  ${outputdir}/${sample}.indel --tumor-purity 0.01 --output-vcf
    21. 

  21. 

  22. ##对原始的结果进行合并并注释
    23. 

  23. less ${outputdir}/${sample}.snp.vcf | grep -vE "^#"  > ${outputdir}/${sample}.merge.Somatic.vcf
    24. 

  24. less ${outputdir}/${sample}.indel.vcf | grep -vE "^#" >> ${outputdir}/${sample}.merge.Somatic.vcf
    25. 

  25. #3.3 将邻近的位点进行合并#3.3 将邻近的位点进行合并
    26. 

  26. python3 /cgdata/liuxiangqiong/work62pancancer/Client/v0/script/complex_variant_new.py -i  ${outputdir}/${sample}.merge.Somatic.vcf -b  ${bam_dir}/${tumor}_clean.bam -f ${hg19} -mode 2 -n 10 -p 0.5 -o ${outputdir}/${sample}.merge.Somatic.combind1.vcf
    27. 

  27. grep "^#" ${outputdir}/${sample}.snp.vcf >${outputdir}/${sample}.merge.Somatic.combind.vcf
    28. 

  28. less  ${outputdir}/${sample}.merge.Somatic.combind1.vcf >>${outputdir}/${sample}.merge.Somatic.combind.vcf
    29. 

  29. rm -rf  ${outputdir}/${sample}.merge.Somatic.combind1.vcf
    30. 

  30. #perl /cgdata/CGTools/soft/src/annovar_20200608/table_annovar.pl ${outputdir}/${sample}.merge.Somatic.combind.vcf /cgdata/soft/src/annovar.archive/DB/ -buildver hg19 -out ${outputdir}/${sample}.merge.Somatic.combind.anno -remove -protocol refGene,popfreq_max_20150413,gnomad_genome,avsnp150,snp138NonFlagged,dbnsfp33a,cosmic94,clinvar_20210501,intervar_20170202,icgc28 -operation g,f,f,f,f,f,f,f,f,f -nastring . -vcfinput --thread 40
    31. 

  31. #5.2先拆分再标准化对齐
    32. 

  32. vt decompose ${outputdir}/${sample}.merge.Somatic.combind.vcf |vt normalize - -r ${hg19} -o ${outputdir}/${sample}.merge.Somatic.combind.norm.vcf
    33. 

  33. conda deactivate
    34. 

  34. #5.3对经过标准化的进行vep注释
    35. 

  35. source  /home/liuxiangqiong/miniconda3/bin/activate  ensemble-vep
    36. 

  36. PERL5LIB=/home/liuxiangqiong/miniconda3/envs/ensemble-vep/share/ensembl-vep-105.0-1/Bio/EnsEMBL/Variation/Utils/DB/HTS/Faidx
    37. 

  37. export PERL5LIB=$PERL5LIB
    38. 

  38. inputdir=${outputdir}/${sample}.merge.Somatic.combind.norm.vcf
    39. 

  39. outdir1=${outputdir}/${sample}.merge.Somatic.combind.norm.vepanno.vcf
    40. 

  40. vep --cache --fork 16 \
    41. 

  41. -i ${inputdir} \
    42. 

  42. -o ${outdir1} \
    43. 

  43. --assembly GRCh37 \
    44. 

  44. --cache_version 104 \
    45. 

  45. --dir_cache /home/liuxiangqiong/miniconda3/share/ensembl-vep-88.9-0/hg19 \
    46. 

  46. --refseq \
    47. 

  47. --canonical \
    48. 

  48. --offline \
    49. 

  49. --hgvs \
    50. 

  50. --fasta /cgdata/Database/GATK/b37/human_g1k_v37_decoy.fasta \
    51. 

  51. --hgvsg \
    52. 

  52. --max_af \
    53. 

  53. --check_existing \
    54. 

  54. --numbers \
    55. 

  55. --symbol \
    56. 

  56. --canonical \
    57. 

  57. --offline \
    58. 

  58. -vcf \
    59. 

  59. --use_given_ref \
    60. 

  60. --af_gnomad \
    61. 

  61. --check_existing
    62. 

  62. 

  63. #5.4对vep的再进行annovar注释
    64. 

  64. outdir2=${outputdir}/${sample}.merge.Somatic.combind.norm.vepanno.annovar
    65. 

  65. perl /cgdata/CGTools/soft/src/annovar_20200608/table_annovar.pl ${outdir1} /cgdata/soft/src/annovar.archive/DB/ -buildver hg19 -out ${outdir2} -remove -protocol refGene,popfreq_max_20150413,gnomad_genome,avsnp150,snp138NonFlagged,dbnsfp33a,cosmic94,clinvar_20210501,intervar_20170202,icgc28 -operation g,f,f,f,f,f,f,f,f,f -nastring . -vcfinput --thread 40
    66. 

  66. conda deactivate
    67. 

  67. source  /home/liuxiangqiong/miniconda3/bin/activate base
    68. 

  68. ########################3.Run the processSomatic subcommand to divide the output into separate files based on somatic status and confidence########
    69. 

  69. #3.1获得高置信的突变
    70. 

  70. varscan processSomatic ${outputdir}/${sample}.snp.vcf --min-tumor-freq 0.001 --max-normal-freq 0.05 --p-value 0.5
    71. 

  71. varscan processSomatic ${outputdir}/${sample}.indel.vcf --min-tumor-freq 0.001 --max-normal-freq 0.05 --p-value 0.5
    72. 

  72. #3.2merge the hc
    73. 

  73. less ${outputdir}/${sample}.snp.Somatic.hc.vcf | grep -vE "^#"> ${outputdir}/${sample}.merge.Somatic.hc.nohead.vcf
    74. 

  74. less ${outputdir}/${sample}.indel.Somatic.hc.vcf | grep -vE "^#" >> ${outputdir}/${sample}.merge.Somatic.hc.nohead.vcf
    75. 

  75. #3.3 将邻近的位点进行合并#3.3 将邻近的位点进行合并
    76. 

  76. python3 /cgdata/liuxiangqiong/work62pancancer/Client/v0/script/complex_variant_new.py -i  ${outputdir}/${sample}.merge.Somatic.hc.nohead.vcf -b ${bam_dir}/${tumor}_clean.bam -f ${hg19} -mode 2 -n 10 -p 0.5 -o ${outputdir}/${sample}.merge.Somatic.hc.nohead.combind.vcf
    77. 

  77. grep "^#" ${outputdir}/${sample}.snp.Somatic.hc.vcf >${outputdir}/${sample}.merge.Somatic.hc.combind.vcf
    78. 

  78. less ${outputdir}/${sample}.merge.Somatic.hc.nohead.combind.vcf >>${outputdir}/${sample}.merge.Somatic.hc.combind.vcf
    79. 

  79. #rm -rf  ${outputdir}/${sample}.merge.Somatic.hc.nohead.combind.vcf
    80. 

  80. 

  81. 

  82. #########################4.filter the false positive########################################################################
    83. 

  83. grep -v "#" ${outputdir}/${sample}.merge.Somatic.hc.combind.vcf | awk '{a=$2+10;print $1"\t"$2"\t"a}' > ${outputdir}/${sample}.merge.Somatic.hc.combind.bed
    84. 

  84. conda deactivate
    85. 

  85. source  /home/liuxiangqiong/miniconda3/bin/activate  python2
    86. 

  86. bam-readcount -q 20 -b 20 -w 100 -f ${hg19} -l ${outputdir}/${sample}.merge.Somatic.hc.combind.bed ${bam_dir}/${tumor}_clean.bam >${outputdir}/${sample}.merge.Somatic.hc.combind.readcount
    87. 

  87. conda deactivate
    88. 

  88. source  /home/liuxiangqiong/miniconda3/bin/activate base
    89. 

  89. varscan fpfilter \
    90. 

  90. ${outputdir}/${sample}.merge.Somatic.hc.combind.vcf \
    91. 

  91. ${outputdir}/${sample}.merge.Somatic.hc.combind.readcount \
    92. 

  92. --min-var-count 5 \
    93. 

  93. --min-var-basequal 20 \
    94. 

  94. --min-ref-basequal 20 \
    95. 

  95. --min-var-freq ${min_af} \
    96. 

  96. --output-file ${outputdir}/${sample}.merge.Somatic.hc.combind.pass.vcf
    97. 

  97. 

  98. ##########################5.注释########################################################################
    99. 

  99. #5.1对没有过滤的进行注释对高置信的过滤
    100. 

  100. #perl /cgdata/CGTools/soft/src/annovar_20200608/table_annovar.pl ${outputdir}/${sample}.merge.Somatic.hc.head.combind.vcf /cgdata/soft/src/annovar.archive/DB/ -buildver hg19 -out ${outputdir}/${sample}.merge.Somatic.hc.head.combind.anno -remove -protocol refGene,popfreq_max_20150413,gnomad_genome,avsnp150,snp138NonFlagged,dbnsfp33a,cosmic94,clinvar_20210501,intervar_20170202,icgc28 -operation g,f,f,f,f,f,f,f,f,f -nastring . -vcfinput --thread 40
    101. 

  101. #perl /cgdata/CGTools/soft/src/annovar_20200608/table_annovar.pl ${outputdir}/${sample}.snp.Somatic.hc.vcf /cgdata/soft/src/annovar.archive/DB/ -buildver hg19 -out ${outputdir}/${sample}.merge.Somatic.hc.anno -remove -protocol refGene,popfreq_max_20150413,gnomad_genome,avsnp150,snp138NonFlagged,dbnsfp33a,cosmic94,clinvar_20210501,intervar_20170202,icgc28 -operation g,f,f,f,f,f,f,f,f,f -nastring . -vcfinput --thread 40
    102. 

  102. 

  103. 

  104. #5.2先拆分再标准化对齐
    105. 

  105. vt decompose ${outputdir}/${sample}.merge.Somatic.hc.combind.pass.vcf |vt normalize - -r ${hg19} -o ${outputdir}/${sample}.merge.Somatic.hc.combind.pass.norm.vcf
    106. 

  106. conda deactivate
    107. 

  107. #5.3对经过标准化的进行vep注释
    108. 

  108. source  /home/liuxiangqiong/miniconda3/bin/activate  ensemble-vep
    109. 

  109. PERL5LIB=/home/liuxiangqiong/miniconda3/envs/ensemble-vep/share/ensembl-vep-105.0-1/Bio/EnsEMBL/Variation/Utils/DB/HTS/Faidx
    110. 

  110. export PERL5LIB=$PERL5LIB
    111. 

  111. inputdir0=${outputdir}/${sample}.merge.Somatic.hc.combind.pass.norm.vcf
    112. 

  112. outdir1_0=${outputdir}/${sample}.merge.Somatic.hc.combind.pass.norm.vepanno.vcf
    113. 

  113. vep --cache --fork 16 \
    114. 

  114. -i ${inputdir0} \
    115. 

  115. -o ${outdir1_0} \
    116. 

  116. --assembly GRCh37 \
    117. 

  117. --cache_version 104 \
    118. 

  118. --dir_cache /home/liuxiangqiong/miniconda3/share/ensembl-vep-88.9-0/hg19 \
    119. 

  119. --refseq \
    120. 

  120. --canonical \
    121. 

  121. --offline \
    122. 

  122. --hgvs \
    123. 

  123. --fasta /cgdata/Database/GATK/b37/human_g1k_v37_decoy.fasta \
    124. 

  124. --hgvsg \
    125. 

  125. --max_af \
    126. 

  126. --check_existing \
    127. 

  127. --numbers \
    128. 

  128. --symbol \
    129. 

  129. --canonical \
    130. 

  130. --offline \
    131. 

  131. -vcf \
    132. 

  132. --use_given_ref \
    133. 

  133. --af_gnomad \
    134. 

  134. --check_existing
    135. 

  135. 

  136. #5.4对vep的再进行annovar注释
    137. 

  137. outdir2_0=${outputdir}/${sample}.merge.Somatic.hc.combind.pass.norm.vepanno.annovar
    138. 

  138. perl /cgdata/CGTools/soft/src/annovar_20200608/table_annovar.pl ${outdir1_0} /cgdata/soft/src/annovar.archive/DB/ -buildver hg19 -out ${outdir2_0} -remove -protocol refGene,popfreq_max_20150413,gnomad_genome,avsnp150,snp138NonFlagged,dbnsfp33a,cosmic94,clinvar_20210501,intervar_20170202,icgc28 -operation g,f,f,f,f,f,f,f,f,f -nastring . -vcfinput --thread 40
    139. 

  139. 

  140. ###############################6.删除中间文件
    141. 

  141. #rm -rf ${outputdir}/${normal}.pileup 
    142. 

  142. #rm -rf ${outputdir}/${tumor}.pileup
    143. 

  143. rm -rf ${outputdir}/${sample}.merge.Somatic.hc.readcount
    144.