第一步 将测序原始的序列除去宿主的基因组信息需要软件是minimap2，samtools, seqtk
minimap2给宿主基因组fasta建立index
minimap2 -I 80g -t 10 -d Host.ref.mmi Pangenome_ALLmerge.fa
minimap2 除去宿主序列信息生成过滤后的序列，并压缩
minimap2 -t 56 -N 50 -ax sr Host.ref.mmi xx.fq1.gz xx.fq2.gz |samtools view | awk '{if ($12~/NM:i:/) {split($12,a,":"); if(a[3]>=15) print $1 } else print $1 }' >xx.lis; sort xx.lis | uniq >xx.nr.lis; seqtk subseq xx.fq1.gz xx.nr.lis >xx.filtered.fq1; seqtk subseq xx.fq2.gz xx.nr.lis >xx.filtered.fq2; gzip xx.filtered.fq1; gzip xx.filtered.fq2
第二步 gatk FastqToSam命令 带着测序的平台等信息，将FASTQ格式的测序数据转换为未对齐的BAM文件
gatk FastqToSam -F1 xx.filtered.fq1.gz -F2 xx.filtered.fq2.gz -O xx.unaligned.bam -SM xx -RG rg0013 -PL ILLUMINA
第三步 对细菌的参考基因组建立index
samtools faidx ref.fa
gatk BwaMemIndexImageCreator -I ref.fa -O ref.fa.img
gatk4 CreateSequenceDictionary ref.fa
第四步 BwaAndMarkDuplicatesPipelineSpark 把未比对的unaligned.BAM 文件，在本地 Spark 环境下，并行完成 BWA-MEM 比对 + MarkDuplicates 去重，生成可直接用于 BQSR/变异检测的干净 BAM
gatk --java-options "-Xmx20G" BwaAndMarkDuplicatesPipelineSpark --spark-master local[4] --conf "spark.executor.memory=8g" --conf "spark.driver.memory=10g" -R ../04.snpcalling/dr_mox_cs.fa -I xx.unaligned.bam -O xx.derep.bam
第五步 SortSamSpark
gatk --java-options "-Xmx50G -XX:+UseParallelGC -XX:ParallelGCThreads=2" SortSamSpark --spark-runner LOCAL --spark-master local[*] --conf spark.executor.cores=6 --conf spark.executor.memory=90g --conf spark.num.executors=10 -I xx.derep.bam -O xx.sorted.bam
第六步 RealignerTargetCreator + IndelRealigner
gatk3 -T RealignerTargetCreator -R ../04.snpcalling/dr_mox_cs.fa -I xx.sorted.bam --minReadsAtLocus 4 -o xx.intervals
gatk3 -T IndelRealigner -R ../04.snpcalling/dr_mox_cs.fa -I xx.sorted.bam -targetIntervals xx.intervals -LOD 3.0 -o xx.realigned.bam
第七步 计算样品的深度 和所有样品的平均深度
samtools depth xx.realigned.bam | awk '{sum+=$3;cnt++}END{print sum/cnt}' > xx.depth.txt
scripts/cov_mean.py
