专栏首页liu_ll的生信学习笔记GATK4-germline-mut-study Day-1

抱歉,你查看的文章已删除

GATK4-germline-mut-study 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的学习可以进步一点点,这样即使是蜗牛爬也不怕了

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • Matlab 飞机航向INS仿真

    close all; clear all; clc; format long %初始化参数 %vx(1)=0.000048637; %vy(1)=0.00020...

    Pulsar-V
  • .Net Core微服务入门全纪录(六)——EventBus-事件总线

    上一篇【.Net Core微服务入门全纪录(五)——Ocelot-API网关(下)】中已经完成了Ocelot + Consul的搭建,这一篇简单说一下Event...

    xhznl
  • 一线大厂十年架构师总结整理出的微服务详解「PDF」

    服务注册就是维护一个登记簿,它管理系统内所有的服务地址。当新的服务启动后,它会向登记 簿交待自己的地址信息。服务的依赖方直接向登记簿要Service Provi...

    Java知音
  • 垂直电商架构进化之路

    作者:张增、邓良驹,分别为乐视云计算电商云团队负责人,乐视云计算高级开发工程师 来自:高效运维 1. 电商系统发展过程 电商网站在不同时期的架构复杂度有所不同:...

    架构师小秘圈
  • springboot使用hibernate validator校验

    一、参数校验  在开发中经常需要写一些字段校验的代码,比如字段非空,字段长度限制,邮箱格式验证等等,写这些与业务逻辑关系不大的代码个人感觉有两个麻烦: 验证代码...

    庞小明
  • T4随记

    MSDN的介绍 https://msdn.microsoft.com/zh-cn/library/bb126478.aspx

    冰封一夏
  • CUDA PTX ISA阅读笔记(一)

    不知道这是个啥的看这里:Parallel Thread Execution ISA Version 5.0. 简要来说,PTX就是.cu代码编译出来的一种东西...

    用户1148523

扫码关注云+社区

领取腾讯云代金券