天真的我准备把全部流程迁移到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
纯干货代码，谁学谁获益。
■ ■ ■
生信基础知识大全系列：生信基础知识100讲
史上最强的生信自学环境准备课来啦！！ 7次改版，11节课程，14K的讲稿，30个夜晚打磨，100页PPT的课程。
如果需要组装自己的服务器；代办生物信息学服务器
如果需要帮忙下载海外数据(GEO/TCGA/GTEx等等)，点我？
如果需要线下辅导及培训，看招学徒
如果需要个人电脑：个人计算机推荐
如果需要置办生物信息学书籍，看：生信人必备书单
如果需要实习岗位：实习职位发布
如果需要售后：点我
如果需要入门资料大全：点我