GATK4-germline-mut-study Day-1

作者: liu_ll | 来源:发表于2018-12-11 23:22 被阅读35次

GATK4-germline-mut-study Day-1
头部效应——站位比努力更重要
尽管我们的手中空无一物
User Experience
Day1---Dr.leng
尼泊尔之旅-加德满都
python学习目录
Day-1
Day-1
Day-1

这几天一直都在忙着规划自己的学习安排，GATK的学习也是间断的，中间还搁置了一段的时间。我发现，不能推进自己的学习计划简直是太可怕了。。。。。
用代码记录一下自己学习的进度

#!/bin/sh
#PBS -q middleq
#PBS -l mem=40gb,walltime=168:00:00,nodes=1:ppn=1
#PBS -N wgs_gatk4
#HSCHED -s ctDNA+wgs-gatk4+Human
# Author: Liu linli

#this pipeline is to preprocess the WGS data befote varients calling 

#prepare reference genome
#all data in /share_bio/unisvx1/sunyl_group/liull/Database/bundlehg38/
#clean data in /share_bio/unisvx1/sunyl_group/liull/twins_WGS/BCA0106/2_bwa  /share_bio/unisvx1/sunyl_group/liull/twins_WGS/BCA0106-2

genome=/share_bio/unisvx1/sunyl_group/liull/Database/bundlehg38/Homo_sapiens_assembly38.fasta
1000G_ref=/share_bio/unisvx1/sunyl_group/liull/Database/bundlehg38/1000G_phase1.snps.high_confidence.hg38.vcf 
Mill_1000G=/share_bio/unisvx1/sunyl_group/liull/Database/bundlehg38/Mills_and_1000G_gold_standard.indels.hg38.vcf 
Hapmap=/share_bio/unisvx1/sunyl_group/liull/Database/bundlehg38/hapmap_3.3.hg38.vcf 
Ommin=/share_bio/unisvx1/sunyl_group/liull/Database/bundlehg38/1000G_omni2.5.hg38.vcf 
dbsnp=/share_bio/unisvx1/sunyl_group/liull/Database/bundlehg38/dbsnp_146.hg38.vcf 
sample1=/share_bio/unisvx1/sunyl_group/liull/twins_WGS/BCA0106/2_bwa/BCA0106_R1.fq.gz 
sample2=/share_bio/unisvx1/sunyl_group/liull/twins_WGS/BCA0106/2_bwa/BCA0106_R2.fq.gz
sample3=/share_bio/unisvx1/sunyl_group/liull/twins_WGS/BCA0106-2/BCA0106-2_R1.fq.gz
sample4=/share_bio/unisvx1/sunyl_group/liull/twins_WGS/BCA0106-2/BCA0106-2_R2.fq.gz 
outdir=/share_bio/unisvx1//sunyl_group/liull/twins_WGS/outdir
#bwa mapping to reference 
bwa mem -t 8 -M -R "@RG\tID:L1\tSM:BCA0106\tLB:WGS\tPL:Illumina" $genome $sample1 $sample2 |samtools view -Sb - > $outdir/bwa/BCA0106.paired.bam 
bwa mem -t 8 -M -R  "@RG\tID:L1\tSM:BCA010-2\tLB:WGS\tPL:Illumina" $genome $sample3 $sample4 |samtools view -Sb - > $outdir/bwa/BCA0106-2.paired.bam 
#samtools sort bam 
samtools sort -@ 6 -m 6G  -o $outdir/BCA0106.paired.sorted.bam $outdir/BCA0106.paired.bam 
samtools sort -@ 6 -m 6G  -o $outdir/BCA0106.paired-2.sorted.bam $outdir/BCA0106-2.paired.bam 

#markduplication
gatk --java-options"-Xmx10G" MarkDuplicates \
    -I $outdir/BCA0106.paired.sorted.bam \
    -O $outdir/BCA0106.paired.sorted.MarkDuplicates.bam \
    -M $outdir/BCA0106.sorted.bam.metrics \
    --CREAT_INDEX
gatk --java-options"-Xmx10G" MarkDuplicates \
    -I $outdir/BCA0106-2.paired.sorted.bam \
    -O $outdir/BCA0106-2.paired.sorted.MarkDuplicates.bam \
    -M $outdir/BCA0106-2.sorted.bam.metrics \
    --CREAT_INDEX
#GATK BaseRecalibrator 
gatk BaseRecalibrator \
    -R $genome \
    -I  $outdir/BCA0106.paired.sorted.MarkDuplicates.bam \
    -O $outdir/BCA0106.paired.sorted.MarkDuplicates.recal.table \
   --known-sites  $Mill_1000G \
   --known-sites $dbsnp \
   --konwn-sites $Mill_1000G 
gatk BaseRecalibrator \
    -R $genome \
    -I  $outdir/BCA0106-2.paired.sorted.MarkDuplicates.bam \
    -O $outdir/BCA0106-2.paired.sorted.MarkDuplicates.recal.table \
   --known-sites  $Mill_1000G \
   --known-sites $dbsnp \
   --konwn-sites $Mill_1000G 
# GATK ApplyBQSR
gatk ApplyBQSR \
    -R $genome \
    -I $outdir/BCA0106.paired.sorted.MarkDuplicates.bam \
    -bqsr $outdir/BCA0106.paired.sorted.MarkDuplicates.recal.table \
    -O $outdir/BCA0106.paired.sorted.MarkDuplicates.BQSR.bam
gatk ApplyBQSR \
    -R $genome \
    -I $outdir/BCA0106-2.paired.sorted.MarkDuplicates.bam \
    -bqsr $outdir/BCA0106-2.paired.sorted.MarkDuplicates.recal.table \
    -O $outdir/BCA0106-2.paired.sorted.MarkDuplicates.BQSR.bam

# samtools Bam index 
samtools index $outdir/BCA0106.paired.sorted.MarkDuplicates.BQSR.bam
samtools index $outdir/BCA0106-2.paired.sorted.MarkDuplicates.BQSR.bam

#gatk HaplotypeCaller
gatk HaplotypeCaller \
    -R $genome \
    -I $outdir/BCA0106.paired.sorted.MarkDuplicates.BQSR.bam \
    -O $outdir/BCA0106.paired.sorted.MarkDuplicates.BQSR.vcf \
    -D $dbsnp 
gatk HaplotypeCaller \
    -R $genome \
    -I $outdir/BCA0106-2.paired.sorted.MarkDuplicates.BQSR.bam \
    -O $outdir/BCA0106-2.paired.sorted.MarkDuplicates.BQSR.vcf \
    -D $dbsnp
#对snp和indel分别进行分析，首先是snp的模式
gatk VariantRecalibrator \
    -R $genome \
    -V $outdir/BCA0106.paired.sorted.MarkDuplicates.BQSR.vcf \
    -resource hapmap,known=false,training=true,truth=true,prior=15.0:$Hapmap \
    -resource omini,known=false,training=true,truth=false,prior=12.0:$Ommin \
    -resource 1000G,known=false,training=true,truth=false,prior=10.0:$1000G_ref \
    -resource dbsnp,known=true,training=false,truth=false,prior=6.0:$dbsnp \
    -an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR -an DP \
    -mode SNP \
    -O $outdir/BCA0106.GATK.snps.recal.vcf \
    --rscript-file $outdir/wes.GATK.snps.plots.R \
    --tranches-file $outdir/BCA0106.gatk.snps.recal.tranches 
gatk  ApplyRecalibration  \
    -R $genome 
    -V $outdir/BCA0106.GATK.snps.recal.vcf \
    -O BCA0106.GATK_VQSR.vcf \
    --recal-file wes.GATK.snps.recal.vcf  --tranches-file wes.GATK.snps.tranches  -mode SNP

到这一步有点卡克了，没太明白是什么意思。希望DAY2的学习可以进步一点点，这样即使是蜗牛爬也不怕了