天真的我准备把全部流程迁移到GATK4

我在生信技能树上面发布的GATK4教程也有不少了 本着尽量使用最新版软件的原则,也准备把之前的gatk对RNA-seq数据找变异的流程进行转换:

$GATK  --java-options "-Xmx25G -Djava.io.tmpdir=./" AddOrReplaceReadGroups \
    -I $id -O ${sample}_right.bam -SO coordinate -ID ${sample}  -LB rna \
    -PL illumina -PU hiseq  -SM ${sample}

    $GATK  --java-options "-Xmx25G -Djava.io.tmpdir=./"  MarkDuplicates \
    -I ${sample}_right.bam -O ${sample}_marked.bam -M $sample.metrics --REMOVE_DUPLICATES TRUE

    $GATK  --java-options "-Xmx25G -Djava.io.tmpdir=./"  FixMateInformation \
    -I ${sample}_marked.bam -O ${sample}_marked_fixed.bam  -SO coordinate 

    $GATK  --java-options "-Xmx25G -Djava.io.tmpdir=./"  SplitNCigarReads \
    -R $GENOME  -I ${sample}_marked_fixed.bam  -O ${sample}_marked_fixed_split.bam \
    -rf ReassignOneMappingQuality -RMQF 255 -RMQT 60 -U ALLOW_N_CIGAR_READS

    #--fix_misencoded_quality_scores
    ## --fix_misencoded_quality_scores only if phred 64

但是走到了 SplitNCigarReads 才发现,这个命令当初学的太久了,忘记各个参数啥意思了,就想搜索看看如何转换。

还真发现了有人问同样的问题,GATK4: How to reassign STAR mapping quality from 255 to 60 with SplitNCigarReads ,而且GATK4开发团队也回答了:EDIT: Geraldine responded here.

但是这是一个否定回答,开发团队让我们回去用GATK3来跑流程。

One risk that I see is that using the STAR --outSAMmapqUnique 60 option maybe fixes the issue with GATK, but that other downstream tools maybe still depend on the (still default) STAR mapping quality value of 255 (e.g. cufflinks).

The mapping quality MAPQ (column 5) is 255 for uniquely mapping reads, and int(-10*log10(1- 1/Nmap)) for multi-mapping reads. This scheme is same as the one used by TopHat and is com- patible with Cuffinks. The default MAPQ=255 for the unique mappers maybe changed with --outSAMmapqUnique parameter (integer 0 to 255) to ensure compatibility with downstream tools such as GATK.

https://github.com/alexdobin/STAR/blob/master/doc/STARmanual.pdf

好吧,回去就回去,gatk3代码是:

module load java/1.8.0_91 
GATK=/home/jianmingzeng/biosoft/GATK/GenomeAnalysisTK.jar
PICARD=/home/jianmingzeng/biosoft/picardtools/2.9.2/picard.jar 
GENOME=/home/jianmingzeng/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.fasta
INDEX=/home/jianmingzeng/biosoft/GATK/resources/bundle/hg38/bwa_index/gatk_hg38 
DBSNP=/home/jianmingzeng/biosoft/GATK/resources/bundle/hg38/dbsnp_146.hg38.vcf.gz
kgSNP=/home/jianmingzeng/biosoft/GATK/resources/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz
kgINDEL=/home/jianmingzeng/biosoft/GATK/resources/bundle/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz

TMPDIR=/home/jianmingzeng/tmp/software
## samtools and bwa are in the environment
## samtools Version: 1.3.1 (using htslib 1.3.1)
## bwa Version: 0.7.15-r1140
cat $1 |while read id
do
    echo $id
    file=$(basename $id )
    sample=${file%%_*}
    echo $sample
    java -Djava.io.tmpdir=$TMPDIR    -Xmx25g -jar $PICARD  AddOrReplaceReadGroups \
    I=$id O=${sample}.bam SO=coordinate RGID=${sample}  RGLB=rna \
    RGPL=illumina RGPU=hiseq  RGSM=${sample}
    java -Djava.io.tmpdir=$TMPDIR    -Xmx25g -jar $PICARD MarkDuplicates \
    INPUT=${sample}.bam OUTPUT=${sample}_marked.bam METRICS_FILE=$sample.metrics REMOVE_DUPLICATES=TRUE
    java -Djava.io.tmpdir=$TMPDIR    -Xmx25g -jar $PICARD FixMateInformation \
    INPUT=${sample}_marked.bam OUTPUT=${sample}_marked_fixed.bam SO=coordinate

    samtools index ${sample}_marked_fixed.bam
    java -Djava.io.tmpdir=$TMPDIR   -Xmx25g -jar $GATK -T SplitNCigarReads \
    -R $GENOME  -I ${sample}_marked_fixed.bam  -o ${sample}_marked_fixed_split.bam \
    -rf ReassignOneMappingQuality -RMQF 255 -RMQT 60 -U ALLOW_N_CIGAR_READS
    #--fix_misencoded_quality_scores
    ## --fix_misencoded_quality_scores only if phred 64
    java -Djava.io.tmpdir=$TMPDIR   -Xmx25g -jar $GATK -T HaplotypeCaller  \
    -R $GENOME -I ${sample}_marked_fixed_split.bam --dbsnp $DBSNP  \
    -stand_emit_conf 10 -o  ${sample}_raw.vcf
    rm ${sample}.bam ${sample}_marked.bam ${sample}_marked_fixed.bam ${sample}_marked_fixed_split.bam
done

纯干货代码,谁学谁获益。

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2018-11-23

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

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏HansBug's Lab

3298: [USACO 2011Open]cow checkers

3298: [USACO 2011Open]cow checkers Time Limit: 10 Sec  Memory Limit: 128 MB Subm...

28360
来自专栏CreateAMind

autoware 代码阅读 sensor

ROS driver to parse NMEA strings and publish standard ROS NavSat message types. ...

24920
来自专栏阿杜的世界

【小马哥】Spring Boot系列讲座系列套餐讲座大纲

20710
来自专栏杨建荣的学习笔记

hotspare的copyback(r7笔记第30天)

最近做硬件巡检,发现一部分硬盘出现了坏块,同事就帮忙去协调处理这个事情,晚些时候接到了现场工程师的电话,问我可以不可以换,简单确认是raid5的盘。所以只能一个...

43350
来自专栏耕耘实录

国外程序员推荐的免费编程书籍资源

我正试着搜集整理一个可在网上免费阅读的计算机编程书籍列表。这些书可以是某种特定编程语言,也可以计算机方面通用书籍。网上有哪些免费可用的书籍呢?

43540
来自专栏HansBug's Lab

1615: [Usaco2008 Mar]The Loathesome Hay Baler麻烦的干草打包机

1615: [Usaco2008 Mar]The Loathesome Hay Baler麻烦的干草打包机 Time Limit: 5 Sec  Memory ...

342110
来自专栏HansBug's Lab

1610: [Usaco2008 Feb]Line连线游戏

1610: [Usaco2008 Feb]Line连线游戏 Time Limit: 5 Sec  Memory Limit: 64 MB Submit: 13...

33460
来自专栏ml

HDUOJ-----取(m堆)石子游戏

取(m堆)石子游戏 Time Limit : 3000/1000ms (Java/Other)   Memory Limit : 32768/32768K (J...

37170
来自专栏HansBug's Lab

1578: [Usaco2009 Feb]Stock Market 股票市场

1578: [Usaco2009 Feb]Stock Market 股票市场 Time Limit: 10 Sec  Memory Limit: 64 MB S...

276110
来自专栏张高兴的博客

张高兴的 Windows 10 IoT 开发笔记:FM 电台模块 KT0803L

39460

扫码关注云+社区

领取腾讯云代金券